Develop state class

examples/state_specl_ops/compare_eq_to_specl.py

@@ -10,6 +10,7 @@
 from src.hsp2.hsp2tools.commands import import_uci, run
 from src.hsp2.hsp2tools.HDF5 import HDF5
 from src.hsp2.hsp2tools.HBNOutput import HBNOutput
+import tabulate
 
 def test_h5_file_exists():
     assert os.path.exists('test10.h5')
@@ -24,15 +25,15 @@ def test_h5_file_exists():
 # 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])
+quantile_silt_hspf = np.quantile(ts_silt_hspf,[0.0,0.25,0.5,0.75,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])
+quantile_silt_nospecl_hspf = np.quantile(ts_silt_nospecl_hspf,[0.0,0.25,0.5,0.75,1.0])
 
 ############# hsp2 SPECL 
 # Run and Analyze hsp2 WITH SPECL actions
@@ -54,19 +55,19 @@ def test_h5_file_exists():
 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])
+quantile_silt_hsp2 = np.quantile(hsp2_specl_rsed5,[0.0,0.25,0.5,0.75,1.0])
+quantile_ro_hsp2 = np.quantile(hsp2_specl_hydr5['RO'],[0.0,0.25,0.5,0.75,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"
+nospecl_root_hsp2 = Path(nospecl_root) / "HSP2results"
+hsp2_nospecl_uci = nospecl_root_hsp2.resolve() / "test10.uci"
 hsp2_nospecl_uci.exists()
-temp_nospecl_h5file = nospecl_root_hspf / "nospecl_case.h5"
+temp_nospecl_h5file = nospecl_root_hsp2 / "test10.h5"
 
 # IF we want to run it from python, do this:
     #if temp_nospecl_h5file.exists():
@@ -79,8 +80,8 @@ def test_h5_file_exists():
 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])
+quantile_silt_nospecl_hsp2 = np.quantile(hsp2_nospecl_rsed5,[0.0,0.25,0.5,0.75,1.0])
+np.quantile(hsp2_nospecl_hydr5['RO'],[0.0,0.25,0.5,0.75,1.0])
 total_silt_nospecl_hsp2 = hsp2_nospecl_rsed5.mean()
 dstore_nospecl.close()
 
@@ -105,8 +106,8 @@ def test_h5_file_exists():
 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])
+quantile_silt_eq_hsp2 = np.quantile(hsp2_eq_rsed5,[0.0,0.25,0.5,0.75,1.0])
+quantile_ro_eq_hsp2 = np.quantile(hsp2_eq_hydr5['RO'],[0.0,0.25,0.5,0.75,1.0])
 total_silt_eq_hsp2 = hsp2_eq_rsed5.mean()
 dstore_eq_specl.close()
 
@@ -133,3 +134,19 @@ def test_h5_file_exists():
 print("Total SiltHSP2 vs. HSPF, % difference = ", pct_dif_specl, "%")
 print("No SPECL: Total SiltHSP2 vs. HSPF, % difference = ", pct_dif_nospecl, "%")
 
+
+a = pd.DataFrame(
+    {
+        'HSP2 no specl':quantile_silt_nospecl_hsp2,
+        'HSPF no specl':quantile_silt_nospecl_hspf,
+        'HSPF w/specl':quantile_silt_hspf,
+        'HSP2 w/specl':quantile_silt_hsp2,
+        'HSP2 w/EQ':quantile_silt_eq_hsp2
+    }
+)
+
+
+tabulate(a, headers='keys', tablefmt='markdown')
+tabulate(a, headers='keys', tablefmt='psql')
+# Thi is better, it USES the tabulate lib but has usable format
+a.to_markdown(index=False)

examples/state_specl_ops/debug_pyt_est.py

@@ -0,0 +1,77 @@
+from pathlib import Path
+
+import pytest
+from hsp2.hsp2tools.commands import import_uci, run
+from hsp2.hsp2tools.HDF5 import HDF5
+
+from convert.regression_base import RegressTest as RegressTestBase
+
+
+class RegressTest(RegressTestBase):
+    def _get_hsp2_data(self, test_root) -> None:
+        test_root_hspf = Path(test_root) / "HSPFresults"
+        hspf_uci = test_root_hspf.resolve() / f"{self.compare_case}.uci"
+        assert hspf_uci.exists()
+
+        temp_h5file = test_root_hspf / f"{self.compare_case}.h5"
+        if temp_h5file.exists():
+            temp_h5file.unlink()
+
+        self.temp_h5file = temp_h5file
+
+        import_uci(str(hspf_uci), str(self.temp_h5file))
+        run(self.temp_h5file, saveall=True, compress=False)
+        self.hsp2_data = HDF5(str(self.temp_h5file))
+
+    def _init_files(self):
+        test_dir = Path(__file__).resolve().parent
+        assert test_dir.name == "tests"
+
+        test_root = test_dir / self.compare_case
+        assert test_root.exists()
+
+        self._get_hbn_data(str(test_root))
+        self._get_hsp2_data(str(test_root))
+
+
+@pytest.mark.parametrize(
+    "case",
+    [
+        # "test05",
+        # "test09",
+        "test10",
+        "test10specl",
+        # "test10b",
+    ],
+)
+def test_case(case):
+    test = RegressTest(case, threads=1)
+    results = test.run_test()
+    test.temp_h5file.unlink()
+
+    found = False
+    mismatches = []
+    for key, results in results.items():
+        no_data_hsp2, no_data_hspf, match, diff = results
+        if any([no_data_hsp2, no_data_hspf]):
+            continue
+        if not match:
+            mismatches.append((case, key, results))
+        found = True
+    assert found
+
+    if mismatches:
+        for case, key, results in mismatches:
+            _, _, _, diff = results
+            print(case, key, f"{diff:0.00%}")
+        raise ValueError("results don't match hspf output")
+
+
+# demo the code in regression_base.py with this:
+# creates a shell object namd "self" that can have attributes just added
+self = type('', (), {})()
+self.compare_case='test10'
+self.html_file = os.path.join(
+    tests_root_dir, f"HSPF_HSP2_{self.compare_case}.html"
+)
+self.html_file

examples/state_specl_ops/demo_land_specl.py

@@ -0,0 +1,38 @@
+# 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
+
+from hsp2.hsp2.main import *
+from hsp2.hsp2.om import *
+from hsp2.hsp2.SPECL import *
+from hsp2.hsp2io.hdf import HDF5
+from hsp2.hsp2io.io import IOManager
+from hsp2.hsp2tools.readUCI import *
+
+from hsp2.state.state import *
+
+fpath = "./tests/land_spec/HSP2Results/hwmA51800.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)
+
+parameter_obj = io_manager.read_parameters()
+opseq = parameter_obj.opseq
+ddlinks = parameter_obj.ddlinks
+ddmasslinks = parameter_obj.ddmasslinks
+ddext_sources = parameter_obj.ddext_sources
+ddgener = parameter_obj.ddgener
+model = parameter_obj.model
+siminfo = parameter_obj.siminfo
+ftables = parameter_obj.ftables
+specactions = parameter_obj.specactions
+monthdata = parameter_obj.monthdata
+conditions = specactions["conditions"]
+actions = specactions["ACTIONS"]
+
+start, stop = siminfo["start"], siminfo["stop"]
+hdf5_instance._store.close()

examples/state_specl_ops/demo_specl_initialize.py

@@ -1,36 +1,90 @@
+# 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
 
-##########################################################################################
-# LOAD HSP2 RUNTIME CODE AND UCI FILE
-##########################################################################################
-import os, numpy
+import os
+import numpy
 from hsp2.hsp2.main import *
+from hsp2.state.state import *
 from hsp2.hsp2.om import *
-from hsp2.hsp2.state import *
+from hsp2.hsp2.SPECL import *
 from hsp2.hsp2io.hdf import HDF5
 from hsp2.hsp2io.io import IOManager
-hdf5_instance = HDF5("./tests/test10specl/HSP2results/test10specl.demo.h5")
+from hsp2.state.state import *
+from hsp2.hsp2.om_timer import timer_class
 
+fpath = "./tests/test10specl/HSP2results/test10specl.h5"
+timer = timer_class()
+# try also:
+# fpath = './tests/testcbp/HSP2results/JL1_6562_6560.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)
+# 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)
+
+parameter_obj = io_manager.read_parameters()
+opseq = parameter_obj.opseq
+ddlinks = parameter_obj.ddlinks
+ddmasslinks = parameter_obj.ddmasslinks
+ddext_sources = parameter_obj.ddext_sources
+ddgener = parameter_obj.ddgener
+model = parameter_obj.model
+siminfo = parameter_obj.siminfo
+ftables = parameter_obj.ftables
+specactions = parameter_obj.specactions
+monthdata = parameter_obj.monthdata
 
-# 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)
+# read user control, parameters, states, and flags parameters and map to local variables
+parameter_obj = io_manager.read_parameters()
+opseq = parameter_obj.opseq
+ddlinks = parameter_obj.ddlinks
+ddmasslinks = parameter_obj.ddmasslinks
+ddext_sources = parameter_obj.ddext_sources
+ddgener = parameter_obj.ddgener
+model = parameter_obj.model
+siminfo = parameter_obj.siminfo
+ftables = parameter_obj.ftables
+specactions = parameter_obj.specactions
+monthdata = parameter_obj.monthdata
 
+start, stop = siminfo["start"], siminfo["stop"]
 
+copy_instances = {}
+gener_instances = {}
+section_timing = {}
+
+section_timing["io_manager.read_parameters() call and config"] = str(timer.split()) + "seconds"
+#######################################################################################
+# initialize STATE dicts
+#######################################################################################
+# Set up Things in state that will be used in all modular activities like SPECL
+state = state_class(
+    state_empty["state_ix"], state_empty["op_tokens"], state_empty["state_paths"], 
+    state_empty["op_exec_lists"], state_empty["model_exec_list"], state_empty["dict_ix"], 
+    state_empty["ts_ix"], state_empty["hsp_segments"]
+)
+om_operations = om_init_state()  # set up operational model specific containers
+state_siminfo_hsp2(state, parameter_obj, siminfo, io_manager)
+state_om_model_root_object(state, om_operations, 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, timer)
+om_init_hsp2_segments(state, om_operations)
+# now initialize all state variables for mutable variables
+hsp2_domain_dependencies(state, opseq, activities, om_operations, False)
+# 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)
+# - finally stash specactions in state, not domain (segment) dependent so do it once
+specl_load_om(om_operations, specactions)  # load traditional special actions
+state_load_dynamics_om(
+    state, io_manager, siminfo, om_operations
+)  # operational model for custom python
+# finalize all dynamically loaded components and prepare to run the model
+state_om_model_run_prep(opseq, activities, state, om_operations, siminfo)
+section_timing["state om initialization()"] = str(timer.split()) + "seconds"
+statenb = state_class_lite(0)
+state_copy(state, statenb)
+#######################################################################################