Develop step

.gitignore

@@ -68,7 +68,7 @@ tests/GLWACSO/HSP2results/hspp007.uci
 tests/test_report_conversion.html
 
 # Omit big files
-tests/**/*.h5
+tests/**/**/*.h5
 tests/testcbp/HSP2results/*.csv
 
 # R files

examples/state_specl_ops/compare_eq_to_specl.py

@@ -0,0 +1,135 @@
+# if testing manually you may need to os.chdir('./tests/test10specl/HSP2results')
+import os
+import pandas as pd
+import numpy as np
+from pandas import read_hdf
+import h5py
+#from pandas.io.pytables import read_hdf
+#from HSP2IO.hdf import HDF5
+from pathlib import Path
+from src.hsp2.hsp2tools.commands import import_uci, run
+from src.hsp2.hsp2tools.HDF5 import HDF5
+from src.hsp2.hsp2tools.HBNOutput import HBNOutput
+
+def test_h5_file_exists():
+    assert os.path.exists('test10.h5')
+
+
+############# HSPF With SPECIAL ACTION 
+test_root = Path("tests/test10specl")
+test_root.exists()
+test_root_hspf = Path(test_root) / "HSPFresults"
+hspf_data = HBNOutput(os.path.join(test_root_hspf, 'test10speclR.hbn'))
+hspf_data.read_data()
+# get RCHRES 5 sedtrn silt
+ts_silt_hspf = hspf_data.get_time_series('RCHRES', 5, 'RSEDTOTSILT', 'SEDTRN', 'full')
+total_silt_hspf = ts_silt_hspf.mean()
+quantile_silt_hspf = np.quantile(ts_silt_hspf,[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+
+############# HSPF NO SPECL
+hspf_nospecl_root = Path("tests/test10/HSPFresults")
+hspf_nospecl_data = HBNOutput(os.path.join(hspf_nospecl_root, 'test10R.hbn'))
+hspf_nospecl_data.read_data()
+ts_silt_nospecl_hspf = hspf_nospecl_data.get_time_series('RCHRES', 5, 'RSEDTOTSILT', 'SEDTRN', 'full')
+total_silt_nospecl_hspf = ts_silt_nospecl_hspf.mean()
+quantile_silt_nospecl_hspf = np.quantile(ts_silt_nospecl_hspf,[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+
+############# hsp2 SPECL 
+# Run and Analyze hsp2 WITH SPECL actions
+specl_root = Path("tests/test10specl")
+specl_root.exists()
+specl_root_hsp2 = Path(specl_root) / "HSP2results"
+hsp2_specl_uci = specl_root_hsp2.resolve() / "test10specl.uci"
+hsp2_specl_uci.exists()
+temp_specl_h5file = specl_root_hsp2 / "test10specl.h5"
+
+# IF we want to run it from python, do this:
+    # if temp_specl_h5file.exists():
+    #    temp_specl_h5file.unlink()
+    # load the UCI into the h5 then run it
+    # import_uci(str(hsp2_specl_uci), str(temp_specl_h5file))
+    # run(temp_specl_h5file, saveall=True, compress=False)
+# Load Data from hdf5 & Analyze
+dstore_specl = pd.HDFStore(str(temp_specl_h5file), mode='r')
+hsp2_specl_hydr5 = read_hdf(dstore_specl, '/RESULTS/RCHRES_R005/HYDR')
+hsp2_specl_sedtrn5 = read_hdf(dstore_specl, '/RESULTS/RCHRES_R005/SEDTRN')
+hsp2_specl_rsed5 = hsp2_specl_sedtrn5['RSED5']
+quantile_silt_hsp2 = np.quantile(hsp2_specl_rsed5,[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+quantile_ro_hsp2 = np.quantile(hsp2_specl_hydr5['RO'],[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+total_silt_hsp2 = hsp2_specl_rsed5.mean()
+dstore_specl.close()
+
+############# hsp2 w/out SPECL 
+# Run and Analyze hsp2 without SPECL actions
+nospecl_root = Path("tests/test10")
+nospecl_root.exists()
+nospecl_root_hspf = Path(nospecl_root) / "HSPFresults"
+hsp2_nospecl_uci = nospecl_root_hspf.resolve() / "test10.uci"
+hsp2_nospecl_uci.exists()
+temp_nospecl_h5file = nospecl_root_hspf / "nospecl_case.h5"
+
+# IF we want to run it from python, do this:
+    #if temp_nospecl_h5file.exists():
+    #    temp_nospecl_h5file.unlink()
+    # load the UCI into the h5 then run it
+    #import_uci(str(hsp2_nospecl_uci), str(temp_nospecl_h5file))
+    #run(temp_nospecl_h5file, saveall=True, compress=False)
+# Load Data from hdf5 & Analyze
+dstore_nospecl = pd.HDFStore(str(temp_nospecl_h5file), mode='r')
+hsp2_nospecl_hydr5 = read_hdf(dstore_nospecl, '/RESULTS/RCHRES_R005/HYDR')
+hsp2_nospecl_sedtrn5 = read_hdf(dstore_nospecl, '/RESULTS/RCHRES_R005/SEDTRN')
+hsp2_nospecl_rsed5 = hsp2_nospecl_sedtrn5['RSED5']
+np.quantile(hsp2_nospecl_rsed5,[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+np.quantile(hsp2_nospecl_hydr5['RO'],[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+total_silt_nospecl_hsp2 = hsp2_nospecl_rsed5.mean()
+dstore_nospecl.close()
+
+
+############# hsp2 w/equation replicating SPECL 
+# Run and Analyze hsp2 WITH SPECL actions
+eq_root = Path("tests/test10eq")
+eq_root.exists()
+eq_root_hsp2 = Path(eq_root) / "HSP2results"
+hsp2_eq_uci = eq_root_hsp2.resolve() / "test10eq.uci"
+hsp2_eq_uci.exists()
+temp_eq_h5file = eq_root_hsp2 / "test10eq.h5"
+# IF we want to run it from python, do this:
+    #if temp_eq_h5file.exists():
+    #    temp_eq_h5file.unlink()
+
+    # load the UCI into the h5 then run it
+    #import_uci(str(hsp2_eq_uci), str(temp_eq_h5file))
+    #run(temp_eq_h5file, saveall=True, compress=False)
+# Load Data from hdf5 & Analyze
+dstore_eq_specl = pd.HDFStore(str(temp_eq_h5file), mode='r')
+hsp2_eq_hydr5 = read_hdf(dstore_eq_specl, '/RESULTS/RCHRES_R005/HYDR')
+hsp2_eq_sedtrn5 = read_hdf(dstore_eq_specl, '/RESULTS/RCHRES_R005/SEDTRN')
+hsp2_eq_rsed5 = hsp2_eq_sedtrn5['RSED5']
+quantile_silt_eq_hsp2 = np.quantile(hsp2_eq_rsed5,[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+quantile_ro_eq_hsp2 = np.quantile(hsp2_eq_hydr5['RO'],[0.0,0.10,0.5,0.5,0.75,0.9,0.95,1.0])
+total_silt_eq_hsp2 = hsp2_eq_rsed5.mean()
+dstore_eq_specl.close()
+
+# Calculate the % difference
+print("Total Silt ")
+print([
+  {'HSPF no specl': total_silt_nospecl_hspf},
+  {'HSPF specl': total_silt_hspf},
+  {'HSP2 specl': total_silt_hsp2}
+])
+print(
+  {'HSPF no specl': total_silt_nospecl_hspf}
+)
+print([
+  {'HSPF no specl': total_silt_nospecl_hspf},
+  {'HSP2 no specl': total_silt_nospecl_hsp2},
+  {'HSP2 eq': total_silt_eq_hsp2},
+  {'HSP2 specl': total_silt_hsp2}
+])
+print([
+])
+pct_dif_specl = round(100.0 * (total_silt_hsp2.mean() - total_silt_hspf.mean()) / total_silt_hsp2.mean(), 3)
+pct_dif_nospecl = round(100.0 * (total_silt_nospecl_hsp2.mean() - total_silt_nospecl_hspf.mean()) / total_silt_nospecl_hspf.mean(), 3)
+print("Total SiltHSP2 vs. HSPF, % difference = ", pct_dif_specl, "%")
+print("No SPECL: Total SiltHSP2 vs. HSPF, % difference = ", pct_dif_nospecl, "%")
+

examples/state_specl_ops/demo_equation_cli.py

@@ -0,0 +1,75 @@
+# Must be run from the HSPsquared source directory, the h5 file has already been setup with hsp import_uci test10.uci
+# bare bones tester - must be run from the HSPsquared source directory
+import os
+import numpy
+from hsp2.hsp2.main import *
+from hsp2.hsp2.om import *
+from hsp2.hsp2.state import *
+from hsp2.hsp2io.hdf import HDF5
+from hsp2.hsp2io.io import IOManager
+
+fpath = "./tests/test10eq/HSP2results/test10eq.h5"
+# try also:
+# fpath = './tests/testcbp/HSP2results/JL1_6562_6560.h5'
+# sometimes when testing you may need to close the file, so try:
+# f = h5py.File(fpath,'a') # use mode 'a' which allows read, write, modify
+# # f.close()
+hdf5_instance = HDF5(fpath)
+
+
+io_manager = IOManager(hdf5_instance)
+uci_obj = io_manager.read_uci()
+siminfo = uci_obj.siminfo
+opseq = uci_obj.opseq
+# Note: now that the UCI is read in and hdf5 loaded, you can see things like:
+# - hdf5_instance._store.keys() - all the paths in the UCI/hdf5
+# - finally stash specactions in state, not domain (segment) dependent so do it once
+# now load state and the special actions
+state = init_state_dicts()
+state_siminfo_hsp2(uci_obj, siminfo, io_manager, state)
+# Add support for dynamic functions to operate on STATE
+# - Load any dynamic components if present, and store variables on objects
+state_load_dynamics_hsp2(state, io_manager, siminfo)
+# Iterate through all segments and add crucial paths to state
+# before loading dynamic components that may reference them
+state_init_hsp2(state, opseq, activities)
+# - finally stash specactions in state, not domain (segment) dependent so do it once
+state["specactions"] = uci_obj.specactions  # stash the specaction dict in state
+om_init_state(state)  # set up operational model specific state entries
+specl_load_state(state, io_manager, siminfo)  # traditional special actions
+state_load_dynamics_om(
+    state, io_manager, siminfo
+)  # operational model for custom python
+# finalize all dynamically loaded components and prepare to run the model
+state_om_model_run_prep(state, io_manager, siminfo)
+
+# check dependencies 
+# get context first (sets up domain in state) 
+# operation = RCHRES, segment = RCHRES_004, activity = SEDTRN
+state_context_hsp2(state, "RCHRES", "R005", "SEDTRN")
+ep_list = ["RSED4", "RSED5", "RSED6"]
+domain = state['domain']
+model_exec_list = model_domain_dependencies(state, state['domain'], ep_list)
+rsed4 = state['model_object_cache'][domain + "/" + 'RSED4']
+test10eq = state['model_object_cache']['/STATE/test10eq']
+RCHRESR005 = state['model_object_cache']['/STATE/test10eq/RCHRES_R005']
+RCHRESR005.inputs
+hydr_get_ix(state['state_ix'], state['state_paths'], domain)
+hvars = hydr_state_vars()
+get_domain_state(state['state_paths'], state['state_ix'], domain, hvars)
+dep,ivol,o1,o2,o3,ovol1,ovol2,ovol3,prsupy,ro,rovol,sarea,tau,ustar,vol,volev = get_domain_state(state['state_paths'], state['state_ix'], domain, hvars)
+sedvars = sedtrn_state_vars()
+rsed4,rsed5,rsed6 = get_domain_state(state['state_paths'], state['state_ix'], domain, sedvars)
+start = time.time()
+model_domain_dependencies(state, state["domain"], ep_list)
+iterate_models(
+    model_exec_list, state["op_tokens"], state["state_ix"], state["dict_ix"], state["ts_ix"], siminfo["steps"], -1)
+end = time.time()
+print(
+    len(model_exec_list),
+    "components iterated over state_ix",
+    siminfo["steps"],
+    "time steps took",
+    end - start,
+    "seconds",
+)

examples/state_specl_ops/demo_specl_initialize.py

@@ -0,0 +1,36 @@
+
+##########################################################################################
+# LOAD HSP2 RUNTIME CODE AND UCI FILE
+##########################################################################################
+import os, numpy
+from hsp2.hsp2.main import *
+from hsp2.hsp2.om import *
+from hsp2.hsp2.state import *
+from hsp2.hsp2io.hdf import HDF5
+from hsp2.hsp2io.io import IOManager
+hdf5_instance = HDF5("./tests/test10specl/HSP2results/test10specl.demo.h5")
+
+
+iomanager = IOManager(hdf5_instance)
+uci_obj = iomanager.read_uci()
+siminfo = uci_obj.siminfo
+opseq = uci_obj.opseq
+state = init_state_dicts()
+state_siminfo_hsp2(uci_obj, siminfo, iomanager, state)
+
+# Add support for dynamic functions to operate on STATE
+state_load_dynamics_hsp2(state, iomanager, siminfo)
+state_init_hsp2(state, opseq, activities)
+state["specactions"] = uci_obj.specactions  # stash the specaction dict in state
+om_init_state(state)  # set up operational model specific state entries
+specl_load_state(state, iomanager, siminfo)  # traditional special actions
+state_load_dynamics_om(state, iomanager, siminfo) 
+state_om_model_run_prep(state, iomanager, siminfo)
+state_context_hsp2(state, "RCHRES", "R005", "SEDTRN")
+
+domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens = state["domain"], state["state_paths"], state["state_ix"], state["dict_ix"], state["ts_ix"], state["op_tokens"]
+ep_list = np.asarray(["RSED1", "RSED2", "RSED3", "RSED4", "RSED5", "RSED6"], dtype='U')
+model_exec_list = model_domain_dependencies(state, domain, ep_list, True)
+get_domain_state(state_paths, state_ix, domain, ep_list)
+
+

examples/state_specl_ops/demo_step_specl_state.py

@@ -0,0 +1,31 @@
+# MUST RUN CODE FROM demo_specl_iniitialize.py before running this.
+from hsp2.hsp2.sedtrn_step import step_sedtrn
+##########################################################################################
+# SAMPLE MAIN LOOP (does not include timestep changes)
+##########################################################################################
+start = time.time()
+numsteps = siminfo['steps'] * 40
+for step in range(numsteps):
+    #step_hydr(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+    #step_adcalc(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+    #step_cons(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+    #step_htrch(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+    step_sedtrn(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list, step, ep_list)
+    #step_gqual(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+    #step_rqual(domain, state_paths, state_ix, dict_ix, ts_ix, op_tokens, model_exec_list)
+
+end = time.time()
+print(
+    len(model_exec_list), "components iterated over state_ix", numsteps,
+    "time steps took", end - start, "seconds",
+)
+
+state['model_root_object'].get_state("/STATE/test10specl.demo/RCHRES_R005/RSED4")
+state['model_root_object'].get_state("/STATE/test10specl.demo/RCHRES_R005/RSED5")
+state['model_root_object'].get_state("/STATE/test10specl.demo/RCHRES_R005/RSED6")
+state['model_root_object'].inputs
+state['model_object_cache']['/STATE/test10specl.demo/RCHRES_R005'].inputs
+state['model_object_cache']['/STATE/test10specl.demo/RCHRES_R005/RSED6'].inputs
+state['model_object_cache']['/STATE/test10specl.demo/SPECACTION2'].inputs
+# Show all state var paths and indices
+##    print(f"{key}: {value}")

src/hsp2/hsp2/om.py

@@ -599,6 +599,9 @@ def model_order_recursive(
 
 
 def model_input_dependencies(state, exec_list, only_runnable=False):
+    # TODO: is this redundant to model_domain_dependencies?
+    # Cmment in github suggest it is not, and has specific utility
+    # for timeseries values? https://github.com/HARPgroup/HSPsquared/issues/60#issuecomment-2231668979
     mello = exec_list
     mtl = []
     mel = []
@@ -633,9 +636,10 @@ def model_domain_dependencies(state, domain, ep_list, only_runnable=False):
                 endpoint = state["model_object_cache"][domain + "/" + ep]
                 model_order_recursive(endpoint, state["model_object_cache"], mel, mtl)
                 mello = mello + mel
-
+    # TODO: stash the runnable list (mellorun) as a element in dict_ix for cached access during runtime
+    mellorun = ModelObject.runnable_op_list(state["op_tokens"], mello)
     if only_runnable == True:
-        mello = ModelObject.runnable_op_list(state["op_tokens"], mello)
+        mello = mellorun
     return mello
 
 

src/hsp2/hsp2/om_special_action.py

@@ -196,28 +196,31 @@ def hdf5_load_all(hdf_source):
 
 # njit functions for runtime
 
+"""
+# these indices must be adjusted to reflect the number of common op tokens
+# SpecialAction has:
+# - type of condition (+=, -=, ...)
+# - operand 1 (left side)
+# - operand 2 (right side)
+# @tbd: check if time ops have been set and enable/disable accordingly
+#     - 2 ops for each time matching switch: state_ix of the time element (year, month, ...) and the state_ix of the constant to match
+#     - if (state_ix[tix1] <> state_ix[vtix1]): return state_ix[ix1] (don't modify the value)
+#     - alternative: save the integer timestamp or timestep of the start, and if step/stamp > value, enable
+# @tbd: add number of repeats, and save the value of repeats in a register
+"""
+
 
 @njit(cache=True)
 def step_special_action(op, state_ix, dict_ix, step):
     ix = op[1]  # ID of this op
-    # these indices must be adjusted to reflect the number of common op tokens
-    # SpecialAction has:
-    # - type of condition (+=, -=, ...)
-    # - operand 1 (left side)
-    # - operand 2 (right side)
-    # @tbd: check if time ops have been set and enable/disable accordingly
-    #     - 2 ops will be added for each time matching switch, the state_ix of the time element (year, month, ...) and the state_ix of the constant to match
-    #     - matching should be as simple as if (state_ix[tix1] <> state_ix[vtix1]): return state_ix[ix1] (don't modify the value)
-    #     - alternative: save the integer timestamp or timestep of the start, and if step/stamp > value, enable
-    # @tbd: add number of repeats, and save the value of repeats in a register
     ix1 = op[2]  # ID of source of data and destination of data
-    sop = op[3]
-    ix2 = op[4]
-    tix = op[5]  # which slot is the time comparison in?
-    if tix in state_ix and step < state_ix[tix]:
-        return
+    sop = op[3]  # type of conditional comparator
+    ix2 = op[4]  # ID of the other operand
+    tix = op[5]  # which slot is the time to start in?
     ctr_ix = op[6]  # id of the counter variable
     num_ix = op[7]  # max times to complete
+    if tix in state_ix and step < state_ix[tix]:
+        return
     num_done = state_ix[ctr_ix]
     num = state_ix[num_ix]  # num to complete
     if tix in state_ix and num_done >= num:
@@ -234,9 +237,7 @@ def step_special_action(op, state_ix, dict_ix, step):
         elif sop == 5:
             result = state_ix[ix1] / state_ix[ix2]
 
-    # set value in target
-    # tbd: handle this with a model linkage? cons: this makes a loop since the ix1 is source and destination
-
+    # set value in target. tbd: handle this with a model linkage? pros: dependencies. cons: this makes a loop since the ix1 is source and destination
     state_ix[ix1] = result
     state_ix[op[1]] = result
     return result