text:
# Generate static graphs
import os
import sys
import json
import csv
import plotly.graph_objects as go
import plotly.express as px
import numpy as np

from plotly.subplots import make_subplots
from datetime import *
from utils import *

# Generate a graph (cases) for Maryland Zip Code
# Here we pass all the data
def generate_MD_zip_graph_with_avg(data,name,folder,_color):

   if(len(data['stats'])>0):
      all_x = []
      all_y = []
      all_x_avg = []
      all_y_avg = []

      # 7 days average
      first_val = -1
      total_day = 0
      max_day = 7 # Average based on max_day days

      _type = "cases"

      tempValForAvg = []
      tempValFormax_day = []

      for d in data['stats']:
         for day in d:

            # Warning - it looks like in Maryland, they put a lot of cases on the first day in the stats
            # so we ignore the data in the graphs to have something legible
            d_day = day.split('-')
            b1 = date(int(d_day[0]), int(d_day[1]), int(d_day[2]))

            if(b1 > date(2020, 4, 12)):

               # Org Data
               all_x.append(day)
               all_y.append(d[day][_type])

               # For average of _type
               tempValForAvg.append(float(d[day][_type]))

               if(len(tempValForAvg) <  max_day):
                  tempValFormax_day = tempValForAvg
               else:
                  tempValFormax_day = tempValForAvg[len(tempValForAvg)-max_day:len(tempValForAvg)]

               # We have strings...
               tempValFormax_day = [float(i) for i in tempValFormax_day]

               all_x_avg.append(day)
               all_y_avg.append(np.mean(tempValFormax_day))

      if(_color=="r"):
         _color = "red"
      elif(_color=="g"):
         _color = "green"
      elif(_color=="o"):
         _color = "orange"
      else:
         _color = "black"

      print("Generating graph for zip: " + name + " color: " + _color)

      fig = go.Figure()
      fig.add_trace(go.Bar(x=all_x, y=all_y, marker_color='rgba(158,158,158,.4)' ))
      fig.add_trace(go.Scatter(x=all_x_avg, y=all_y_avg, marker_color=_color))

      # Add line to every 1s & 15th of all months
      for d in all_x:
         if(d.endswith('15') or d.endswith('01')):
            fig.add_shape(
               type="line",
               x0=d,
               y0=0,
               x1=d,
               y1=np.max(all_y),
               opacity=0.4,
               line=dict(
                  color="rgba(0,0,0,.7)",
                  width=1,
                  dash="dot",
               )
            )

      fig.update_xaxes(rangemode="nonnegative")
      fig.update_yaxes(rangemode="nonnegative")

      fig.update_layout(
         width=350,
         height=350,
         margin=dict(l=30, r=20, t=0, b=20),   # Top 0 with no title
         paper_bgcolor='rgba(255,255,255,1)',
         plot_bgcolor='rgba(255,255,255,1)',
         showlegend= False
      )

      #print(folder + name + " > created")
      fig.write_image(folder + name + ".png")

# Generate Large Graph for the State Detail page
# with 3day average line value, new cases & tests
# WARNING THIS IS THE DUAL AXIS VERSION WITH CASES & TESTS
def generate_dual_graph_X_and_cases(state, _color, folder, dataset1, dataset2, output_ext, large = False):

   # Daily Data Source File
   cur_json_file = open(PATH_TO_STATES_FOLDER + os.sep + state + os.sep + state + ".json", 'r')
   data = json.load(cur_json_file)

   # Structure for datasets1
   all_x1 = []
   all_y1 = []

   all_x2=[]
   all_y2=[]

   # Get raw data
   for d in data['stats']:

      for day in d:

         # Dataset1
         all_x1.append(day)
         all_y1.append(d[day][dataset1['_type']])

         # Dataset2
         all_x2.append(day)
         all_y2.append(d[day][dataset2['_type']])

   if(_color=="r"):
      _color = "red"
      _3dcolor = "rgba(255,0,0,0.5)"
   elif(_color=="g"):
      _color = "green"
      _3dcolor = "rgba(34,139,34,0.5)"
   elif(_color=="o"):
      _color = "orange"
      _3dcolor = "rgba(255,165,0,0.5)"
   else:
      _color = "black"
      _3dcolor = "rgba(0,0,0,0.5)"

    # Create Fig with secondary ax
   fig = make_subplots(specs=[[{"secondary_y": True}]])

   fig.update_xaxes(rangemode="nonnegative")
   fig.update_yaxes(rangemode="nonnegative")

   # Add the bars (always dataset2)
   fig.add_trace(go.Bar(x=all_x2, y=all_y2, marker_color='rgba(158,158,158,.4)', name= dataset2['name'] ),  secondary_y=True )

   # AVGs in dataset1
   maxY_avg_dataset1 = [0]
   if('avg' in dataset1):
      line_width = len(dataset1['avg'])
      if(line_width==1):
         line_width = 2
      for avg in dataset1['avg']:
         avgx, avgy, delta =  get_avg_data(avg,state,dataset1['_type'])
         fig.add_trace(go.Scatter(x=avgx, y=avgy, name= str(avg) +"-Day Avg " + dataset1['name'],  line=dict(color="purple",width=line_width)))
         maxY_avg_dataset1.append(np.max(avgy))
         if('raw_line' in dataset1):
            if(dataset1['raw_line'] is False):
               line_width -= 1

   # AVGs in dataset2
   maxY_avg_dataset2 = [0]
   if('avg' in dataset2):
      line_width = len(dataset2['avg'])
      if(line_width==1):
         line_width = 2
      for avg in dataset2['avg']:
         avgx, avgy, delta =  get_avg_data(avg,state,dataset2['_type'])
         fig.add_trace(go.Scatter(x=avgx, y=avgy, name= str(avg) +"-Day Avg " + dataset2['name'],  line=dict(color= _color,width=line_width)),secondary_y=True)
         maxY_avg_dataset2.append(np.max(avgy))
         if('raw_line' in dataset2):
            if(dataset2['raw_line'] is False):
               line_width -= 1

   # Dataset1
   if('raw_line' in dataset1):
      if(dataset1['raw_line'] is True):
         fig.add_trace(go.Scatter(x=all_x1, y=all_y1, name= dataset1['name'],  line=dict(color="purple",width=2)))
         maxY_avg_dataset1 = [np.max(all_y1)]
   else:
      maxY_avg_dataset1 = [0]

   # Dataset2
   if('raw_line' in dataset2):
      if(dataset2['raw_line'] is True):
         fig.add_trace(go.Scatter(x=all_x2, y=all_y2, name= dataset2['name'],  line=dict(color=_color,width=2)),secondary_y=True)
         maxY_avg_dataset2 = [np.max(all_y2)]
   else:
      maxY_avg_dataset2 = [0]

   # Get MAX Y for drawing the 1st & 15th lines
   # + the lockdown period
   if(len(maxY_avg_dataset2)==1 and len(maxY_avg_dataset1)==1):
      max_y = np.max(all_y1)
   elif(len(maxY_avg_dataset2)==1 and len(maxY_avg_dataset1)!=1):
      max_y = np.max(maxY_avg_dataset1)
   elif(len(maxY_avg_dataset1)==1 and len(maxY_avg_dataset2)!=1):
      max_y = np.max(all_y1)
   else:
      max_y = np.max(maxY_avg_dataset1)

   # Add line to every 1s & 15th of all months
   for date in all_x1:
      if(date.endswith('15') or date.endswith('01')):
         fig.add_shape(
            type="line",
            x0=date,
            y0=0,
            x1=date,
            y1=max_y,
            opacity=0.4,
            line=dict(
               color="rgba(0,0,0,.5)",
               width=1,
            ),
            layer="below"
         )

   # Do we have a period in key-dates.txt
   # for the current state?
   key_dates = open(KEY_DATES,'r')
   csv_reader = csv.DictReader(key_dates)
   rows = list(csv_reader)
   start_lockdown_date  = -1
   end_lockdown_date    = -1

   for key_date_row in rows:
      if(key_date_row['state']==state):
         start_lockdown_date = key_date_row['start']

         if(key_date_row['end'] is not None):
            end_lockdown_date = key_date_row['end']

   # Add lockdown period
   if(start_lockdown_date!=-1 and start_lockdown_date is not None ):

      if(end_lockdown_date!=-1 and end_lockdown_date is not None):

         fig.add_shape(
            type="rect",
            x0=start_lockdown_date,
            y0=0,
            x1=end_lockdown_date,
            y1=max_y,
            fillcolor="LightSalmon",
            opacity=0.1,
            line_width=0,
            layer="below"
         )

      elif(start_lockdown_date is not None):

         fig.add_shape(
            type="rect",
            x0=start_lockdown_date,
            y0=0,
            x1=all_x1[len(all_x1)-1],
            y1=max_y,
            fillcolor="LightSalmon",
            opacity=0.1,
            line_width=0,
            layer="below"
         )

   if(large is True):

      fig.update_layout(
         width=1000,
         height=450,
         title = US_STATES[state] + " Tests and New Cases",
         margin=dict(l=30, r=20, t=45, b=30),   # Top Title
         paper_bgcolor='rgba(255,255,255,1)',
         plot_bgcolor='rgba(255,255,255,1)',
         showlegend= True,
         yaxis2=dict(  showgrid=False,
             titlefont=dict(
               color=_color
            )
         ),
          yaxis1=dict(
             showgrid=False,
             titlefont=dict(
               color="purple"
            )
         ),
         legend_orientation="h"
      )

   else:

      fig.update_layout(
         width=455,
         height=290,
         margin=dict(l=30, r=20, t=0, b=20),   # Top 0 with no title
         paper_bgcolor='rgba(255,255,255,1)',
         plot_bgcolor='rgba(255,255,255,1)',
         showlegend= False,
         yaxis2=dict(
            showgrid=False,
            titlefont=dict(   color=_color   )
         ),
         yaxis1=dict(
            showgrid=False,
            titlefont=dict(
               color="purple"
            )
         ),
         legend_orientation="h"
      )

   fig.update_yaxes(title_text=""+dataset1['name']+"", secondary_y=False)
   fig.update_yaxes(title_text=""+dataset2['name']+"", secondary_y=True)

   # Create file
   fig.write_image(folder + os.sep + state + output_ext)
   print(folder + os.sep + state + output_ext + " Created")

# Generate a graph based on state, type (like deaths, cases, mortality etc.) & color
# For states & county
def generate_graph_with_avg(state, _type, _color, folder, county, large=False):

   # Get JSON Data for current state or county
   if(county != '' and 'for_a_state' not in county):
      cur_json_file = open(PATH_TO_STATES_FOLDER + os.sep + state + os.sep + "counties" + os.sep +  county + ".json", 'r')
   else:
      cur_json_file = open(PATH_TO_STATES_FOLDER + os.sep + state + os.sep + state + ".json", 'r')

   data = json.load(cur_json_file)

   # Do we have a period in key-dates.txt
   # for the current state?
   key_dates = open(KEY_DATES,'r')
   csv_reader = csv.DictReader(key_dates)
   rows = list(csv_reader)
   start_lockdown_date  = -1
   end_lockdown_date    = -1

   for key_date_row in rows:
      if(key_date_row['state']==state):
         start_lockdown_date = key_date_row['start']

         if(key_date_row['end'] is not None):
            end_lockdown_date = key_date_row['end']

   if(county=="" or 'for_a_state' in county):
      all_data = data['stats']
   else:
      all_data = data['stats']
      # We sort the data by inverse date for counties
      all_data = list(reversed(all_data))

   # All Data
   all_x = []
   all_y = []

   # Special to get % on the graphs
   if(_type == 'test_pos_p'):

      for d in all_data:
         for day in d:
            # Org Data
            all_x.append(day)
            all_y.append(d[day][_type]/100) # To get the % in the graphs

   elif(_type != 'mortality'):
      # 7day Average Data
      all_x_avg, all_y_avg, delta = get_X_day_avg(7,all_data,_type)

      for d in all_data:
         for day in d:
            # Org Data
            all_x.append(day)
            all_y.append(d[day][_type])

   else:

      # We compute the Mortality rate
      # (total_dead/total_cases) *100
      for d in all_data:
         for day in d:
            all_x.append(day)
            if(d[day]['total_c']>0):
               all_y.append(d[day]['total_d']/d[day]['total_c'])
            else:
               all_y.append(0)

   # Compute the 7d avg for mortality
   tempValForAvg = []
   tempValFormax_day = []
   all_x_avg = []
   all_y_avg = []
   max_day = 7

   for c,y in enumerate(all_y):

      # For average of _type
      tempValForAvg.append(y)

      if(len(tempValForAvg) <  max_day):
         tempValFormax_day = tempValForAvg
      else:
         tempValFormax_day = tempValForAvg[len(tempValForAvg)-max_day:len(tempValForAvg)]

      all_x_avg.append(all_x[c])
      all_y_avg.append(np.mean(tempValFormax_day))

   if(_color=="r"):
      _color = "red"
   elif(_color=="g"):
      _color = "green"
   elif(_color=="o"):
      _color = "orange"
   else:
      _color = "black"

   fig = go.Figure()
   fig.add_trace(go.Bar(x=all_x, y=all_y, marker_color='rgba(158,158,158,.4)'))
   fig.add_trace(go.Scatter(x=all_x_avg, y=all_y_avg, marker_color=_color))

   # Add line to every 1s & 15th of all months
   for date in all_x:
      if(date.endswith('15') or date.endswith('01')):
         fig.add_shape(
            type="line",
            x0=date,
            y0=0,
            x1=date,
            y1=np.max(all_y),
            opacity=0.4,
            line=dict(
                color="rgba(0,0,0,.5)",
                width=1,
            )
         )

   # Add lockdown period
   if(start_lockdown_date!=-1 and start_lockdown_date is not None ):

      if(end_lockdown_date!=-1 and end_lockdown_date is not None):

         fig.add_shape(
            type="rect",
            x0=start_lockdown_date,
            y0=0,
            x1=end_lockdown_date,
            y1=np.max(all_y),
            fillcolor="LightSalmon",
            opacity=0.1,
            layer="below",
            line_width=0,
         )

      elif(start_lockdown_date is not None):

         fig.add_shape(
            type="rect",
            x0=start_lockdown_date,
            y0=0,
            x1=all_x[len(all_x)-1],
            y1=np.max(all_y),
            fillcolor="LightSalmon",
            opacity=0.1,
            layer="below",
            line_width=0,
         )

   fig.update_xaxes(rangemode="nonnegative")
   fig.update_yaxes(rangemode="nonnegative")

   fig.update_layout(
      margin=dict(l=30, r=20, t=5, b=20),   # Top 0 with no title
      paper_bgcolor='rgba(255,255,255,1)',
      plot_bgcolor='rgba(255,255,255,1)',
      showlegend= False,
      autosize=False
   )

   if(_type == 'test_pos_p' or _type == 'mortality'):
      fig.update_layout(yaxis=dict(tickformat=".2%"))

   if(county ==""):
      fig.write_image(folder + os.sep + state + ".png", width=350, height=350)
      print("Graph for " + state + ' Cases (' +  _color + ') created')
   elif('for_a_state' in county):
      tmp = county.split('|')[1]
      fig.write_image(folder + os.sep + tmp + ".png", width=350, height=350)
      print("Graph for " + state + '  ' +  tmp + '  created')
   else:
      fig.write_image(folder + os.sep + county + ".png", width=350, height=350)
      print("Graph for " + county + ", " + state + ' Cases (' +  _color + ') created')

   if(large is True):
      # We also save the larger version of the graph
      fig.update_layout(
         margin=dict(l=30, r=20, t=0, b=20),   # Top 0 with no title
         paper_bgcolor='rgba(255,255,255,1)',
         plot_bgcolor='rgba(255,255,255,1)',
         showlegend= False,
      )

      if(county ==""):
         fig.write_image(folder + os.sep + state + "_lg.png", width=1084, height=450)
         print("Graph for " + state + ' Cases (' +  _color + ') Larger Version created')
      elif('for_a_state' in county):
         tmp = county.split('|')[1]
         fig.write_image(folder + os.sep + tmp + "_lg.png", width=1084, height=450)
         print("Graph for " + state + '  ' +  tmp + ' Larger Version created')
      else:
         fig.write_image(folder + os.sep + county + "_lg.png", width=1084, height=450)
         print("Graph for " + county + ", " + state + ' Cases (' +  _color + ') Larger Version created')

def main_menu():
   print("---------------")
   print(" Enter the attributes of the graph ")
   print("---------------")
   state = input("State Code (ex: AK): ")
   _color = input("Color ('r' for red ,'g' for green, 'o' for orange, 'b' for black):")
   generate_graph_with_avg(state, 'cases', _color, PATH_TO_STATES_FOLDER + os.sep + state, '')

if __name__ == "__main__":
   #os.system("clear")
   #main_menu()
   #generate_graph_with_avg("FL", 'test_pos_p', "r", PATH_TO_STATES_FOLDER + os.sep + "FL"  + os.sep , 'for_a_state|test_pos_p')
   #generate_graph_with_avg("FL", 'test_pos_p', "r", PATH_TO_STATES_FOLDER + os.sep + "FL"  + os.sep , 'for_a_state|test_pos_p'  , True)

   generate_graph_with_avg("TX", 'mortality', "r", PATH_TO_STATES_FOLDER + os.sep + "TX"  + os.sep , 'for_a_state|mortaliy')

   generate_graph_with_avg("TX", 'mortality', "r", PATH_TO_STATES_FOLDER + os.sep + "TX"  + os.sep , 'for_a_state|mortaliy', True)
   #generate_graph_with_avg("CA", 'mortality', "r", PATH_TO_STATES_FOLDER + os.sep + "CA"  + os.sep , 'for_a_state|mortaliy', True)

   #generate_large_graph_with_avg("CA",  "r", PATH_TO_STATES_FOLDER + os.sep + "CA"  + os.sep)

   #generate_dual_graph_X_and_cases("FL", "r", PATH_TO_STATES_FOLDER + os.sep + "FL" ,
   #          {'_type':'test','name':'Tests','raw_line':True},
   #          {'_type':'cases','name':'Cases','avg':[7,3], 'raw_line':False},
   #          output_ext = '_blg.png', # Just for Mike
   #          large=True)

   #generate_dual_graph_X_and_cases("CA", "r", PATH_TO_STATES_FOLDER + os.sep + "CA"  + os.sep,
   #         {'_type':'test','name':'7D. Avg Tests','avg':[7],'raw_line':False},
   #         {'_type':'cases','name':'7D. Avg Cases','avg':[7],'raw_line':False},
   #         output_ext = '_tac.png', # TAC = test & case
   #         large=False)

   # generate_dual_graph_X_and_cases("CA", "r", PATH_TO_STATES_FOLDER + os.sep + "CA",
   #          {'_type':'act_hosp','name':'7D. Avg Hospitalization','avg':[7],'raw_line':False},
   #          {'_type':'cases','name':'7D. Avg Cases','avg':[7],'raw_line':False},
   #          output_ext = '_hospi_and_cases.png',
   #          large=False)

ext: py

max_stars_repo_name: mikehankey/mikehankey

max_issues_repo_name: mikehankey/mikehankey

max_forks_repo_name: mikehankey/mikehankey

reason: import numpy


text:
"""
giasim.py
Author: Samuel B. Kachuck
Date  : 08 01 2015

    Main module for the giapy package. It provides the class GiaSimGlobal,
    which contains the method for computing the glacial isostatic adjustment to
    an earth with a given rheological decay spectrum to an ice history model.

Methods
-------
configure_giasim

Classes
-------
GiaSimGlobal
GiaSimOutput

"""
from __future__ import division

import numpy as np
import spharm
import subprocess
try:
    from progressbar import ProgressBar, Percentage, Bar, ETA
except:
    pass

from giapy.map_tools import GridObject, sealevelChangeByMelt,\
                    volumeChangeLoad, sealevelChangeByUplift, oceanUpliftLoad,\
                    floatingIceRedistribute

from giapy import GITVERSION, timestamp, MODPATH, call, os

class GiaSimGlobal(object):
    def __init__(self, earth, ice, grid=None, topo=None):
        """
        Compute glacial isostatic adjustment on a globe.

        Paramaters
        ----------
        earth : <giapy.code.earth_tools.earthSpherical.SphericalEarth>
        ice   : <giapy.code.icehistory.IceHistory / PersistentIceHistory>
        grid  : <giapy.code.map_tools.GridObject>
        topo  : numpy.ndarray

        Methods
        -------
        performConvolution

        """

        self.earth = earth

        self.ice = ice
        self.nlat, self.nlon = ice.shape

        # The grid used is a cylindrical projection (equispaced lat/lon grid
        # unless otherwise specified)
        if grid is not None:
            self.grid = grid
        else:
            self.grid = GridObject(mapparam={'projection': 'cyl'},
                                    shape=ice.shape)

        self.topo = topo

        # Precompute and store harmonic transform coefficients, for
        # computational efficiency, but at a memory cost.
        self.harmTrans = spharm.Spharmt(self.nlon, self.nlat, legfunc='stored')

    def performConvolution(self, out_times=None, ntrunc=None, topo=None,
                            verbose=False, eliter=5, nrem=1, massconerr=1e-2):
        """Convolve an ice load and an earth response model in fft space.
        Calculate the uplift associated with stored earth and ice model.

        Parameters
        ----------
        out_times : an array of times at which to caluclate the convolution.
                    (default is to use previously stored values).
        ntrunc : int
           The truncation number of the spherical harmonic expansion. Default
           is from the earth model, must be <= ice model's and < grid.nlat
        topo   : array
            Topography on which to compute. If None (default), assumes a flat
            topography. Must be the same shapt as ice.
        verbose : boolean
            Display progress on computation. Depends on progressbar module.
            Default is False.
        eliter  : int
            The maximum number of iterations allowed to compute initial elastic
            response to redistributed load at each stage. If 0, instantaneous
            elastic response is not computed. Default 5.
        nrem   : int
            Number of removal stages between the provided ice stages
            (intermediate steps are interpolated linearly). Default 1.

        Results
        -------
        observerDict : GiaSimOutput
            A dictionary whose keys are fields of interest, such as uplift
            ('upl'), geoid ('geo'), and solid surface topography ('sstopo'),
            computed on the input grid at out_times.
        """

        DENICE      = 931.   #934.           # kg/m^3
        DENWAT      = 1000.  #999.           # kg/m^3
        DENSEA      = 1000.  #1029.          # kg/m^3
        GSURF       = 9.815          # m/s^2
        #PAperM      = DENSEA*GSURF
        NREM        = nrem           # number of intermediate steps

        earth = self.earth
        ice = self.ice
        grid = self.grid
        if topo is None and self.topo is not None:
            topo = self.topo
        if topo is not None:
            assert topo.shape == ice.shape, 'Topo and Ice must have the same shape'

        # Resolution
        ntrunc = ntrunc or min(earth.nmax, ice.nlat-1)
        assert ntrunc <= ice.nlat-1, 'ntrunc > ice.nlat-1'
        ms, ns = spharm.getspecindx(ice.nlat-1)
        # npad is the indices in the larger (padded) array of spherical
        # harmonics that correspond to the smaller (response) array.
        npad = (ns <= ntrunc)

        # Store out_times
        if out_times is None:
            out_times = self.out_times
        else:
            out_times = out_times
        self.out_times = out_times
        assert out_times is not None, 'out_times is not set'

        # Calculate times of intermediate removal stages.
        diffs = np.diff(ice.times)
        addRemovalTimes = []
        for i in range(1, NREM+1):
            addRemovalTimes.append(ice.times[:-1]+i*diffs/NREM)
        addRemovalTimes = np.array(addRemovalTimes).flatten()
        remTimes = np.union1d(ice.times, addRemovalTimes)[::-1]
        calcTimes = np.union1d(remTimes, out_times)[::-1]

        # Initialize output observer
        observerDict = initialize_output(self, out_times, calcTimes, ice.nlat-1,
                                            ntrunc, ns, ice.shape)

        for o in observerDict:
            o.loadStageUpdate(ice.times[0], sstopo=topo)

        esl = 0                 # Equivalent sea level assumed to start at 0.

        elRespArray = earth.getResp(0.)
        ssResp = np.zeros_like(ns)
        ssResp[npad] = observerDict['SS'].isolateRespArray(elRespArray)

        # Convolve each ice stage to the each output time.
        # Primary loop: over ice load changes.
        for icea, ta, iceb, tb in ice.pairIter():
            ################### LOAD STAGE CALCULATION ###################
            # Determine the water load redistribution for ice, uplift, and
            # geoid changes between ta and tb,
            if topo is not None:
                # Get index for starting time.
                nta = observerDict['SS'].locateByTime(ta)
                # Collect the solid-surface topography at beginning of step.
                Ta = observerDict['sstopo'].array[nta]

                # Redistribute the ocean by change in ocean floor / surface.
                ssa, ssb = observerDict['SS'].array[[nta, nta+1]]
                dSS = self.harmTrans.spectogrd(ssb-ssa)
                dhwBarU = sealevelChangeByUplift(dSS, Ta+DENICE/DENSEA*icea,
                                                        grid)
                dhwU = oceanUpliftLoad(dhwBarU, Ta+DENICE/DENSEA*icea, dSS)

                # Update the solid-surface topography with uplift / geoid.
                Tb = Ta + dSS - dhwBarU
                esl += dhwBarU
                dLoad = dhwU.copy()
                dwLoad = dhwU.copy()                # Save the water load

                # Redistribute ice, consistent with current floating ice.
                dILoad, dhwBarI = floatingIceRedistribute(icea, iceb, Tb, grid,
                                                            DENICE/DENSEA)

                # Combine loads from ocean changes and ice volume changes.
                dLoad += dILoad
                esl += dhwBarI
                dwLoad += volumeChangeLoad(dhwBarI, Tb+DENICE/DENSEA*iceb)
                Tb -= dhwBarI

                # Calculate instantaneous (elastic and gravity) responses to
                # the load shift and redistribute ocean accordingly.
                # Note: WE DO NOT CURRENTLY RECHECK FOR FLOATING ICE LOADS.
                if eliter:
                    # Get elastic and geoid response to the water load.
                    # Find the elastic uplift in response to stage's load
                    # redistribution.
                    dSSel = self.harmTrans.spectogrd((ssResp)*\
                                self.harmTrans.grdtospec(dLoad))

                    dhwBarUel = sealevelChangeByUplift(dSSel,
                                                        Tb+DENICE/DENSEA*iceb, grid)
                    dhwUel = oceanUpliftLoad(dhwBarUel,
                                                Tb+DENICE/DENSEA*iceb, dSSel)

                    Tb = Tb + dSSel - dhwBarUel
                    esl += dhwBarUel
                    dLoad = dLoad + dhwUel
                    dwLoad += dhwUel

                    # Iterate elastic responses until they are sufficiently small.
                    for i in range(eliter):
                        # Need to save elastic uplift and geoid at each iteration
                        # to compare to previous steps for convergence.
                        dSSelp = self.harmTrans.spectogrd((ssResp)*\
                                    self.harmTrans.grdtospec(dhwUel))

                        dhwBarUel = sealevelChangeByUplift(dSSelp,
                                                            Tb+DENICE/DENSEA*iceb, grid)
                        dhwUel = oceanUpliftLoad(dhwBarUel,
                                                    Tb+DENICE/DENSEA*iceb, dSSelp)

                        # Correct topography
                        Tb = Tb + dSSelp - dhwBarUel
                        esl += dhwBarUel
                        dLoad = dLoad + dhwUel
                        dwLoad += dhwUel

                        # Truncation error from further iteration
                        err = np.mean(np.abs(dSSelp))/np.mean(np.abs(dSSel))
                        if err <= massconerr:
                            break
                        else:
                            dSSel = dSSel + dSSelp

                            continue

                observerDict['SS'].array[nta+1] += self.harmTrans.grdtospec(dSSel)

                for o in observerDict:
                    # Topography and load for time tb are updated and saved.
                    o.loadStageUpdate(tb, dLoad=dLoad,
                                      topo=Tb+iceb*(Tb + DENICE/DENSEA*iceb>=0),
                                      esl=esl, dwLoad=dwLoad, sstopo=Tb)

            else:
                dLoad = (iceb-icea)*DENICE/DENSEA
                Tb = None

                for o in observerDict:
                    # Topography and load for time tb are updated and saved.
                    o.loadStageUpdate(tb, dLoad=dLoad)

            # Transform load change into spherical harmonics.
            loadChangeSpec = self.harmTrans.grdtospec(dLoad)/NREM

            # Check for mass conservation.
            massConCheck = np.abs(loadChangeSpec[0])/np.abs(loadChangeSpec.max())
            if  verbose and massConCheck >= massconerr:
                print("Load at {0} doesn't conserve mass: {1}.".format(ta,
                                                                massConCheck))
            # N.B. the n=0 load should be zero in cases of glacial isostasy, as
            # mass is conserved during redistribution.

            ################# RESPONSE STAGE CALCULATION #################
            # Secondary loop: over output times.
            for inter_time in np.linspace(tb, ta, NREM, endpoint=False)[::-1]:
                # Perform the time convolution for each output time
                for t_out in calcTimes[calcTimes < inter_time]:
                    respArray = earth.getResp(inter_time-t_out)
                    for o in observerDict:
                        o.respStageUpdate(t_out, respArray,
                                            DENSEA*loadChangeSpec)

        # Don't keep the intermediate uplift stages for water redistribution
        #observerDict.removeObserver('eslUpl', 'eslGeo')

        return observerDict

def configure_giasim(configdict=None):
    """
    Convenience function for setting up a GiaSimGlobal object.

    The function uses inputs stored in giapy/data/inputs/, which can be
    downloaded from pages.physics.cornell.edu/~skachuck/giainputs.tar.gz or
    using the script giapy/dldata.

    Parameters
    ----------
    configdict : dict
        A dictionary with keys 'earth' and 'ice' and values of the
        names of these inputs. The key 'topo' is optional, and defaults to a
        flat initial topography.

    sim        : GiaSimGlobal object
    """

    DEFAULTCONFIG = {'earth': '75km0p04Asth_4e23Lith',
                     'ice'  : 'AA2_Tail_nochange5_hightres_Pers_288_square',
                     'topo' : 'sstopo288'}

    configdict = configdict or DEFAULTCONFIG

    assert configdict.has_key('earth'), 'GiaSimGlobal needs earth specified'
    assert configdict.has_key('ice'), 'GiaSimGlobal needs ice specified'

    #ppath = os.path.dirname(os.path.split(__file__)[0])
    dpath = MODPATH + '/data/inputs/'
    ename = configdict['earth']+'.earth'
    iname = configdict['ice']+'.ice'

    filecheck = os.path.isfile(dpath+ename) and os.path.isfile(dpath+iname)

    topo = configdict.get('topo', None)
    if topo is not None:
        tname = topo+'.topo'
        filecheck = filecheck and os.path.isfile(dpath+tname)

    if not filecheck:
        with open(MODPATH+'/data/inputfilelist', 'r') as f:
            filelist = f.read().splitlines()
        if ename in filelist and iname in filelist and tname in filelist:
            print('Downloading inputs')
            p = subprocess.Popen(os.path.join(MODPATH+'/dldata'), shell=True,
                                cwd=ppath)
            p.wait()
        else:
            raise IOError('One of the inputs you specified does not exist.')

    earth = np.load(dpath+ename)
    ice = np.load(dpath+iname)
    if topo is not None:
        topo = np.load(dpath+tname)[2]

    sim = GiaSimGlobal(earth=earth, ice=ice, topo=topo)

    return sim

def initialize_output(sim, out_times, calcTimes, nmax, ntrunc, ns, shape):
    earth = sim.earth
    # Initialize the return object to include...
    # ... values desired at output times
    #   [1] Uplift
    uplObserver = earth.TotalUpliftObserver(out_times, nmax, ntrunc, ns)
    #   [2] Horizontal deformation
    horObserver = earth.TotalHorizontalObserver(out_times, nmax, ntrunc, ns)
    #   [3] Uplift velocities
    velObserver = earth.VelObserver(out_times, nmax, ntrunc, ns)
    #   [4] Geoid perturbations
    geoObserver = earth.GeoidObserver(out_times, nmax, ntrunc, ns)
    #   [5] Gravitational acceleration perturbations
    gravObserver = earth.GravObserver(out_times, nmax, ntrunc, ns)

    # ... and values needed to perform the convolution
    #   [1] Uplift for ocean redistribution
    SeaSurfaceObserver = earth.SeaSurfaceObserver(calcTimes, nmax, ntrunc, ns)
    #   [2] Geoid for ocean redistribution
    eslGeoObserver = earth.GeoidObserver(calcTimes, nmax, ntrunc, ns)
    #   [3] Topography (to top of ice) to find floating ice
    topoObserver = HeightObserver(calcTimes, shape, 'topo')
    #   [4] Load (total water + ice load in water equivalent)
    loadObserver = HeightObserver(calcTimes, shape, 'dLoad')
    #   [5] Water load
    wloadObserver = HeightObserver(calcTimes, shape, 'dwLoad')
    #   [6] Solid surface topography for ocean redistribution
    rslObserver = HeightObserver(calcTimes, shape, 'sstopo')
    #   [7] Eustatic sea level, with average uplift and geoid over oceans.
    eslObserver = EslObserver(calcTimes)

    observerDict = GiaSimOutput(sim)
    observerDict.addObserver('upl'   , uplObserver)
    observerDict.addObserver('hor'   , horObserver)
    observerDict.addObserver('grav'  , gravObserver)
    observerDict.addObserver('load'  , loadObserver)
    observerDict.addObserver('wload' , wloadObserver)
    observerDict.addObserver('vel'   , velObserver)
    observerDict.addObserver('geo'   , geoObserver)

    observerDict.addObserver('SS', SeaSurfaceObserver)
    observerDict.addObserver('topo'  , topoObserver)
    observerDict.addObserver('esl'   , eslObserver)
    observerDict.addObserver('sstopo', rslObserver)
    return observerDict

class GiaSimOutput(object):
    """A container object for computations of glacial isostatic adjustment.

    The results of the GIA computation can be accessed via object attributes or
    via a dictionary (e.g., if iteration is desired). The object's __repr__
    gives the computed attributes.

    Parameters
    ----------
    inputs : anything
        anything desired to be stored as inputs. Typically, this will be a
        GiaSimGlobal object that contains the ice load, earth model, and other
        necessary objects for reproducing the computation.

    Methods
    -------
    addObserver - add an observer to the watchlist
    removeObserver - remove an observer from the watchlist
    transformObservers - transform each observer in the watchlist, using each's
        own transform function (must be provided by observer).

    Data
    ----
    GITVERSION : the hash of the current giapy git HEAD
    TIMESTAMP : the datetime of creation (calculation)
    """
    def __init__(self, inputs):
        self.GITVERSION = GITVERSION
        self.TIMESTAMP = timestamp()
        self.inputs = inputs
        self._observerDict = {}

    def __getitem__(self, key):
        return self._observerDict.__getitem__(key)

    def __iter__(self):
        return self._observerDict.itervalues()

    def __repr__(self):
        retstr = ''
        retstr = 'GIA computed on {}\n'.format(self.TIMESTAMP)
        retstr += 'Has observers: '+', '.join(self._observerDict)
        return retstr

    def addObserver(self, name, observer):
        self._observerDict[name] = observer
        setattr(self, name, observer)

    def removeObserver(self, *names):
        for name in names:
            del self._observerDict[name]
            delattr(self, name)

    def transformObservers(self, inverse=False):
        for obs in self:
            obs.transform(self.inputs.harmTrans, inverse=inverse)

class AbstractGiaSimObserver(object):
    """The GiaSimObserver mediates between an earth model's response function
    and the convolved result. It needs to store the harmonic response

    initialize and update must be overwritten for class to work, all others
    simply pass.
    """
    def __init__(self):
        pass

    def __getitem__(self, key):
        return self.array.__getitem__(key)

    def __iter__(self):
        return self.array.__iter__()

    @property
    def shape(self):
        return self.array.shape

    def initialize(self, outTimes, ntrunc):
        raise NotImplemented

    def loadStageUpdate(self, *args, **kwargs):
        pass

    def respStageUpdate(self, *args, **kwargs):
        pass

    def update(self, tout, tdur, dLoad):
        raise NotImplemented

    def transform(self, trans, inverse=True):
        pass

    def locateByTime(self, time):
        if time not in self.outTimes:
            raise ValueError('time not in self.outTimes')
        return np.argwhere(self.outTimes == time)[0][0]

    def nearest_to(self, time):
        """Return a field from outTimes nearest to time.
        """
        idx = (np.abs(self.outTimes-time)).argmin()
        return self[idx]

class AbstractEarthGiaSimObserver(AbstractGiaSimObserver):
    """Abstract class for results in spherical harmonic space, updated during
    the response stage.

    Must implement isolateRespArray to pull proper response curve from the
    computed earth model.
    """
    def __init__(self, outTimes, nmax, ntrunc, ns):
        self.initialize(outTimes, nmax, ns)
        self.npad = (ns <= ntrunc)
        self.ns = ns[self.npad]
        self.spectral = True

    def initialize(self, outTimes, ntrunc, ns):
        self.array = np.zeros((len(outTimes),
                               int((ntrunc+1)*(ntrunc+2)/2)), dtype=complex)
        self.outTimes = outTimes
        self.ns = ns

    def respStageUpdate(self, *args, **kwargs):
        self.update(*args, **kwargs)

    def update(self, tout, respArray, dLoad):
        if tout not in self.outTimes:
            return
        n = self.locateByTime(tout)
        resp = self.isolateRespArray(respArray)
        self.array[n][self.npad] += resp * dLoad[self.npad]

    def transform(self, trans, inverse=True):
        if not inverse and self.spectral:
            self.array = trans.spectogrd(self.array.T).T
            self.spectral = False
        elif inverse and not self.spectral:
            self.array = trans.grdtospec(self.array.T).T
            self.spectral = True

    def isolateRespArray(self, respArray):
        raise NotImplemented()

class HeightObserver(AbstractGiaSimObserver):
    """General observer for heights computed on the real-space grid and updated
    during the loadStage.
    """
    def __init__(self, outTimes, iceShape, name):
        self.initialize(outTimes, iceShape)
        self.name = name

    def initialize(self, outTimes, iceShape):
        self.array = np.zeros((len(outTimes),
                                iceShape[0], iceShape[1]))
        self.outTimes = outTimes

    def loadStageUpdate(self, tout, **kwargs):
        if self.name in kwargs.keys():
            self.update(tout, kwargs[self.name])

    def update(self, tout, load):
        if tout not in self.outTimes:
            return
        n = self.locateByTime(tout)
        self.array[n] = load

class EslObserver(AbstractGiaSimObserver):
    def __init__(self, outTimes):
        self.array = np.zeros(len(outTimes))
        self.outTimes = outTimes

    def loadStageUpdate(self, tout, **kwargs):
        if 'esl' in kwargs.keys():
            self.update(tout, kwargs['esl'])

    def update(self, tout, esl):
        if tout not in self.outTimes:
            return
        n = self.locateByTime(tout)
        self.array[n] = esl

max_stars_repo_name: cameronbook/giapy

max_issues_repo_name: cameronbook/giapy

max_forks_repo_name: cameronbook/giapy


text:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cbook as cbook
from matplotlib_scalebar.scalebar import ScaleBar

with cbook.get_sample_data("s1045.ima.gz") as dfile:
    im = np.frombuffer(dfile.read(), np.uint16).reshape((256, 256))

fig, ax = plt.subplots()
ax.axis("off")

ax.imshow(im, cmap="gray", extent=[0, 20.48, 0, 20.48])

scalebar = ScaleBar(1, "cm", length_fraction=0.25)
ax.add_artist(scalebar)

fig.savefig("argument_dx.png", dpi=60, bbox_inches="tight")


text:
# -*- coding: utf-8 -*-
#from pyramid.arima import auto_arima
import numpy as np
import logging
import sys
from statsmodels.tsa.api import ExponentialSmoothing, SimpleExpSmoothing, Holt
from statsmodels.tsa.arima_model import ARIMA
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import acf
import matplotlib.pylab as plt
#from fbprophet import Prophet
#from tbats import BATS, TBATS
from datetime import datetime, timedelta
from dateutil.relativedelta import relativedelta
import pandas as pd

def getClosest_DateCode(ref_data_seg, current_dateCode, last_year_date):
    """
    This function returns the closest day with the same code according to last year reference data
    :param current_date: this is the current date. The starting point to make forecast
    :param current_dateCode:  code into forecast date
    :param last_year_date: last year date (reference data)
    :return: day (date): closest day
    """
    i = 0
    first = True
    code1 = 0
    code2 = 0
    day_plus = None
    day_minus = None

    if last_year_date.year == ref_data_seg['FECHA'].iloc[0].year:
        while((code1!=current_dateCode) & (code2 != current_dateCode)):
            if first: # TODO: refractor this part of code and put at the begining of the function
                code1 = ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == last_year_date]['COD_LABORALIDAD'].item()
                first = False
                i = i+1
                day_plus = pd.to_datetime(last_year_date)
                day_minus = pd.to_datetime(last_year_date)
            else:
                try:
                    day_plus = day_plus + timedelta(days=i)
                    if (day_plus.year == last_year_date.year):
                        if len(ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == day_plus]['COD_LABORALIDAD'].index) == 0:
                            code1 = current_dateCode
                        else:
                            code1 = ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == day_plus]['COD_LABORALIDAD'].iloc[0]

                    day_minus =  day_minus - timedelta(days=i)
                    if(day_minus.year == last_year_date.year):
                        if len(ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == day_minus]['COD_LABORALIDAD'].index)==0:
                            code2=current_dateCode
                        else:
                            code2 = ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == day_minus]['COD_LABORALIDAD'].iloc[0]
                except OverflowError as err:
                    print("Not found day with same code in last year data, using the same day")
                    return last_year_date

        if code2==current_dateCode:
            return day_minus
        else:
            return day_plus
    else:
        logging.error("No reference data available for the descriptive model")
        sys.exit(0)

def descriptive_model(ref_data_seg, current_date, segment_id, current_dateCode, total_vehiclesBuffer, smooth_mode, forecast_horizon):
    """
    :param reference_data: it could be the collected anual data, the collected anual data using the monthly code average
    :param current_date: this is the current date. The starting point to make forecast
    :param segment_id:  segment_id selected in forecast (forecast are made by segment)
    :param current_dateCode:  code into forecast date
    :param total_vehiclesBuffer: n old values (before forecast date)
    :param smooth_mode: type of smoothing method
    :param forecast_horizon: number of predictions to be made
    :return: predictions (list): a list with the forecast values
    """
    last_year_date= current_date - relativedelta(years=1)
    if len(ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) == last_year_date].index) == 0: # NO DATA AVAILABLE FOR THAT DAY
        closest_ref_date=ref_data_seg.loc[pd.to_datetime(ref_data_seg['FECHA']) >= last_year_date]['FECHA'].iloc[0] # get first available
    else:
        closest_ref_date=getClosest_DateCode(ref_data_seg, current_dateCode, last_year_date)

    predictions = ref_data_seg.loc[(pd.to_datetime(ref_data_seg['FECHA'])>=closest_ref_date)&(pd.to_datetime(ref_data_seg['FECHA'])<=(closest_ref_date+ timedelta(minutes=forecast_horizon)))]
    return predictions['TOTAL_VEHICULOS'].values

def stats_model(reference_data, current_date, segment_id, current_dateCode, total_vehiclesBuffer, smooth_mode, forecast_horizon, algo_name):
    """
    :param reference_data: it could be the collected anual data, the collected anual data using the monthly code average
    :param current_date: this is the current date. The starting point to make forecast
    :param segment_id:  segment_id selected in forecast (forecast are made by segment)
    :param current_dateCode:  code into forecast date
    :param total_vehiclesBuffer: n old values (before forecast date)
    :param smooth_mode: type of smoothing method
    :param forecast_horizon: number of predictions to be made
    :return: predictions (list): a list with the forecast values
    """
    predictions=[]
    if algo_name=='fbprophet':
        predictions = reference_data.loc[(reference_data.MES == current_date.month) & (reference_data.ID_SEGMENT == segment_id)  & ((reference_data.DAY*24*60)+(reference_data.HOUR*60 + reference_data.MINUTE) >=((current_date.day*24*60)+current_date.hour*60 +current_date.minute))]['TOTAL_VEHICULOS'].values
    else:
        ref_data = reference_data.loc[(reference_data.MES==current_date.month)&(reference_data.ID_SEGMENT==segment_id)&(reference_data.COD_LABORALIDAD==current_dateCode)]
        if current_date.hour >=22: # circle hour
            predictions=[]
            predictions.extend(ref_data.loc[(pd.to_datetime(ref_data['TIME']).dt.hour * 60 + (pd.to_datetime(ref_data['TIME']).dt.minute)) >= (current_date.hour * 60 + current_date.minute)]['TOTAL_VEHICULOS'].values)
            predictions.extend(ref_data.loc[(pd.to_datetime(ref_data['TIME']).dt.hour >= 0) & (pd.to_datetime(ref_data['TIME']).dt.hour <= 2)]['TOTAL_VEHICULOS'].values) # take 2 hours more to complete
        else:
            predictions = ref_data.loc[(pd.to_datetime(ref_data['TIME']).dt.hour*60+(pd.to_datetime(ref_data['TIME']).dt.minute))>=(current_date.hour*60+current_date.minute)]['TOTAL_VEHICULOS'].values
    return predictions

def auto_arima_model(data):
    """ This method train an ARIMA model using the maximum likelihood method to find automatically the best paramaters.
        See pkg documentation, for implementation details.

    Args:
        data (dataframe): The training dataset used to train the model

    Returns:
        model: The trained model

    """
    """   stepwise_model = auto_arima(data, start_p=2, start_q=2,
                           max_p=3, max_q=3,m=96,
                           start_P=0, seasonal=True,
                           d=1, D=1, trace=True,
                           error_action='ignore',
                           suppress_warnings=True,
                           stepwise=True)
    """
    return None #stepwise_model

def arima_model(train):
    """ This method train an ARIMA model
    Args:
        data (dataframe): The training dataset used to train the model

    Returns:
        model: The trained model

    """
    #model = ARIMA(train, order=(96, 0, 1))
    model = sm.tsa.statespace.SARIMAX(train, trend='ct', order=(1, 1, 1), seasonal_order=(1, 1, 1, 96), enforce_stationarity=True, enforce_invertibility=False)
    return model

def prophet_model(df_training):
    """ This method train an FBPROPHET model using the well-known developed pkg by facebook called fbprophet.
        See https://facebook.github.io/prophet/docs/quick_start.html for implementation details
    Args:
        df_training (dataframe): The training dataset used to train the model

    Returns:
        model: The trained model
    """
    df_training['ds'] = df_training.index
    df_training['y'] = df_training['TOTAL_VEHICULOS']
    train_cols = df_training[['ds', 'y']]
    model = None#Prophet()  # instantiate Prophet
    model.fit(train_cols);  # fit the model with your dataframe
    return model

def holtwinters_exp_smoothing_model(train):
    """ This method train an Exponential Smoothing model (Holt-Winters)
        See https://otexts.com/fpp2/holt-winters.html for implementation details
    Args:
        df_training (dataframe): The training dataset used to train the model

    Returns:
        model: The trained model
    """
    ses_model = ExponentialSmoothing(np.asarray(train), seasonal_periods=int(96), seasonal='add').fit()
    return ses_model

def tbats_model():
    """ This method train an TBATs model
    Args:
        df_training (dataframe): The training dataset used to train the model

    Returns:
        model: The trained model
    """
    model = None#TBATS(seasonal_periods=[14, 30.5])
    return model

def seasonal_decompose_model(train):
    """ Performns a seasonal decomposition of the training dataset passed as argument program
    Args:
        df_training (dataframe): The training dataset used to train the model

    Returns:
        bool: The return value. True for success, False otherwise
    """
    sm.tsa.seasonal_decompose(train.Count).plot()
    result = sm.tsa.stattools.adfuller(train.Count)
    plt.show()
    return True

def autocorrelation_function(train):
    """ Plots the autocorrelation from the training dataset
    Args:
        df_training (dataframe): The training dataset used to train the model

    Returns:
        bool: The return value. True for success, False otherwise
    """
    #https://www.vinta.com.br/blog/2018/understanding-time-series-forecasting-python/
    data = np.log(train)
    lac_acf = acf(data, nlags=40)

    plt.figure(figsize=(15, 5))
    plt.subplot(121)
    plt.stem(lac_acf)
    plt.axhline(y=0, linestyle='-', color='black')
    plt.xlabel('Lag')
    plt.ylabel('ACF')
    plt.show()
    return True

reason: import numpy,import statsmodels,from statsmodels


text:
import os
import shutil
from contextlib import ContextDecorator
from functools import wraps
from pathlib import Path
from tempfile import mkdtemp
from typing import Dict, List, Optional, Union

import numpy as np
from skimage.io import imsave

from .shapes import random_shapes

def get_dataset_path(dataset: str = "MVTec") -> str:
    """Selects path based on tests in local system or docker image.

    Local install assumes datasets are located in anomaly/datasets/.
    In either case, if the location is empty, the dataset is downloaded again.
    This speeds up tests in docker images where dataset is already stored in /tmp/anomalib

    Example:
    Assume that `datasets directory exists in ~/anomalib/,

    >>> get_dataset_path(dataset="MVTec")
    './datasets/MVTec'

    """
    # Initially check if `datasets` directory exists locally and look
    # for the `dataset`. This is useful for local testing.
    path = os.path.join("./datasets", dataset)

    # For docker deployment or a CI that runs on server, dataset directory
    # may not necessarily be located in the repo. Therefore, check anomalib
    # dataset path environment variable.
    if not os.path.isdir(path):
        path = os.path.join(os.environ["ANOMALIB_DATASET_PATH"], dataset)
    return path

def generate_random_anomaly_image(
    image_width: int,
    image_height: int,
    shapes: List[str] = ["triangle", "rectangle"],
    max_shapes: Optional[int] = 5,
    generate_mask: Optional[bool] = False,
) -> Dict:
    """Generate a random image with the corresponding shape entities.

    Args:
        image_width (int): Width of the image
        image_height (int): Height of the image
        shapes (List[str]): List of shapes to draw in the image. Make sure these are different from `anomalous_shapes`
        max_shapes (int): Maximum shapes of a kind in the image. Defaults to 5.
        max_size (Optional[int], optional): Maximum size of the test shapes. Defaults to 10.
        generate_mask (bool): Toggle between train/test split. Train images are restricted to first quadrant.
                    Also generates the mask for the image. Defaults to False.
    Returns:
        Tuple: image if `train` is False. Else return image, mask
    """

    image: np.ndarray = np.full((image_height, image_width, 3), 255, dtype=np.uint8)

    input_region = [0, 0, image_width - 1, image_height - 1]

    for shape in shapes:
        shape_image = random_shapes(input_region, (image_height, image_width), max_shapes=max_shapes, shape=shape)
        image = np.minimum(image, shape_image)  # since white is 255

    result = {"image": image}

    if generate_mask:
        mask = np.zeros((image_height, image_width), dtype=np.uint8)
        # if color exists in any channel turn the mask bit white.
        # not sure if there is a better way to do this.
        mask[image[..., 0] < 255] = 255
        mask[image[..., 1] < 255] = 255
        mask[image[..., 2] < 255] = 255
        result["mask"] = mask

    return result

class TestDataset:
    def __init__(
        self,
        num_train: int = 1000,
        num_test: int = 100,
        img_height: int = 128,
        img_width: int = 128,
        max_size: int = 10,
        train_shapes: List[str] = ["triangle", "rectangle"],
        test_shapes: List[str] = ["hexagon", "star"],
        path: Union[str, Path] = "./datasets/MVTec",
        use_mvtec: bool = False,
        seed: int = 0,
    ) -> None:
        """Creates a context for Generating Dummy Dataset. Useful for wrapping test functions.
        NOTE: for MVTec AD dataset it does not return a category.
        It is adviced to use a default parameter in the function

        Args:
            num_train (int, optional): Number of training images to generate. Defaults to 1000.
            num_test (int, optional): Number of testing images to generate per category. Defaults to 100.
            img_height (int, optional): Height of the image. Defaults to 128.
            img_width (int, optional): Width of the image. Defaults to 128.
            max_size (Optional[int], optional): Maximum size of the test shapes. Defaults to 10.
            train_shapes (List[str], optional): List of good shapes. Defaults to ["circle", "rectangle"].
            test_shapes (List[str], optional): List of anomalous shapes. Defaults to ["triangle", "ellipse"].
            path (Union[str, Path], optional): Path to MVTec AD dataset. Defaults to "./datasets/MVTec".
            use_mvtec (bool, optional): Use MVTec AD dataset or dummy dataset. Defaults to False.
            seed (int, optional): Fixes seed if any number greater than 0 is provided. 0 means no seed. Defaults to 0.

        Example:
            >>> @TestDataset
            >>> def test_some_function(path, category="leather"):
            >>>     do something
        """
        self.num_train = num_train
        self.num_test = num_test
        self.img_height = img_height
        self.img_width = img_width
        self.max_size = max_size
        self.train_shapes = train_shapes
        self.test_shapes = test_shapes
        self.path = path
        self.use_mvtec = use_mvtec
        self.seed = seed

    def __call__(self, func):
        @wraps(func)
        def inner(*args, **kwds):
            # If true, will use MVTech AD dataset for testing.
            # Useful for nightly builds
            if self.use_mvtec:
                return func(*args, path=self.path, **kwds)
            else:
                with GeneratedDummyDataset(
                    num_train=self.num_train,
                    num_test=self.num_test,
                    img_height=self.img_height,
                    img_width=self.img_width,
                    train_shapes=self.train_shapes,
                    test_shapes=self.test_shapes,
                    max_size=self.max_size,
                    seed=self.seed,
                ) as dataset_path:
                    kwds["category"] = "shapes"
                    return func(*args, path=dataset_path, **kwds)

        return inner

class GeneratedDummyDataset(ContextDecorator):
    """Context for generating dummy shapes dataset.
    Example:
        >>> with GeneratedDummyDataset(num_train=1000,num_test=100) as dataset_path:
        >>>     some_function()

        Args:
            num_train (int, optional): Number of training images to generate. Defaults to 1000.
            num_test (int, optional): Number of testing images to generate per category. Defaults to 100.
            img_height (int, optional): Height of the image. Defaults to 128.
            img_width (int, optional): Width of the image. Defaults to 128.
            max_size (Optional[int], optional): Maximum size of the test shapes. Defaults to 10.
            train_shapes (List[str], optional): List of good shapes. Defaults to ["circle", "rectangle"].
            test_shapes (List[str], optional): List of anomalous shapes. Defaults to ["triangle", "ellipse"].
            seed (int, optional): Fixes seed if any number greater than 0 is provided. 0 means no seed. Defaults to 0.
    """

    def __init__(
        self,
        num_train: int = 1000,
        num_test: int = 100,
        img_height: int = 128,
        img_width: int = 128,
        max_size: Optional[int] = 10,
        train_shapes: List[str] = ["triangle", "rectangle"],
        test_shapes: List[str] = ["star", "hexagon"],
        seed: int = 0,
    ) -> None:
        self.root_dir = mkdtemp()
        self.num_train = num_train
        self.num_test = num_test
        self.train_shapes = train_shapes
        self.test_shapes = test_shapes
        self.image_height = img_height
        self.image_width = img_width
        self.max_size = max_size
        self.seed = seed

    def _generate_dataset(self):
        """Generates dummy dataset in a temporary directory using the same
        convention as MVTec AD."""
        # create train images
        train_path = os.path.join(self.root_dir, "shapes", "train", "good")
        os.makedirs(train_path, exist_ok=True)
        for i in range(self.num_train):
            result = generate_random_anomaly_image(
                image_width=self.image_width,
                image_height=self.image_height,
                shapes=self.train_shapes,
                generate_mask=False,
            )
            image = result["image"]
            imsave(os.path.join(train_path, f"{i:03}.png"), image, check_contrast=False)

        # create test images
        for test_category in self.test_shapes:
            test_path = os.path.join(self.root_dir, "shapes", "test", test_category)
            mask_path = os.path.join(self.root_dir, "shapes", "ground_truth", test_category)
            os.makedirs(test_path, exist_ok=True)
            os.makedirs(mask_path, exist_ok=True)
            # anomaly and masks. The idea is to superimpose anomalous shapes on top of correct ones
            for i in range(self.num_test):
                correct_shapes = generate_random_anomaly_image(
                    image_width=self.image_width,
                    image_height=self.image_height,
                    shapes=self.train_shapes,
                    generate_mask=False,
                )
                result = generate_random_anomaly_image(
                    image_width=self.image_width,
                    image_height=self.image_height,
                    shapes=[test_category],
                    generate_mask=True,
                )
                correct_shapes = correct_shapes["image"]
                image, mask = result["image"], result["mask"]
                image = np.minimum(image, correct_shapes)  # since 255 is white
                imsave(os.path.join(test_path, f"{i:03}.png"), image, check_contrast=False)
                imsave(os.path.join(mask_path, f"{i:03}_mask.png"), mask, check_contrast=False)
        # good test
        test_good = os.path.join(self.root_dir, "shapes", "test", "good")
        os.makedirs(test_good, exist_ok=True)
        for i in range(self.num_test):
            result = generate_random_anomaly_image(
                image_width=self.image_width,
                image_height=self.image_height,
                shapes=self.train_shapes,
            )
            image = result["image"]
            imsave(os.path.join(test_good, f"{i:03}.png"), image, check_contrast=False)

    def __enter__(self):
        """Creates the dataset in temp folder."""
        if self.seed > 0:
            np.random.seed(self.seed)
        self._generate_dataset()
        return self.root_dir

    def __exit__(self, _exc_type, _exc_value, _exc_traceback):
        """Cleanup the directory."""
        shutil.rmtree(self.root_dir)


text:
import os
import h5py
from click import progressbar
import numpy as np

def read_skeleton_file(file_path: str):
    with open(file_path, 'r') as f:
        skeleton_sequence = {'numFrame': int(f.readline()), 'frameInfo': []}
        for _ in range(skeleton_sequence['numFrame']):
            frame_info = {'numBody': int(f.readline()), 'bodyInfo': []}
            for _ in range(frame_info['numBody']):
                body_info_key = [
                    'bodyID', 'clippedEdges', 'handLeftConfidence',
                    'handLeftState', 'handRightConfidence', 'handRightState',
                    'isResticted', 'leanX', 'leanY', 'trackingState'
                ]
                body_info = {
                    k: float(v)
                    for k, v in zip(body_info_key, f.readline().split())
                }
                body_info['numJoint'] = int(f.readline())
                body_info['jointInfo'] = []
                for _ in range(body_info['numJoint']):
                    joint_info_key = [
                        'x', 'y', 'z', 'depthX', 'depthY', 'colorX', 'colorY',
                        'orientationW', 'orientationX', 'orientationY',
                        'orientationZ', 'trackingState'
                    ]
                    joint_info = {
                        k: float(v)
                        for k, v in zip(joint_info_key, f.readline().split())
                    }
                    body_info['jointInfo'].append(joint_info)
                frame_info['bodyInfo'].append(body_info)
            skeleton_sequence['frameInfo'].append(frame_info)

    return skeleton_sequence

def parse_skeleton_data(file_path: str, joints: int = 25, persons = 2) -> None:

    sequence_data = read_skeleton_file(file_path)
    data = np.zeros((3, sequence_data['numFrame'], joints, persons), dtype=np.float32)

    for frame_number, frame in enumerate(sequence_data['frameInfo']):
        for body_number, body in enumerate(frame['bodyInfo']):
            for joint_number, joint in enumerate(body['jointInfo']):
                if body_number < persons and joint_number < joints:
                    data[:, frame_number, joint_number, body_number] = [joint['x'], joint['y'], joint['z']]

    return np.around(data, decimals=3)

def generate_skeleton_dataset(data_path: str, output_path: str) -> None:

    target_file_path = f"{output_path}/skeleton.h5"

    with h5py.File(target_file_path, "w") as target_file:
        progress_bar = progressbar(iterable=None, length=len(next(os.walk(data_path))[2]))
        for file_name in os.listdir(data_path):
            sequence_name = os.path.splitext(file_name)[0]

            skeleton_data = parse_skeleton_data(f"{data_path}/{file_name}")

            f = open(output_path + "log.txt", "w+")
            f.write(sequence_name)
            f.write("\r\n")
            f.close()

            target_file.create_group(sequence_name).create_dataset("skeleton", data=skeleton_data)
            progress_bar.update(1)

max_stars_repo_name: kanmaytacker/heart-net

max_issues_repo_name: kanmaytacker/heart-net

max_forks_repo_name: kanmaytacker/heart-net


text:
#!/usr/bin/python2.7
import os
from prettytable import PrettyTable
from numpy import median

def parse_time(time):
    minute = float(time.split("m")[0])
    second = float(time.split("m")[1].split("s")[0])
    total = minute * 60 + second
    if total > 300:
        return 300
    else:
        return total

def output(a):
    for element in a:
        print element + "\t",
    print "\n",

def get_num(storage, filename):
    index = 0
    lines = os.popen("cat " + filename).read().split("\n")
    for line in lines:
        if index == 0:
            index = 1
        else:
            storage.append(parse_time(line))
            index = 0

def get_name(storage, filename):
    index = 0
    lines = os.popen("cat " + filename).read().split("\n")
    for line in lines:
        if index == 0:
            index = 1
            storage.append(line)
        else:
            index = 0

def get_average(elements):
    total = 0.0
    count = 0
    for element in elements:
        total += element
        count += 1
    return round(total / count, 1)

def get_median(elements):
    return round(median(elements), 1)

def get_timeout(elements):
    count = 0
    for element in elements:
        if element >= 300:
            count += 1
    return count

def main(name, optimize, string):
    x = PrettyTable(["Lines", "**Average (Opt)**", "Average (Str)"
        , "**Median (Opt)**", "Median (Str)"
        , "**Timeout (Opt)**", "Timeout (Str)"])
    index = 0
    opts = []
    strs = []
    line = 1
    while index < len(name):
        if len(name[index]) == 0:
            index += 1
            continue
        target = int(name[index].split("_")[0].replace("test", ""))

        if target > line:
            x.add_row([str(line * 50)
                , "**" + str(get_average(opts)) + "**", str(get_average(strs))
                , "**" + str(get_median(opts)) + "**", str(get_median(strs))
                , "**" + str(get_timeout(opts)) + "**", str(get_timeout(strs))])
            line = target
            opts = []
            strs = []

        opts.append(optimize[index])
        strs.append(string[index])
        index += 1

    x.add_row([str(line * 50)
        , "**" + str(get_average(opts)) + "**", str(get_average(strs))
        , "**" + str(get_median(opts)) + "**", str(get_median(strs))
        , "**" + str(get_timeout(opts)) + "**", str(get_timeout(strs))])

    print x

if __name__ == "__main__":
    name = []
    optimize = []
    string = []

    get_name(name, "Optimize.log")
    get_num(optimize, "Optimize.log")
    get_num(string, "String.log")

    main(name, optimize, string)

max_stars_repo_name: spencerwuwu/j-ReCoVer

max_issues_repo_name: spencerwuwu/j-ReCoVer

max_forks_repo_name: spencerwuwu/j-ReCoVer

reason: from numpy


text:
# ===============================================================================
# Copyright 2011 Jake Ross
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
# ===============================================================================

# ============= enthought library imports =======================
from __future__ import absolute_import
from .data_manager import DataManager

# ============= standard library imports ========================
import csv
from numpy import loadtxt

# import time
# from matplotlib.dates import num2date
import os
from six.moves import map
from six.moves import zip

# ============= local library imports  ==========================
class CSVDataManager(DataManager):
    """ """

    format_str = "{:0.6f}"
    _writer = None
    _file = None
    delimiter = ","

    def load(self, frame_key=None):
        if frame_key is None:
            frame_key = self._current_frame

        frame = self._get_frame(frame_key)
        if frame is not None:
            #            from pylab import datestr2num
            #            converters = {0:datestr2num}
            #            return loadtxt(frame, converters=converters, delimiter=',')
            return loadtxt(frame, delimiter=",")

    def write_metadata(self, md, frame_key=None):

        sline = ["#<metadata>==================================================="]
        eline = ["#</metadata>==================================================="]
        data = [sline] + md + [eline]
        self.write_to_frame(data, frame_key)

    def write_to_frame(self, datum, frame_key=None):
        if frame_key is None:
            frame_key = self._current_frame

        frame = self._get_frame(frame_key)
        #        print frame, frame_key, self.frames.keys()
        if frame is not None:
            self.new_writer(frame, datum)

    def new_writer(self, p, datum, append=True):
        """ """
        mode = "w"
        if append:
            mode = "a"
        with open(p, mode) as f:
            #            writer = self.writer
            writer = csv.writer(f, delimiter=self.delimiter)
            if isinstance(datum[0], (list, tuple)):
                writer.writerows(datum)
            else:
                try:
                    datum = list(map(self.format_str.format, datum))
                #                    datum = [self.format_str.format(d) for d in datum]
                except Exception:
                    pass

                writer.writerow(datum)

    def get_current_path(self):
        return self.frames[self._current_frame]

    def read_data(self, path=None):
        if path is None:
            path = self.get_current_path()

        if not os.path.isfile(path):
            return

        with open(path, "r") as f:
            reader = csv.reader(f)

            data = [row for row in reader]
            return list(zip(*data))

#    @property
#    def writer(self):
#        if self._writer is None:
#            self._writer = csv.writer(self._file)
#
#        return self._writer

# if __name__ == '__main__':
#    d = CSVDataManager()
#    d.new_frame()
#    for i in range(3):
#        d.add_time_stamped_value(i)
#        time.sleep(1)
#    x, y = d.load().transpose()
#    print x
#    for xi in x:
#        print num2date(xi).second

# ============= EOF ====================================

max_stars_repo_name: ASUPychron/pychron

max_issues_repo_name: ASUPychron/pychron

max_forks_repo_name: UIllinoisHALPychron/pychron


text:
"""
Author:         Anonymous
Description:
                Data generation functions
                - static data (synthetic and MNIST)
                - dynamic data, control policy outputs (displacement and NN)

                Here the sampling/exploration strategy is implemented
"""

import logging
import itertools
import numpy as np

from operator import itemgetter

from behaviour_representations.utils.utils import _TAB, _SEED
from behaviour_representations.utils.utils import info_outcome, timing

from behaviour_representations.exploration.data_generate_ps import DisplacementParamSampling
from behaviour_representations.exploration.data_generate_ls import LatentSampling
from behaviour_representations.exploration.data_generate_mape import SampleMAPE

logger = logging.getLogger(__name__)

class GenerateData(object):
    """
        Base class for data generation
    """
    def __init__(self, seed_data, **kwargs):
        np.random.seed(seed_data)
        self.initial = True

    def _shuffle(self, *args):
        idx = np.random.permutation(len(args[0]))
        out_args = []
        for a in args:
            if type(a) == np.ndarray:
                out_args.append(a[idx])
            elif type(a) == list:
                out_args.append(list(itemgetter(*idx)(a)))
        return tuple(out_args)

    def generate_data(self, **kwargs):
        raise NotImplementedError

class GenerateDisplacement(GenerateData, SampleMAPE,
                           DisplacementParamSampling, LatentSampling):
    """
        Wrapper around agent environment evaluations.

        Generate data by running the brick-hitting-the-ball environment,
        with a displacement-based policy, or a nn policy.
    """

    def __init__(self, task_object, exploration, experiment_directory,
                       **kwargs):
        super().__init__(**kwargs)
        self.dirname = experiment_directory
        # Unpack task details
        self.task = task_object
        self.num_outcomes = self.task.num_outcomes
        self.parameter_dims = self.task.controller_obj.parameter_dims
        self.parameter_arch = self.task.controller_obj.parameter_arch
        self.traj_dim = self.task.traj_dim

        self.metric_dim = self.task.environment.bd_metric.metric_dim
        self.metric_name = self.task.environment.bd_metric.metric_name
        self.env_info = self.task.environment.env_info
        self.env_name = self.task.env_name
        self.ctrl_name = self.task.ctrl_name
        # Unpack exploration details
        self.stype = exploration['scale_type']
        self.srange = exploration['limit_range']
        self.irange = exploration['init_range'] \
            if exploration['init_range'] is not None else (-1, 1)
        self.exploration_normal = exploration['normal']
        self.exploration_initial = \
                exploration['initial'] if exploration['initial'] is not None \
                                       else self.exploration_normal
        self.num_samples = exploration['num_samples']
        self.k_init = exploration['k_init']
        if '_mape' in exploration['normal']:
            self.mape_dict = {'niter': exploration['mape_niter'],
                              'decay': exploration['mape_decay'],
                              'sigma': exploration['mape_sigma'],
                              'use_fit': int(exploration['mape_use_fitness'])}

    def _eval_population(self, next_param_population,
                         from_latent, recn_fn=None, convert=True,
                         shuffle=False, balance=False, **kwargs):
        """
            Reconstruct samples from ebedded space and evaluate
            the policy parameters
        """
        assert len(next_param_population)>0, "Empty population!"
        if from_latent and recn_fn is not None:
            if len(next_param_population.shape) == 1:
                next_param_population = next_param_population.reshape(
                                            -1, next_param_population.shape[0])
            next_param_population_hat = recn_fn(next_param_population)
        elif convert:
            next_param_population_hat = \
                self.task.controller_obj.convert_params(
                                            uniforms=next_param_population)
        else:
            next_param_population_hat = next_param_population.copy()
        # Evaluate the params
        out_dict = self.task.evaluate_population(
                                  param_population=next_param_population_hat)
        # # Balance data to avoid disproportionate labels
        # if balance:
        #     out_tuple = self.task.balance_outcomes(*out_tuple)
        # if shuffle:
        #     out_tuple = self._shuffle(*out_tuple)
        if from_latent:
            out_dict.update({'param_embedded': next_param_population})
        # else:
        #     out_dict.update({'param_embedded': None})
        return out_dict

    @timing
    def generate_data(self, **kwargs):
        """
            Generate a parameter population using a given exploration_strategy
        """
        if kwargs['nloop']==0:
            self.exploration_type = self.exploration_initial
            _k, _name = self.k_init, 'INITIAL'
            _mape_iter = ''
        else:
            self.exploration_type = self.exploration_normal
            _k, _name = 1, 'NORMAL'
            _mape_iter = '{} iters'.format(self.mape_dict['niter']) \
                            if 'mape' in self.exploration_normal else ''
        logging.info("EXPLORATION START - {} '{}' ({}@{:.0f} samples);".format(
                                          _name, self.exploration_type,
                                          _mape_iter, _k*self.num_samples))

        # Call function to generate data population
        low, high = self.irange
        trial_outputs = getattr(self, self.exploration_type)(k=_k, low=low,
                                                                   high=high,
                                                                   **kwargs)

        # Organise the returned data in a dict for main dataset
        if isinstance(trial_outputs, tuple):
            trial_outputs = dict(param_original=trial_outputs[0],
                                 outcomes=trial_outputs[1],
                                 metric_bd=trial_outputs[2],
                                 traj_main=trial_outputs[3],
                                 traj_aux=trial_outputs[4])
            if len(trial_outputs)==6:
                trial_outputs.update({'param_embedded': trial_outputs[5]})
            else:
                trial_outputs.update({'param_embedded': None})

        # Log exploration progress
        trial_outs = trial_outputs['outcomes']
        logger.info("EXPLORATION DONE '{}' generated {} [{}-dim] new samples:"
                    "\n{}{}".format(self.exploration_type,
                    len(trial_outs) if trial_outs is not None else 0,
                    self.parameter_dims, _TAB, info_outcome(trial_outs)))
        return trial_outputs


text:
import random as rn
import numpy as np

class Triangle(object):
    '''
        Triangle Class

        A triangle is defined by its vertex coordinates,
        background color and alpha component.
    '''

    def __init__(self, *args):

        # Define parameters for triangle generation

        if len(args) == 3:

            #  Explicit definition with parameters: height , width , max_opacity

            self.height = args[0]
            self.width = args[1]
            self.maxopacity = args[2]

            self.points = ((rn.randint(0, self.width), rn.randint(0, self.height)),
                           (rn.randint(0, self.width), rn.randint(0, self.height)),
                           (rn.randint(0, self.width), rn.randint(0, self.height)))

            self.color = (rn.randint(0, 255), rn.randint(0, 255), rn.randint(0, 255))

            self.alpha = self.maxopacity * rn.random()

        elif len(args) == 4:

            # Implicit definition by taking: height , width , max_opacity, decoding gene string

            self.height = args[0]
            self.width = args[1]
            self.maxopacity = args[2]

            info = args[3].split('-')
            self.points = ((int(info[0]), int(info[1])),
                           (int(info[2]), int(info[3])),
                           (int(info[4]), int(info[5])))

            self.color = (int(info[6]), int(info[7]), int(info[8]))
            self.alpha = float(info[9])

    def getInfo(self):

        # Retrive information about the object

        return (self.points, self.color, self.alpha)

    def encodeGene(self):

        # Returns a string encoding the gene information

        return ("{0}-{1}-{2}-{3}-{4}-{5}-{6}-{7}-{8}-{9}\n".format(self.points[0][0], self.points[0][1],
                                                                   self.points[1][0], self.points[1][1],
                                                                   self.points[2][0], self.points[2][1],
                                                                   self.color[0], self.color[1],
                                                                   self.color[2], self.alpha))

    def randomize(self):

        # Produces  variations in gene parameters to a certain degree specified

        new_points_x_0 = int(self.points[0][0] + np.random.normal(0, 30))
        new_points_y_0 = int(self.points[0][1] + np.random.normal(0, 30))

        new_points_x_1 = int(self.points[1][0] + np.random.normal(0, 30))
        new_points_y_1 = int(self.points[1][1] + np.random.normal(0, 30))

        new_points_x_2 = int(self.points[2][0] + np.random.normal(0, 30))
        new_points_y_2 = int(self.points[2][1] + np.random.normal(0, 30))

        new_color_r = int(self.color[0] + np.random.normal(0, 5))
        new_color_g = int(self.color[1] + np.random.normal(0, 5))
        new_color_b = int(self.color[2] + np.random.normal(0, 5))

        new_alpha = self.alpha + np.random.normal(0, 0.5)

        # Update randomized values
        self.points = ( (max(0, min(self.width, new_points_x_0)),
                       max(0, min(self.height, new_points_y_0))),
                       (max(0, min(self.width, new_points_x_1)),
                       max(0, min(self.height, new_points_y_1))),
                       (max(0, min(self.width, new_points_x_2)),
                       max(0, min(self.height, new_points_y_2))))

        self.color = (max(0, min(255, new_color_r)),
                      max(0, min(255, new_color_g)),
                      max(0, min(255, new_color_b)))

        self.alpha = max(0.0, min(self.maxopacity, new_alpha))