Merge pull request #103 from HARPgroup/psa

Psa

.gitignore

@@ -63,3 +63,5 @@ tests/GLWACSO/HSPFresults/hspf006.HBN
 tests/GLWACSO/HSP2results/hspp007.uci
 tests/test_report_conversion.html
 tests/land_spec/hwmA51800.h5
+tests/testcbp/HSP2results/PL3_5250_0001.h5
+tests/testcbp/HSP2results/*.csv

HSP2/HYDR.py

@@ -18,6 +18,7 @@
 from numba import njit
 from numba.typed import List
 from HSP2.utilities import initm, make_numba_dict
+from HSP2.SPECL import specl, _specl_
 
 
 ERRMSGS =('HYDR: SOLVE equations are indeterminate',             #ERRMSG0
@@ -30,15 +31,16 @@
 MAXLOOPS  = 100     # newton method exit tolerance
 
 
-def hydr(io_manager, siminfo, uci, ts, ftables):
+def hydr(io_manager, siminfo, uci, ts, ftables, specactions):
     ''' 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)
+    CALL: hydr(store, general, ui, ts, specactions)
        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
-       ts is a dictionary with RID specific timeseries'''
+       ts is a dictionary with RID specific timeseries
+       specactions is a dictionary with all special actions'''
 
     steps   = siminfo['steps']                # number of simulation points
     uunits  = siminfo['units']
@@ -70,11 +72,33 @@ def hydr(io_manager, siminfo, uci, ts, ftables):
 
     # OUTDGT timeseries might come in as OUTDGT, OUTDGT0, etc. otherwise UCI default
     names = list(sorted([n for n in ts if n.startswith('OUTDG')], reverse=True))
+    print(names)
     df = DataFrame()
     for i,c in enumerate(ODGTF):
         df[i] = ts[names.pop()][0:steps] if c > 0 else zeros(steps)
     OUTDGT = df.to_numpy()
 
+    print(OUTDGT)
+
+    # OUTDGT replacement testing
+    # OUTDGT = ts[names][0:steps]
+    xkeys = ['OUTDGT1', 'OUTDGT2']
+    # xvalues = list(map(ts.get, xkeys))
+    xvalues = array(ts['OUTDGT1'],ts['OUTDGT2'])
+    # xvalues = array(map(ts.get, xkeys))
+    # xvalues = list(map(ts.get, names))
+    # print(names)
+    print(xvalues)
+    # print(xvalues.shape)
+    # print(xvalues[0][0])
+    # xvalues[0][0] = xvalues[0][0] * 2
+    # print("ts['OUTDGT1']", ts['OUTDGT1'][0])
+    # print("ts['OUTDGT2']", ts['OUTDGT2'][0])
+
+    # List all names in ts, for jk testing purposes only
+    # ts_names = list(sorted([n for n in ts], reverse=True))
+    # print(ts_names)
+
     # generic SAVE table doesn't know nexits for output flows and rates
     if nexits > 1:
         u = uci['SAVE']
@@ -111,15 +135,21 @@ def hydr(io_manager, siminfo, uci, ts, ftables):
     ui['nodfv']  = any(ODFVF)
     ui['uunits'] = uunits
 
+    # List all names in ui, for jk testing purposes only
+    # ui_names = list(sorted([n for n in ui], reverse=True))
+    # print(ui_names)
+
     # Numba can't do 'O' + str(i) stuff yet, so do it here. Also need new style lists
     Olabels = List()
     OVOLlabels = List()
     for i in range(nexits):
         Olabels.append(f'O{i+1}')
         OVOLlabels.append(f'OVOL{i+1}')
 
+    specactions = make_numba_dict(specactions) # Note: all values coverted to float automatically
     ###########################################################################
-    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels)                  # run reaches simulation code
+    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels, specactions)                  # run reaches simulation code
+#    errors = _hydr_(ui, ts, COLIND, OUTDGT, rchtab, funct, Olabels, OVOLlabels)                  # run reaches simulation code
     ###########################################################################
 
     if 'O'    in ts:  del ts['O']
@@ -134,7 +164,7 @@ def hydr(io_manager, siminfo, uci, ts, ftables):
 
 
 @njit(cache=True)
-def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels):
+def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, specactions):
     errors = zeros(int(ui['errlen'])).astype(int64)
 
     steps  = int(ui['steps'])            # number of simulation steps
@@ -262,12 +292,36 @@ def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels):
 
     # HYDR (except where noted)
     for step in range(steps):
+        # print('\n', 'step: ', step, ' of: ', steps, ' steps')
+
+        # ------------------------------------------------------------------------
+        # print('Trying specl')   
+        # OUTDGT2_save = ts['OUTDGT2'][step - 1] # save before calling specl()
+        # OUTDGT1_save = ts['OUTDGT1'][step - 1] 
+        # print("OUTDGT[step, :]", OUTDGT[step, :])
+        # print("ro", ro)
+
+        # call specl
+        specl(ui, ts, step, specactions)
+        # print("ts['OUTDGT2'][step]", ts['OUTDGT2'][step])
+        # print("ts['OUTDGT1'][step]", ts['OUTDGT1'][step])
+        # print("OUTDGT[step, :]", OUTDGT[step, :])
+
+        # set OUTDGT using the values in the ts object which were set inside specl()
+        OUTDGT[step, :] = [ts['OUTDGT1'][step], ts['OUTDGT2'][step], 0.0]
+        # print("OUTDGT[step, :]", OUTDGT[step, :])
+        # ------------------------------------------------------------------------
+
         convf  = CONVF[step]
         outdgt[:] = OUTDGT[step, :]
         colind[:] = COLIND[step, :]
         roseff = ro
         oseff[:] = o[:]
 
+        # jk test prints
+        # print("roseff", roseff)
+        # print("outdgt[:]", outdgt[:])
+
         # 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
@@ -606,4 +660,5 @@ def expand_HYDR_masslinks(flags, uci, dat, recs):
         rec['TMEMSB2'] = dat.TMEMSB2
         rec['SVOL'] = dat.SVOL
         recs.append(rec)
-    return recs
+    return recs
+

HSP2/SPECL.py

@@ -0,0 +1,45 @@
+''' 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 ts[step - 1]
+    specl_actions is a dictionary with all SPEC-ACTIONS entries
+'''
+
+from numba import njit
+
+@njit
+# def specl(io_manager, siminfo, uci, ts, step, specl_actions):
+def specl(ui, ts, step, specactions):
+
+    # print('Made it to specl()')
+    ts = _specl_(ui, ts, step, specactions)
+
+    # return errors, ERRMSGS
+    # return ts
+
+
+
+# def _specl_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels):
+@njit
+def _specl_(ui, ts, step, specactions):
+
+    # print('Made it to _specl_()')
+    # ts['VOL'][step - 1] = ts['VOL'][step - 1] * 5.0
+    # ts['VOL'][step - 1] = ts['VOL'][step - 1] - specactions['test_wd']
+
+    # ts['OUTDGT'][step - 1] = ts['OUTDGT'][step - 1]  - specactions['test_wd']
+    # ts['OUTDGT2'][step - 1] = 99
+    # ts['OUTDGT2'][step - 1] = ts['OUTDGT2'][step - 1]
+    # ts['OUTDGT2'][step] = ts['OUTDGT2'][step - 1] + 99
+
+    ts['OUTDGT2'][step] = 99
+    # ts['OUTDGT2'][step, :] = 99 # this resulted in errors
+
+    # print(specactions['outdgt'])
+    # return errors
+    # return ts
+

HSP2/main.py

@@ -35,7 +35,7 @@ def main(io_manager:IOManager, saveall:bool=False, jupyterlab:bool=True) -> None
         raise FileNotFoundError(f'{hdfname} HDF5 File Not Found')
 
     msg = messages()
-    msg(1, f'Processing started for file {hdfname}; saveall={saveall}')
+    msg(1, f'(specl v0.1) Processing started for file {hdfname}; saveall={saveall}')
 
     # read user control, parameters, states, and flags uci and map to local variables
     uci_obj = io_manager.read_uci()
@@ -49,7 +49,8 @@ def main(io_manager:IOManager, saveall:bool=False, jupyterlab:bool=True) -> None
     ftables = uci_obj.ftables
     specactions = uci_obj.specactions
     monthdata = uci_obj.monthdata
-
+    specactions = {} # placeholder till added to uci parser
+
     start, stop = siminfo['start'], siminfo['stop']
 
     copy_instances = {}
@@ -208,7 +209,7 @@ def main(io_manager:IOManager, saveall:bool=False, jupyterlab:bool=True) -> None
                 ############ calls activity function like snow() ##############
                 if operation not in ['COPY','GENER']:
                     if (activity == 'HYDR'):
-                        errors, errmessages = function(io_manager, siminfo, ui, ts, ftables)
+                        errors, errmessages = function(io_manager, siminfo, ui, ts, ftables, specactions)
                     elif (activity != 'RQUAL'):
                         errors, errmessages = function(io_manager, siminfo, ui, ts)
                     else:                    

HSP2tools/readUCI.py

@@ -464,7 +464,7 @@ def specactions(info, llines):
     #     column to match with conditional expression if applicable
     #   - The value of CURLVL matches with an expression
     sa_actions = [] # referred to as "classic" in old HSPF code comments 
-    head_actions = ['OPERATION','RANGE1','RANGE2','DC','DS','YR','MO','DA','HR','MN','D','T','VARI', 'S1','S2','AC','VALUE','TC','TS','NUM', 'CURLVL']
+    head_actions = ['OPTYP','RANGE1','RANGE2','DC','DS','YR','MO','DA','HR','MN','D','T','VARI', 'S1','S2','AC','VALUE','TC','TS','NUM', 'CURLVL']
     sa_mult = []
     head_mult = []
     sa_uvquan = []