Develop hdf5

environment.yml

@@ -10,7 +10,7 @@ dependencies:
   # Running HSP2
   - scipy  # Scipy also installs numpy
   # Pandas installs most scientific Python modules, such as Numpy, etc.
-  - pandas
+  - pandas <3.0.0
   - numba
   - numpy
   - hdf5

environment_dev.yml

@@ -12,7 +12,7 @@ dependencies:
   # Running HSP2
   - scipy  # Scipy also installs numpy
   # Pandas installs most scientific Python modules, such as Numpy, etc.
-  - pandas >=2.0
+  - pandas >=2.0, <3.0.0
   - numba
   - numpy
   - hdf5

examples/state_specl_ops/compare_eq_to_specl.py

@@ -1,4 +1,5 @@
 # if testing manually you may need to os.chdir('./tests/test10specl/HSP2results')
+# todo: make sure thsi is path independent.
 import os
 import pandas as pd
 import numpy as np

examples/state_specl_ops/update_pandas.py

@@ -0,0 +1,99 @@
+
+from numpy import float64, ones
+from pandas import DataFrame, date_range
+from pandas.tseries.offsets import Minute
+from datetime import datetime as dt
+from typing import Union
+import os
+# must import Category, Modelfirst to prevent circular error with HDF5
+# this does not happen in real runtime...
+from hsp2.hsp2.model import Model
+from hsp2.hsp2io.protocols import Category
+from hsp2.hsp2io.hdf import HDF5
+from hsp2.hsp2.utilities import (
+    versions,
+    get_timeseries,
+    expand_timeseries_names,
+    save_timeseries,
+    get_gener_timeseries,
+    hoursval
+)
+from hsp2.hsp2.configuration import activities, noop, expand_masslinks
+from hsp2.state.state import (
+    init_state_dicts,
+    state_siminfo_hsp2,
+    state_load_dynamics_hsp2,
+    state_init_hsp2,
+    state_context_hsp2,
+)
+from hsp2.hsp2.om import (
+    om_init_state,
+    state_om_model_run_prep,
+    state_load_dynamics_om,
+    state_om_model_run_finish,
+)
+from hsp2.hsp2.SPECL import specl_load_state
+
+from hsp2.hsp2io.io import IOManager, SupportsReadTS, Category
+
+
+fpath = "./tests/test10/HSP2results/test10.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)
+
+# Begin code from main.py
+
+# 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 = {}
+
+#######################################################################################
+# initialize STATE dicts
+#######################################################################################
+# Set up Things in state that will be used in all modular activities like SPECL
+state = init_state_dicts()
+state_siminfo_hsp2(parameter_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"] = 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)
+#######################################################################################
+
+perland_ext = ddext_sources[('PERLND', 'P001')]
+perlnd_ts = get_timeseries(io_manager, perland_ext, siminfo)
+precip = perlnd_ts['PREC']
+# go another way and get the raw TS to compare to the perlnd_ts which has been converted to simulation
+# intervals
+precip_df = io_manager.read_ts(
+    category=Category.INPUTS, segment='TS039'
+)

pyproject.toml

@@ -14,7 +14,7 @@ dependencies = [
     "cltoolbox",
     "numba",
     "numpy<2.0",
-    "pandas",
+    "pandas<3.0.0",
     "tables",
     "pyparsing"
 ]

src/hsp2/hsp2/utilities.py

@@ -14,7 +14,7 @@
 from numba import types
 from numba.typed import Dict
 from numpy import float64, full, tile, zeros
-from pandas import Series, date_range, Timedelta
+from pandas import Series, date_range, Timedelta, to_timedelta
 from pandas.tseries.offsets import Minute
 
 from hsp2.hsp2io.protocols import Category, SupportsReadTS, SupportsWriteTS
@@ -213,35 +213,35 @@ def transform(ts, name, how, siminfo):
     NOTE: these routines work for both regular and sparse timeseries input
     """
 
-    tsfreq = Timedelta("1 " + ts.index.freqstr)
+    tsfreq = ts.index.freq
     fmins = Minute(siminfo["delt"])
-    freq = Timedelta(fmins).to_timedelta64()
+    freq = fmins.nanos
     stop = siminfo["stop"]
 
     # append duplicate of last point to force processing last full interval
     if ts.index[-1] < stop:
         ts[stop] = ts.iloc[-1]
 
-    if freq == tsfreq:
+    if freq == tsfreq.nanos:
         pass
     elif tsfreq is None:  # Sparse time base, frequency not defined
         ts = ts.reindex(siminfo["tbase"]).ffill().bfill()
     elif how == "SAME":
-        ts = ts.resample(fmins).ffill()  # tsfreq >= freq assumed, or bad user choice
+        ts = ts.resample(fmins).ffill()  # tsfreq.nanos >= freq assumed, or bad user choice
     elif not how:
         if name in flowtype:
-            if "Y" in str(tsfreq) or "M" in str(tsfreq) or tsfreq > freq:
+            if "Y" in str(tsfreq) or "M" in str(tsfreq) or tsfreq.nanos > freq:
                 if "M" in str(tsfreq):
                     ratio = 1.0 / 730.5
                 elif "Y" in str(tsfreq):
                     ratio = 1.0 / 8766.0
                 else:
-                    ratio = freq / tsfreq
+                    ratio = freq / tsfreq.nanos
                 ts = (ratio * ts).resample(fmins).ffill()  # HSP2 how = div
             else:
                 ts = ts.resample(fmins).sum()
         else:
-            if "Y" in str(tsfreq) or "M" in str(tsfreq) or tsfreq > freq:
+            if "Y" in str(tsfreq) or "M" in str(tsfreq) or tsfreq.nanos > freq:
                 ts = ts.resample(fmins).ffiio_managerll()
             else:
                 ts = ts.resample(fmins).mean()
@@ -268,10 +268,10 @@ def transform(ts, name, how, siminfo):
             elif "Y" in str(tsfreq):
                 ratio = 1.0 / (8766.0 * mult)
             else:
-                ratio = freq / tsfreq
+                ratio = freq / tsfreq.nanos
             ts = (ratio * ts).resample(fmins).ffill()  # HSP2 how = div
         else:
-            ts = (ts * (freq / tsfreq)).resample(fmins).ffill()
+            ts = (ts * (freq / tsfreq.nanos)).resample(fmins).ffill()
     elif how == "ZEROFILL":
         ts = ts.resample(fmins).fillna(0.0)
     elif how == "INTERPOLATE":
@@ -289,7 +289,7 @@ def hoursval(siminfo, hours24, dofirst=False, lapselike=False):
     start = siminfo["start"]
     stop = siminfo["stop"]
     fmins = Minute(siminfo["delt"])
-    freq = Timedelta(fmins).to_timedelta64()
+    freq = fmins.nanos
 
     dr = date_range(
         start=f"{start.year}-01-01", end=f"{stop.year}-12-31", freq=Minute(60)
@@ -299,16 +299,16 @@ def hoursval(siminfo, hours24, dofirst=False, lapselike=False):
         hours[0] = 1
 
     ts = Series(hours[0 : len(dr)], dr)
-    tsfreq = Timedelta("1 " + ts.index.freqstr)
+    tsfreq = ts.index.freq
     if lapselike:
-        if tsfreq > freq:  # upsample
+        if tsfreq.nanos > freq:  # upsample
             ts = ts.resample(fmins).asfreq().ffill()
-        elif tsfreq < freq:  # downsample
+        elif tsfreq.nanos < freq:  # downsample
             ts = ts.resample(fmins).mean()
     else:
-        if tsfreq > freq:  # upsample
+        if tsfreq.nanos > freq:  # upsample
             ts = ts.resample(fmins).asfreq().fillna(0.0)
-        elif tsfreq < freq:  # downsample
+        elif tsfreq.nanos < freq:  # downsample
             ts = ts.resample(fmins).max()
     return ts.truncate(start, stop).to_numpy()
 
@@ -325,16 +325,16 @@ def monthval(siminfo, monthly):
     start = siminfo["start"]
     stop = siminfo["stop"]
     fmins = Minute(siminfo["delt"])
-    freq = Timedelta(fmins).to_timedelta64()
+    freq = fmins.nanos
 
     months = tile(monthly, stop.year - start.year + 1).astype(float)
     dr = date_range(start=f"{start.year}-01-01", end=f"{stop.year}-12-31", freq="MS")
     ts = Series(months, index=dr).resample("D").ffill()
-    tsfreq = Timedelta("1 " + ts.index.freqstr)
+    tsfreq = ts.index.freq
 
-    if tsfreq > freq:  # upsample
+    if tsfreq.nanos > freq:  # upsample
         ts = ts.resample(fmins).asfreq().ffill()
-    elif tsfreq < freq:  # downsample
+    elif tsfreq.nanos < freq:  # downsample
         ts = ts.resample(fmins).mean()
     return ts.truncate(start, stop).to_numpy()
 
@@ -345,17 +345,17 @@ def dayval(siminfo, monthly):
     start = siminfo["start"]
     stop = siminfo["stop"]
     fmins = Minute(siminfo["delt"])
-    freq = Timedelta(fmins).to_timedelta64()
+    freq = fmins.nanos
 
     months = tile(monthly, stop.year - start.year + 1).astype(float)
     dr = date_range(start=f"{start.year}-01-01", end=f"{stop.year}-12-31", freq="MS")
     ts = Series(months, index=dr).resample("D").interpolate("time")
-    tsfreq = Timedelta("1 " + ts.index.freqstr)
+    tsfreq = ts.index.freq
 
 
-    if tsfreq > freq:  # upsample
+    if tsfreq.nanos > freq:  # upsample
         ts = ts.resample(fmins).ffill()
-    elif tsfreq < freq:  # downsample
+    elif tsfreq.nanos < freq:  # downsample
         ts = ts.resample(fmins).mean()
     return ts.truncate(start, stop).to_numpy()
 

src/hsp2/hsp2io/hdf.py

@@ -13,9 +13,12 @@ def __init__(self, file_path: str) -> None:
         self._store = pd.HDFStore(file_path)
         None
 
-    def __del__(self):
+    def close(self):
         self._store.close()
 
+    def __del__(self):
+        self.close()
+
     def __enter__(self):
         return self
 

src/hsp2/hsp2tools/commands.py

@@ -24,6 +24,7 @@ def run(h5file, saveall=True, compress=True):
     hdf5_instance = HDF5(h5file)
     io_manager = IOManager(hdf5_instance)
     main(io_manager, saveall=saveall, jupyterlab=compress)
+    hdf5_instance.close()
 
 
 def import_uci(ucifile, h5file):

tests/cmd_regression.py

@@ -0,0 +1,88 @@
+# Note: This code must be run from the tests dir due to importing
+# the `convert` directory where the regression_base file lives
+# todo: make this convert path passable by argument 
+import numpy as np
+from convert.regression_base import RegressTest
+
+case = "test10specl"
+base_case = "test10"
+#case = "testcbp"
+tdir = "/opt/model/HSPsquared/tests"
+#import_uci(str(hsp2_specl_uci), str(temp_specl_h5file))
+#run(temp_specl_h5file, saveall=True, compress=False)
+
+############# RegressTest data loader
+## TEST CASE
+test = RegressTest(case, threads=1)
+# test object hydr
+rchres_hydr_hsp2_test_table = test.hsp2_data._read_table('RCHRES', '001', 'HYDR')
+perlnd_pwat_hsp2_test_table = test.hsp2_data._read_table('PERLND', '001', 'PWATER')
+rchres_hydr_hsp2_test = test.hsp2_data.get_time_series('RCHRES', '001', 'RO', 'HYDR')
+rchres_hydr_hspf_test = test.get_hspf_time_series( ('RCHRES', 'HYDR', '001', 'RO', 2)) #Note: the order of arguments is wonky in Regressbase
+rchres_hydr_hsp2_test_mo = rchres_hydr_hsp2_test.resample('MS').mean()
+rchres_hydr_hspf_test_mo = rchres_hydr_hspf_test.resample('MS').mean()
+## BASE CASE
+base = RegressTest(base_case, threads=1)
+# base object hydr
+rchres_hydr_hsp2_base_table = base.hsp2_data._read_table('RCHRES', '001', 'HYDR')
+perlnd_pwat_hsp2_base_table = base.hsp2_data._read_table('PERLND', '001', 'PWATER')
+rchres_hydr_hsp2_base = base.hsp2_data.get_time_series('RCHRES', '001', 'RO', 'HYDR')
+rchres_hydr_hspf_base = base.get_hspf_time_series( ('RCHRES', 'HYDR', '001', 'RO', 2))  #Note: the order of arguments is wonky in Regressbase
+rchres_hydr_hsp2_base_mo = rchres_hydr_hsp2_base.resample('MS').mean()
+rchres_hydr_hspf_base_mo = rchres_hydr_hspf_base.resample('MS').mean()
+
+# Show quantiles
+print("hsp2", case, np.quantile(rchres_hydr_hsp2_test, [0,0.25,0.5,0.75,1.0]))
+print("hspf", case, np.quantile(rchres_hydr_hspf_test, [0,0.25,0.5,0.75,1.0]))
+print("hsp2", base_case, np.quantile(rchres_hydr_hsp2_base, [0,0.25,0.5,0.75,1.0]))
+print("hspf", base_case, np.quantile(rchres_hydr_hspf_base, [0,0.25,0.5,0.75,1.0]))
+# Monthly mean value comparisons
+rchres_hydr_hsp2_test_mo
+rchres_hydr_hspf_test_mo
+rchres_hydr_hsp2_base_mo
+rchres_hydr_hspf_base_mo
+
+# Compare ANY arbitrary timeseries, not just the ones coded into the RegressTest object
+# 3rd argument is tolerance to use
+tol = 10.0
+test.compare_time_series(rchres_hydr_hsp2_base_table['RO'], rchres_hydr_hsp2_test_table['RO'], tol)
+# Example:  (True, 8.7855425)
+
+# Compare inflows and outflows
+np.mean(perlnd_pwat_hsp2_test_table['PERO']) * 6000 * 0.0833
+np.mean(rchres_hydr_hsp2_test_table['IVOL'])
+np.mean(rchres_hydr_hsp2_test_table['ROVOL'])
+np.mean(perlnd_pwat_hsp2_base_table['PERO']) * 6000 * 0.0833
+np.mean(rchres_hydr_hsp2_base_table['IVOL'])
+np.mean(rchres_hydr_hsp2_base_table['ROVOL'])
+
+# now do a comparison
+# HYDR diff should be almost nonexistent
+test.check_con(params = ('RCHRES', 'HYDR', '001', 'ROVOL', 2))
+# this is very large for PWTGAS
+# git: test10specl ('PERLND', 'PWTGAS', '001', 'POHT', '2') 1163640%
+test.check_con(params = ('PERLND', 'PWTGAS', '001', 'POHT', '2'))
+# Other mismatches in PERLND
+test.check_con(params = ('PERLND', 'PWATER', '001', 'AGWS', '2'))
+test.check_con(params = ('PERLND', 'PWATER', '001', 'PERO', '2'))
+# Now run the full test
+test.quiet = True # this lets us test without overwhelming the console
+results = test.run_test()
+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
+
+print(mismatches)
+
+if mismatches:
+    for case, key, results in mismatches:
+        diff = results
+        print(case, key, f"{diff:0.00%}")
+else:
+    print("No mismatches found. Success!")