allows computation of summaries, function updates

Author: GJHSimmons

.gitignore

@@ -26,6 +26,9 @@ share/python-wheels/
 *.egg
 MANIFEST
 
+# Examples WIP
+examples/pbdm/
+
 # PyInstaller
 #  Usually these files are written by a python script from a template
 #  before PyInstaller builds the exe, so as to inject date/other infos into it.
@@ -159,3 +162,5 @@ cython_debug/
 #  option (not recommended) you can uncomment the following to ignore the entire idea folder.
 #.idea/
 .vscode/settings.json
+examples/nesting.py
+examples/nesting2.py

psymple/custom_functions.py

@@ -5,6 +5,7 @@
 from scipy.optimize import root_scalar
 from sympy.abc import x, a, b
 from sympy import integrate, sympify
+from psymple import read_wx
 
 T = sym.Symbol("T")
 
@@ -26,7 +27,7 @@
     [19.4, 6.7],
     [19.3, 6.7],
     [20.4, 8.3],
-]
+] 
 
 
 def DegreeDays_fun(
@@ -108,6 +109,7 @@ def FFTemperature_fun(
         return np.max((0.0, 1 - ((temp - temp_min - diff) / diff) ** 2))
 
 
+
 class FFTemperature(sym.Function):
     #TODO: this could have a doit() method
     @classmethod
@@ -116,4 +118,121 @@ def eval(cls, model_day, temp_min, temp_max):
             return 0
 
     def _eval_evalf(self, prec):
-        return sym.Float(FFTemperature_fun(*self.args))._eval_evalf(prec)
+        return sym.Float(FFTemperature_fun(*self.args))._eval_evalf(prec)
+    
+
+wx = read_wx.Weather.readwx("c:\\Users\\georg\\OneDrive\\Documents\\IDEMS\\CASAS Global\\Modelling\\Python\\population-modeling-py-3\\examples\\pbdm\\sample_weather.txt","01","01","2000","12","31","2010", [0])[1:]
+wx = np.array(wx)
+
+def temp_max_fun(
+    time: float,
+    temperature_data = wx
+):
+    model_day = math.ceil(time)
+    return temperature_data[model_day][5]
+
+def temp_min_fun(
+    time: float,
+    temperature_data = wx
+):
+    model_day = math.ceil(time)
+    return temperature_data[model_day][6]
+
+def solar_rad_fun(
+    time: float,
+    temperature_data = wx
+):
+    model_day = math.ceil(time)
+    return temperature_data[model_day][7]
+
+def temp_fun(
+    time: float,
+):
+    temp_max = temp_max_fun(time)
+    temp_min = temp_min_fun(time)
+    #return (temp_max - temp_min)/2 * np.sin(np.float64(2*np.pi*(time+1/4))) + (temp_max + temp_min)/2
+    return (temp_max + temp_min)/2
+
+
+class temp_max(sym.Function):
+    #TODO: this could have a doit() method
+    @classmethod
+    def eval(cls, time):
+        if isinstance(time, sym.Float) and time < 0:
+            return 0
+
+    def _eval_evalf(self, prec):
+        return sym.Float(temp_max_fun(*self.args))._eval_evalf(prec)
+
+class temp_min(sym.Function):
+    #TODO: this could have a doit() method
+    @classmethod
+    def eval(cls, time):
+        if isinstance(time, sym.Float) and time < 0:
+            return 0
+
+    def _eval_evalf(self, prec):
+        return sym.Float(temp_min_fun(*self.args))._eval_evalf(prec)
+    
+class solar_rad(sym.Function):
+    #TODO: this could have a doit() method
+    @classmethod
+    def eval(cls, time):
+        if isinstance(time, sym.Float) and time < 0:
+            return 0
+
+    def _eval_evalf(self, prec):
+        return sym.Float(solar_rad_fun(*self.args))._eval_evalf(prec)
+
+class temp(sym.Function):
+    #TODO: this could have a doit() method
+    @classmethod
+    def eval(cls, time):
+        if isinstance(time, sym.Float) and time < 0:
+            return 0
+
+    def _eval_evalf(self, prec):
+        return sym.Float(temp_fun(*self.args))._eval_evalf(prec)
+    
+def DD_fun(time, Del, T_min, T_max = 0, coeff_1 = 0, coeff_2 = 0):
+    T = temp_fun(time)
+    if T_max < T_min:
+        return np.maximum(0.01,T-T_min)
+    else:
+        return np.maximum(0, Del*(coeff_1*(T-T_min)/(1+ coeff_2**(T - T_max))))
+
+class DD(sym.Function):
+    #TODO: this could have a doit() method
+    @classmethod
+    def eval(cls, time, Del, T_min, T_max = None, coeff_1 = None, coeff_2 = None):
+        if isinstance(time, sym.Float) and time < 0:
+            return 0
+
+    def doit(self, deep=False, **hints):
+        time, *args = self.args
+        if deep:
+            time = time.doit(deep=deep, **hints)
+        return sym.Float(DD_fun(time, *args))
+
+    def _eval_evalf(self, prec):
+        return self.doit(deep=False)._eval_evalf(prec)
+        #return sym.Float(DD_fun(*self.args))._eval_evalf(prec) 
+        #T = temp_fun(time)
+        #return sym.Float(DD_fun(time=self.args[0], Del = self.args[1], T_min = self.args[2], T_max = self.args[3], coeff_1 = self.args[4], coeff_2 = self.args[5]))._eval_evalf(prec)
+
+def ind_above_fun(base, comp):
+    return 1 if comp>base else 0
+
+class ind_above(sym.Function):
+    @classmethod
+    def eval(cls, base, comp):
+        return ind_above_fun(base, comp)
+
+def frac0_fun(numerator, denominator, default):
+    return default if denominator == 0 else numerator/denominator
+
+class frac0(sym.Function):
+    @classmethod
+    def eval(cls, numerator, denominator, default):
+        if denominator == 0:
+            return default

psymple/globals.py

@@ -1,13 +1,33 @@
 import sympy as sym
-from psymple.custom_functions import DegreeDays, FFTemperature
+from psymple.custom_functions import (
+    DegreeDays,
+    FFTemperature,
+    temp,
+    DD,
+    temp_fun,
+    DD_fun,
+    solar_rad_fun,
+    ind_above_fun,
+    frac0_fun,
+    temp_min_fun,
+)
 
 # from custom_functions import DegreeDays, FFTemperature
 
 
 # The dictionary passed to sym.sympify and sym.lambdify to convert custom functions
 # TODO: This should not be a global property.
-#sym_custom_ns = {}
-sym_custom_ns = {'DegreeDays': DegreeDays, 'FFTemperature': FFTemperature}
+# sym_custom_ns = {}
+sym_custom_ns = {
+    "DegreeDays": DegreeDays,
+    "FFTemperature": FFTemperature,
+    "temp": temp_fun,
+    "temp_min": temp_min_fun,
+    "DD": DD_fun,
+    "solar_rad": solar_rad_fun,
+    "ind_above": ind_above_fun,
+    "frac0": frac0_fun,
+}
 
 
 T = sym.Symbol("T")

psymple/ported_objects.py

@@ -61,7 +61,7 @@ def substitute_symbol(self, old_symbol, new_symbol):
     #     ]
     #     self.equation = self.equation.subs(substitutions)
 
-    def get_free_symbols(self, global_symbols=set()):
+    def get_free_symbols(self, global_symbols=set([sym.Symbol('T')])):
         assignment_symbols = self.expression.free_symbols
         return assignment_symbols - global_symbols - {self.symbol_wrapper.symbol}
 

psymple/read_wx.py

@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Feb  2 23:00:56 2021
+
+@author: Jose Ricardo Cure
+"""
+class Weather(object):
+    def __init__(self):
+        super().__init__()
+
+    def filterWX(wx,startM,startD,startY,endM,endD,endY,weatherfile):
+        '''locate dates to run in weather file and call main'''
+        start=False
+        end=False
+        cont=3
+
+        keyStart= str(startM) +' '+ str(startD)+' ' + str(startY)
+        keyEnd= str(endM)+' ' + str(endD)+' ' + str(endY)
+        print(keyStart, keyEnd)
+        for sList in (wx):
+            if cont < len(wx):
+                #cont  = cont+1
+                keyWX = str(wx[cont][0])+' ' +  str(wx[cont][1])+' ' + str(wx[cont][2])
+                if keyWX == keyStart:
+                    if start==False:
+                        index_start = cont
+                        start=True
+                        lat   = float(wx[1][0])
+                        long  = float(wx[1][1])
+                        month = int(wx[cont][0])
+                        day   = int(wx[cont][1])
+                        year  = int(wx[cont][2])
+                        tmax  = float(wx[cont][3])
+                        tmin  = float(wx[cont][4])
+                        solar = float(wx[cont][5])
+                        precip= float(wx[cont][6])
+                        rh    = float(wx[cont][7])
+                        wind  = float(wx[cont][8])
+                        #print('   first date to run '+ str(month) +' ' + str(day) + ' ' + str(year))
+                        weatherfile.append([lat,long,month,day,year, tmax,tmin,solar,precip,rh,wind])
+
+                elif keyWX == keyEnd:
+                    if end == False:
+                        index_finish = cont+1
+                        totdays = index_finish - cont
+                        keyWX = str(wx[cont][0])+' ' +  str(wx[cont][1])+' ' + str(wx[cont][2])
+                        end=True
+                        #print('   last date to run '+ str(month) +' ' + str(day) + ' ' + str(year)+ '  ')
+                        #print( '.. total simulation days in this file = '+str(totdays) )
+
+                cont=cont+1
+                if end==False and start == True:
+                    lat   = float(wx[1][0])
+                    long  = float(wx[1][1])
+                    month = int(wx[cont][0])
+                    day   = int(wx[cont][1])
+                    year  = int(wx[cont][2])
+                    tmax  = float(wx[cont][3])
+                    tmin  = float(wx[cont][4])
+                    solar = float(wx[cont][5])
+                    precip= float(wx[cont][6])
+                    rh    = float(wx[cont][7])
+                    wind  = float(wx[cont][8])
+                    weatherfile.append([lat,long,month,day,year, tmax,tmin,solar,precip,rh,wind])
+
+
+
+    def readwx(wx_dir,startM,startD,startY,endM,endD,endY,weatherfile):
+        '''read weather file'''
+        with open (wx_dir , 'r') as f: # use with to open your files, it close them automatically
+            wx = [x.split() for x in f]
+            wxfile= wx[0] # title of the wx file
+            lat   = float(wx[1][0])
+            long  = float(wx[1][1])
+            print()
+            print('weather file ', wxfile)
+            print(' latitude '+ str(lat), ' longitude '+ str(long))
+            print('.. original number of days in WX file = ' + str(len(wx)))
+            Weather.filterWX(wx,startM,startD,startY,endM,endD,endY,weatherfile)
+            return weatherfile
+
+
+
+#in_locations = 'd:/fruit_flies/r_USA_MEX_observed_1980-2010_AgMERRA_coarse_test.txt'
+#weatherfile = Weather.run_locations(in_locations)
+
+
+

psymple/system.py

@@ -79,7 +79,7 @@ def _create_variables(self, variables):
         for variable in variables:
             if variable.initial_value is None:
                 print(f"Warning: Variable {variable.symbol} has no initial value")
-                variable.initial_value = 0
+                variable.initial_value = 0.0
         return Variables([SimVariable(variable) for variable in variables])
 
     def _create_parameters(self, parameters):
@@ -92,6 +92,7 @@ def _create_parameters(self, parameters):
 
     def _assign_update_rules(self, update_rules):
         combined_update_rules = update_rules._combine_update_rules()
+        print(combined_update_rules)
         for rule in combined_update_rules:
             new_rule = SimUpdateRule.from_update_rule(
                 rule, self.variables + self.time, self.parameters
@@ -165,34 +166,37 @@ def _wrap_for_solve_ivp(self, t, y):
             for v in self.variables
         ]
 
-    def _advance_time(self, time_step):
+    def _advance_time(self, time_step, summaries):
         self.time.update_buffer()
         for variable in self.variables:
             variable.update_buffer()
         for variable in self.variables:
             variable.update_time_series(time_step)
+        for parameter in summaries:
+            parameter.update_value()
         self.time.update_time_series(time_step)
 
-    def _advance_time_unit(self, n_steps):
+    def _advance_time_unit(self, n_steps, summaries):
         if n_steps <= 0 or not isinstance(n_steps, int):
             raise ValueError(
                 "Number of time steps in a day must be a positive integer, "
                 f"not '{n_steps}'."
             )
         for i in range(n_steps):
-            self._advance_time(1 / n_steps)
+            self._advance_time(1 / n_steps, summaries)
 
-    def simulate(self, t_end, n_steps, mode="dscr"):
+    def simulate(self, t_end, n_steps, mode="dscr", summaries = []):
         if t_end <= 0 or not isinstance(t_end, int):
             raise ValueError(
                 "Simulation time must terminate at a positive integer, "
                 f"not '{t_end}'."
             )
         self._compute_substitutions()
+        summaries = [self.parameters[s] for s in summaries]
         if mode == "discrete" or mode == "dscr":
             print("dscr")
             for i in range(t_end):
-                self._advance_time_unit(n_steps)
+                self._advance_time_unit(n_steps, summaries)
         elif mode == "continuous" or mode == "cts":
             print("cts")
             sol = solve_ivp(
@@ -220,7 +224,10 @@ def plot_solution(self, variables, t_range=None):
             variables = {v: {} for v in variables}
         legend = []
         for var_name, options in variables.items():
-            variable = self.variables[var_name]
+            try: 
+                variable = self.variables[var_name]
+            except:
+                variable = self.parameters[var_name]
             if isinstance(options, str):
                 plt.plot(t_series[sl], variable.time_series[sl], options)
             else: