Merge pull request #161 from HARPgroup/develop-spec-lean

Develop spec lean

.gitignore

@@ -68,10 +68,8 @@ tests/GLWACSO/HSP2results/hspp007.uci
 tests/test_report_conversion.html
 
 # Omit big files
-tests/land_spec/hwmA51800.h5
-tests/testcbp/HSP2results/PL3_5250_0001.h5
+tests/**/*.h5
 tests/testcbp/HSP2results/*.csv
-tests/test10/HSP2results/test10.h5
 
 # R files
 .Rdata

HSP2/HYDR.py

@@ -12,14 +12,17 @@
 '''
 
 
-from numpy import zeros, any, full, nan, array, int64, arange
+from numpy import zeros, any, full, nan, array, int64, arange, asarray
 from pandas import DataFrame
 from math import sqrt, log10
-from numba import njit
+from numba import njit, types
 from numba.typed import List
 from HSP2.utilities import initm, make_numba_dict
-from HSP2.state import *
-from HSP2.SPECL import specl
+
+# the following imports added by rb to handle dynamic code and special actions
+from HSP2.state import hydr_get_ix, hydr_init_ix
+from HSP2.om import pre_step_model, step_model, model_domain_dependencies
+from numba.typed import Dict
 
 
 ERRMSGS =('HYDR: SOLVE equations are indeterminate',             #ERRMSG0
@@ -36,7 +39,7 @@ def hydr(io_manager, siminfo, uci, ts, ftables, state):
     ''' find the state of the reach/reservoir at the end of the time interval
     and the outflows during the interval
 
-    CALL: hydr(store, general, ui, ts, specactions)
+    CALL: hydr(store, general, ui, ts, state)
        store is the Pandas/PyTable open store
        general is a dictionary with simulation level infor (OP_SEQUENCE for example)
        ui is a dictionary with RID specific HSPF UCI like data
@@ -122,34 +125,36 @@ def hydr(io_manager, siminfo, uci, ts, ftables, state):
     for i in range(nexits):
         Olabels.append(f'O{i+1}')
         OVOLlabels.append(f'OVOL{i+1}')
-
-    # state_info is some generic things about the simulation 
+
+    #######################################################################################
+    # the following section (1 of 3) added to HYDR by rb to handle dynamic code and special actions
+    #######################################################################################
+    # state_info is some generic things about the simulation
+    # must be numba safe, so we don't just pass the whole state which is not
     state_info = Dict.empty(key_type=types.unicode_type, value_type=types.unicode_type)
     state_info['operation'], state_info['segment'], state_info['activity'] = state['operation'], state['segment'], state['activity']
-    state_info['domain'], state_info['state_step_hydr'] = state['domain'], state['state_step_hydr']
+    state_info['domain'], state_info['state_step_hydr'], state_info['state_step_om'] = state['domain'], state['state_step_hydr'], state['state_step_om']
     hsp2_local_py = state['hsp2_local_py']
     # It appears necessary to load this here, instead of from main.py, otherwise,
     # _hydr_() does not recognize the function state_step_hydr()? 
     if (hsp2_local_py != False):
         from hsp2_local_py import state_step_hydr
     else:
         from HSP2.state_fn_defaults import state_step_hydr
+    # initialize the hydr paths in case they don't already reside here
+    hydr_init_ix(state, state['domain'])
     # must split dicts out of state Dict since numba cannot handle mixed-type nested Dicts
     state_ix, dict_ix, ts_ix = state['state_ix'], state['dict_ix'], state['ts_ix']
     state_paths = state['state_paths']
-    # initialize the hydr paths in case they don't already reside here
-    hydr_init_ix(state_ix, state_paths, state['domain'])
-
-    ###########################################################################
-    # specactions - special actions code TBD
-    ###########################################################################
-    specactions = make_numba_dict(state['specactions']) # Note: all values coverted to float automatically
-
-    ###########################################################################
-    # Do the simulation with _hydr_()
-    ###########################################################################
-    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, specactions, state_step_hydr) # run reaches simulation code
-    ###########################################################################
+    ep_list = ["DEP","IVOL","O1","O2","O3","OVOL1","OVOL2","OVOL3","PRSUPY","RO","ROVOL","SAREA","TAU","USTAR","VOL","VOLEV"]
+    model_exec_list = model_domain_dependencies(state, state_info['domain'], ep_list)
+    model_exec_list = asarray(model_exec_list, dtype="i8") # format for use in numba
+    op_tokens = state['op_tokens']
+    #######################################################################################
+
+    # Do the simulation with _hydr_   (ie run reaches simulation code)
+    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels,
+                    state_info, state_paths, state_ix, dict_ix, ts_ix, state_step_hydr, op_tokens, model_exec_list)
 
     if 'O'    in ts:  del ts['O']
     if 'OVOL' in ts:  del ts['OVOL']
@@ -163,7 +168,7 @@ def hydr(io_manager, siminfo, uci, ts, ftables, state):
 
 
 @njit(cache=True)
-def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, specactions, state_step_hydr):
+def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, state_step_hydr, op_tokens, model_exec_list):
     errors = zeros(int(ui['errlen'])).astype(int64)
 
     steps  = int(ui['steps'])            # number of simulation steps
@@ -292,8 +297,11 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_inf
     # other initial vars
     rovol = 0.0
     volev = 0.0
-    IVOL0   = ts['IVOL']                   # the actual inflow in simulation native units 
-    # prepare for dynamic state
+    IVOL0   = ts['IVOL']                   # the actual inflow in simulation native units
+
+    #######################################################################################
+    # the following section (2 of 3) added by rb to HYDR, this one to prepare for dynamic state including special actions
+    #######################################################################################
     hydr_ix = hydr_get_ix(state_ix, state_paths, state_info['domain'])
     # these are integer placeholders faster than calling the array look each timestep
     o1_ix, o2_ix, o3_ix, ivol_ix = hydr_ix['O1'], hydr_ix['O2'], hydr_ix['O3'], hydr_ix['IVOL']
@@ -306,15 +314,19 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_inf
         out_ix[1] = o2_ix
     if nexits > 2:
         out_ix[2] = o3_ix
+    #######################################################################################
+
     # HYDR (except where noted)
     for step in range(steps):
-        # call specl
-        specl(ui, ts, step, state_info, state_paths, state_ix, specactions)
         convf  = CONVF[step]
         outdgt[:] = OUTDGT[step, :]
         colind[:] = COLIND[step, :]
         roseff = ro
         oseff[:] = o[:]
+
+        #######################################################################################
+        # the following section (3 of 3) added by rb to accommodate dynamic code, operations models, and special actions
+        #######################################################################################
         # set state_ix with value of local state variables and/or needed vars
         # Note: we pass IVOL0, not IVOL here since IVOL has been converted to different units
         state_ix[ro_ix], state_ix[rovol_ix] = ro, rovol
@@ -323,17 +335,30 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_inf
             state_ix[out_ix[oi]] = outdgt[oi] 
         state_ix[vol_ix], state_ix[ivol_ix] = vol, IVOL0[step]
         state_ix[volev_ix] = volev
-        # Execute dynamic code if enabled
+        # - these if statements may be irrelevant if default functions simply return
+        #   when no objects are defined.
+        if (state_info['state_step_om'] == 'enabled'):
+            pre_step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step)
         if (state_info['state_step_hydr'] == 'enabled'):
             state_step_hydr(state_info, state_paths, state_ix, dict_ix, ts_ix, hydr_ix, step)
+        if (state_info['state_step_om'] == 'enabled'):
+            #print("trying to execute state_step_om()")
+            # model_exec_list contains the model exec list in dependency order
+            # now these are all executed at once, but we need to make them only for domain end points
+            step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step)   # traditional 'ACTIONS' done in here
+        if ( (state_info['state_step_hydr'] == 'enabled')
+            or (state_info['state_step_om'] == 'enabled') ):
             # Do write-backs for editable STATE variables
             # OUTDGT is writeable
             for oi in range(nexits):
                 outdgt[oi] = state_ix[out_ix[oi]]
-            # IVOL is writeable. 
-            # Note: we must convert IVOL to the units expected in _hydr_ 
+            # IVOL is writeable.
+            # Note: we must convert IVOL to the units expected in _hydr_
             # maybe routines should do this, and this is not needed (but pass VFACT in state)
             IVOL[step] = state_ix[ivol_ix] * VFACT
+        # End dynamic code step()
+        #######################################################################################
+
         # vols, sas variables and their initializations  not needed.
         if irexit >= 0:             # irrigation exit is set, zero based number
             if rirwdl > 0.0:  # equivalent to OVOL for the irrigation exit

HSP2/SEDTRN.py

@@ -3,12 +3,17 @@
 License: LGPL2
 '''
 
-from numpy import array, zeros, where, int64
+from numpy import array, zeros, where, int64, asarray
 from math import log10, exp
-from numba import njit
+from numba import njit, types
 from HSP2.ADCALC import advect
 from HSP2.utilities  import make_numba_dict
 
+# the following imports added to handle special actions
+from HSP2.state import sedtrn_get_ix, sedtrn_init_ix
+from HSP2.om import pre_step_model, step_model, model_domain_dependencies
+from numba.typed import Dict
+
 ERRMSGS =('SEDTRN: Warning -- bed storage of sediment size fraction sand is empty',                                   #ERRMSG0
           'SEDTRN: Warning -- bed storage of sediment size fraction silt is empty',                                   #ERRMSG1
           'SEDTRN: Warning -- bed storage of sediment size fraction clay is empty',                                   #ERRMSG2
@@ -17,7 +22,7 @@
           'SEDTRN: Simulation of sediment requires all 3 "auxiliary flags" (AUX1FG, etc) in section HYDR must be turned on', #ERRMSG5
           'SEDTRN: When specifying the initial composition of the bed, the fraction of sand, silt, and clay must sum to a value close to 1.0.')  #ERRMSG6
 
-def sedtrn(io_manager, siminfo, uci, ts):
+def sedtrn(io_manager, siminfo, uci, ts, state):
 	''' Simulate behavior of inorganic sediment'''
 
 	# simlen = siminfo['steps']
@@ -68,8 +73,35 @@ def sedtrn(io_manager, siminfo, uci, ts):
 	ui['clay_taucs'] = ui_clay['TAUCS']
 	ui['clay_m']     = ui_clay['M']	* delt60 / 24.0 * 4.880 # convert erodibility coeff from /day to /ivl
 
+	#######################################################################################
+	# the following section (1 of 3) added to SEDTRN by pbd to handle special actions
+	#######################################################################################
+	# state_info is some generic things about the simulation
+	# must be numba safe, so we don't just pass the whole state which is not
+	state_info = Dict.empty(key_type=types.unicode_type, value_type=types.unicode_type)
+	state_info['operation'], state_info['segment'], state_info['activity'] = state['operation'], state['segment'], state['activity']
+	state_info['domain'], state_info['state_step_hydr'], state_info['state_step_om'] = state['domain'], state['state_step_hydr'], state['state_step_om']
+	# hsp2_local_py = state['hsp2_local_py']
+	# # It appears necessary to load this here, instead of from main.py, otherwise,
+	# # _hydr_() does not recognize the function state_step_hydr()?
+	# if (hsp2_local_py != False):
+	# 	from hsp2_local_py import state_step_hydr
+	# else:
+	# 	from HSP2.state_fn_defaults import state_step_hydr
+	# must split dicts out of state Dict since numba cannot handle mixed-type nested Dicts
+	# initialize the sedtrn paths in case they don't already reside here
+	sedtrn_init_ix(state, state['domain'])
+	state_ix, dict_ix, ts_ix = state['state_ix'], state['dict_ix'], state['ts_ix']
+	state_paths = state['state_paths']
+	op_tokens = state['op_tokens']
+	# Aggregate the list of all SEDTRN end point dependencies
+	ep_list = ['RSED4', 'RSED5', 'RSED6']
+	model_exec_list = model_domain_dependencies(state, state_info['domain'], ep_list)
+	model_exec_list = asarray(model_exec_list, dtype="i8") # format for use in 
+	#######################################################################################
+
 	############################################################################
-	errors = _sedtrn_(ui, ts)  # run SEDTRN simulation code
+	errors = _sedtrn_(ui, ts, state_info, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)  # run SEDTRN simulation code
 	############################################################################
 
 	if nexits > 1:
@@ -91,7 +123,7 @@ def sedtrn(io_manager, siminfo, uci, ts):
 	return errors, ERRMSGS
 
 @njit(cache=True)
-def _sedtrn_(ui, ts):
+def _sedtrn_(ui, ts, state_info, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list):
 	''' Simulate behavior of inorganic sediment'''
 	errorsV = zeros(int(ui['errlen'])).astype(int64)
 
@@ -291,8 +323,36 @@ def _sedtrn_(ui, ts):
 
 	#################### END PSED
 
+	#######################################################################################
+	# the following section (2 of 3) added by pbd to SEDTRN, this one to prepare for special actions
+	#######################################################################################
+	sedtrn_ix = sedtrn_get_ix(state_ix, state_paths, state_info['domain'])
+	# these are integer placeholders faster than calling the array look each timestep
+	rsed4_ix, rsed5_ix, rsed6_ix = sedtrn_ix['RSED4'], sedtrn_ix['RSED5'], sedtrn_ix['RSED6']
+	#######################################################################################
+
 	for loop in range(simlen):
 
+		#######################################################################################
+		# the following section (3 of 3) added by pbd to accommodate special actions
+		#######################################################################################
+		# set state_ix with value of local state variables and/or needed vars
+		state_ix[rsed4_ix] = sand_wt_rsed4
+		state_ix[rsed5_ix] = silt_wt_rsed5
+		state_ix[rsed6_ix] = clay_wt_rsed6
+		if (state_info['state_step_om'] == 'enabled'):
+			pre_step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step = loop)
+
+		# (todo) Insert code hook for dynamic python modification of state 
+
+		if (state_info['state_step_om'] == 'enabled'):
+			step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step = loop)  # traditional 'ACTIONS' done in here
+			# Do write-backs for editable STATE variables
+			sand_wt_rsed4 = state_ix[rsed4_ix]
+			silt_wt_rsed5 = state_ix[rsed5_ix]
+			clay_wt_rsed6 = state_ix[rsed6_ix]
+		#######################################################################################
+
 		# perform any necessary unit conversions
 		if uunits == 2:  # uci is in metric units
 			avvele = AVVEL[loop] * 3.28

HSP2/SPECL.py

@@ -1,18 +1,57 @@
-''' process special actions in this domain
-
-CALL: specl(io_manager, siminfo, uci, ts, state, specl_actions)
-    store is the Pandas/PyTable open store
-    siminfo is a dictionary with simulation level infor (OP_SEQUENCE for example)
-    ui is a dictionary with RID specific HSPF UCI like data
-    ts is a dictionary with RID specific timeseries
-    state is a dictionary with value of variables at ts[step - 1]
-    specl_actions is a dictionary with all SPEC-ACTIONS entries
-'''
-
-from numba import njit
-
-@njit
-def specl(ui, ts, step, state_info, state_paths, state_ix, specactions):
-    # ther eis no need for _specl_ because this code must already be njit
-    return
+''' process special actions in this domain
+Notes:
+  - code for parsing UCI SPEC-ACTIONS is in HSP2tools/readUCI.py
+  - code for object classes that transform parsed data into OP codes for OM and STATE support
+    is in this directory tree as om_special_[action type].py,
+    - Ex: om_special_action.py contains object support and runtime functions for classic ACTIONS
+'''
+
+from numba import njit
+from pandas import DataFrame, date_range
+
+def specl_load_actions(state, io_manager, siminfo):
+    if 'ACTIONS' in state['specactions']:
+        dc = state['specactions']['ACTIONS']
+        for ix in dc.index:
+            # add the items to the state['model_data'] dict
+            speca = dc[ix:(ix+1)]
+            # need to add a name attribute
+            opname = 'SPEC' + 'ACTION' + str(ix)
+            state['model_data'][opname] = {}
+            state['model_data'][opname]['name'] = opname
+            for ik in speca.keys():
+                #print("looking for speca key ", ik)
+                state['model_data'][opname][ik] = speca.to_dict()[ik][ix]   # add subscripts?
+                if ik == 'VARI':
+                    if len(speca.to_dict()['S1'][ix]) > 0:
+                        state['model_data'][opname][ik] += speca.to_dict()['S1'][ix]
+                    if len(speca.to_dict()['S2'][ix]) > 0:
+                        state['model_data'][opname][ik] += speca.to_dict()['S2'][ix]
+            state['model_data'][opname]['object_class'] = 'SpecialAction'
+            #print("model_data", ix, " = ", state['model_data'][opname])
+    return
+
+def specl_load_state(state, io_manager, siminfo):
+    specl_load_actions(state, io_manager, siminfo)
+    # others defined below, like:
+    # specl_load_uvnames(state, io_manager, siminfo)
+    # ...
+    return
+
+'''
+# the code specl() is deprecated in favor of execution inside OM
+# see om_special_action.py for example of object support and runtime functions for classic ACTIONS
+CALL: specl(ui, ts, step, state_info, state_paths, state_ix, specactions)
+    store is the Pandas/PyTable open store
+    siminfo is a dictionary with simulation level infor (OP_SEQUENCE for example)
+    ui is a dictionary with RID specific HSPF UCI like data
+    ts is a dictionary with RID specific timeseries
+    state is a dictionary with value of variables at ts[step - 1]
+    specl_actions is a dictionary with all SPEC-ACTIONS entries
+'''
+
+@njit
+def specl(ui, ts, step, state_info, state_paths, state_ix, specactions):
+    # ther eis no need for _specl_ because this code must already be njit
+    return