From 98139e973fdd4ce3a668c15c884e3af6c6ea979e Mon Sep 17 00:00:00 2001 From: George Simmons Date: Thu, 18 Apr 2024 11:58:30 +0100 Subject: [PATCH] allows computation of summaries, function updates --- .gitignore | 5 ++ psymple/custom_functions.py | 123 +++++++++++++++++++++++++++++++++++- psymple/globals.py | 26 +++++++- psymple/ported_objects.py | 2 +- psymple/read_wx.py | 88 ++++++++++++++++++++++++++ psymple/system.py | 21 ++++-- psymple/variables.py | 8 ++- 7 files changed, 258 insertions(+), 15 deletions(-) create mode 100644 psymple/read_wx.py diff --git a/.gitignore b/.gitignore index 40def25..bb1a6f6 100644 --- a/.gitignore +++ b/.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 diff --git a/psymple/custom_functions.py b/psymple/custom_functions.py index dba9dfc..0fdabc7 100644 --- a/psymple/custom_functions.py +++ b/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) \ No newline at end of file + 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 \ No newline at end of file diff --git a/psymple/globals.py b/psymple/globals.py index 0792191..a63f1eb 100644 --- a/psymple/globals.py +++ b/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") diff --git a/psymple/ported_objects.py b/psymple/ported_objects.py index 86ae440..f49146a 100644 --- a/psymple/ported_objects.py +++ b/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} diff --git a/psymple/read_wx.py b/psymple/read_wx.py new file mode 100644 index 0000000..ba99a92 --- /dev/null +++ b/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) + + + diff --git a/psymple/system.py b/psymple/system.py index 37b9b1e..071fb88 100644 --- a/psymple/system.py +++ b/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: diff --git a/psymple/variables.py b/psymple/variables.py index aa707a8..e33d501 100644 --- a/psymple/variables.py +++ b/psymple/variables.py @@ -3,6 +3,7 @@ from typing import List # Deprecated since Python 3.9 import sympy as sym +import numpy as np from psymple.abstract import Container, DependencyError, SymbolWrapper from psymple.globals import T, sym_custom_ns @@ -49,7 +50,7 @@ def update_buffer(self): self.buffer = self.time_series[-1] def update_time_series(self, time_step): - new_value = self.update_rule.evaluate_update(self.buffer, time_step) + new_value = np.maximum(0.0,self.update_rule.evaluate_update(self.buffer, time_step)) self.time_series.append(new_value) @@ -100,6 +101,7 @@ def __init__(self, parameter, computed_value=None): # because it's implicit in the UpdateRule. self.computed_value = computed_value self.update_rule = None + self.time_series = [computed_value] def initialize_update_rule(self, variables, parameters): self.update_rule = SimUpdateRule( @@ -125,6 +127,7 @@ def substitute_parameters(self, variables): def update_value(self): self.computed_value = self.update_rule.evaluate_expression() + self.time_series.append(self.computed_value) class Parameters(Container): @@ -236,12 +239,13 @@ def _lambdify(self): self._equation_lambdified = sym.lambdify( self.variables.get_symbols() + self.parameters.get_symbols(), self.equation, - modules=[sym_custom_ns, "scipy", "numpy"], + modules=[sym_custom_ns, "scipy", "numpy"], cse=True, ) def evaluate_expression(self): if self._equation_lambdified is None: self._lambdify() + print(f"{self.variable.name} lambdified") v_args = [v.buffer for v in self.variables] p_args = [p.computed_value for p in self.parameters] args = v_args + p_args