diff --git a/.gitignore b/.gitignore index 8cad0c9b..7a22225a 100644 --- a/.gitignore +++ b/.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 diff --git a/examples/state_specl_ops/compare_eq_to_specl.py b/examples/state_specl_ops/compare_eq_to_specl.py new file mode 100644 index 00000000..5e21d9b2 --- /dev/null +++ b/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, "%") + diff --git a/examples/state_specl_ops/demo_equation_cli.py b/examples/state_specl_ops/demo_equation_cli.py new file mode 100644 index 00000000..bbc1e910 --- /dev/null +++ b/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", +) diff --git a/examples/state_specl_ops/demo_specl_initialize.py b/examples/state_specl_ops/demo_specl_initialize.py new file mode 100644 index 00000000..2f6b3fbb --- /dev/null +++ b/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) + + diff --git a/examples/state_specl_ops/demo_step_specl_state.py b/examples/state_specl_ops/demo_step_specl_state.py new file mode 100644 index 00000000..38abc465 --- /dev/null +++ b/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}") diff --git a/src/hsp2/hsp2/om.py b/src/hsp2/hsp2/om.py index 31654aef..5cf5eeb2 100644 --- a/src/hsp2/hsp2/om.py +++ b/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 diff --git a/src/hsp2/hsp2/om_special_action.py b/src/hsp2/hsp2/om_special_action.py index 7e03a847..57e1ee01 100644 --- a/src/hsp2/hsp2/om_special_action.py +++ b/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 diff --git a/src/hsp2/hsp2/sedtrn_step.py b/src/hsp2/hsp2/sedtrn_step.py new file mode 100644 index 00000000..69bdb86b --- /dev/null +++ b/src/hsp2/hsp2/sedtrn_step.py @@ -0,0 +1,59 @@ +from numpy import array, zeros, where, int64, asarray +from math import log10, exp +from numba import njit, types +from numba.types import List + +# the following imports added to handle special actions +from hsp2.hsp2.state import ( + sedtrn_get_ix, + sedtrn_init_ix, + get_domain_state, + set_domain_state, +) +from hsp2.hsp2.om import pre_step_model, step_model, model_domain_dependencies +from numba.typed import Dict + + +@njit +def step_sedtrn( + domain, + state_paths, + state_ix, + dict_ix, + ts_ix, + op_tokens, + model_exec_list, + step, + ep_list, +): + # model_exec_list: a list of elements (specl etc.) that influence these SEDTRN end points + # ep_list = np.asarray(["RSED1", "RSED2", "RSED3", "RSED4", "RSED5", "RSED6"], dtype='U') + # NOTE: this could be cached in dict_ix + # call related specl/ops pre-steps, such as loading timeseries values + pre_step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step) + # call related specl/ops steps + step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step) + # get state value at beginning of timestep - python experts will no doubt have a more code efficient method than this + sand_rsed1, silt_rsed2, clay_rsed3, sand_wt_rsed4, silt_wt_rsed5, clay_wt_rsed6 = ( + get_domain_state(state_paths, state_ix, domain, ep_list) + ) + + # now, do sedtrn (simplified for demo purposes) + tsed1 = sand_rsed1 + silt_rsed2 + clay_rsed3 + tsed2 = sand_wt_rsed4 + silt_wt_rsed5 + clay_wt_rsed6 + sand_t_rsed7 = sand_rsed1 + sand_wt_rsed4 + silt_t_rsed8 = silt_rsed2 + silt_wt_rsed5 + clay_t_rsed9 = clay_rsed3 + clay_wt_rsed6 + tsed3 = sand_t_rsed7 + silt_t_rsed8 + clay_t_rsed9 + + # pass values back to state + state_vals = [ + sand_rsed1, + silt_rsed2, + clay_rsed3, + sand_wt_rsed4, + silt_wt_rsed5, + clay_wt_rsed6, + ] + set_domain_state(state_paths, state_ix, domain, ep_list, state_vals) + return diff --git a/src/hsp2/hsp2/state.py b/src/hsp2/hsp2/state.py index d35e6dd8..5d15b8f2 100644 --- a/src/hsp2/hsp2/state.py +++ b/src/hsp2/hsp2/state.py @@ -188,6 +188,42 @@ def state_load_dynamics_hsp2(state, io_manager, siminfo): state["hsp2_local_py"] = hsp2_local_py # Stores the actual function in state +@njit +def get_domain_state(state_paths, state_ix, domain, varkeys): + # get values for a set of variables in a domain + # will not check for the index in state_ix, and will fail if a non-scalar value is needed (like from dict_ix) + # if varkeys = False, assume that we want all the variables + # from the domain, that are predetermined ahead of time, and should save performance + ret_vals = np.zeros(len(varkeys)) + j = 0 + for i in varkeys: + # var_path = f'{domain}/{i}' + var_path = domain + "/" + i + # print(var_path) + ix = state_paths[var_path] + # print("ix",ix) + ret_vals[j] = state_ix[ix] + j += 1 + return ret_vals + + +@njit +def set_domain_state(state_paths, state_ix, domain, varkeys, state_vals): + # get values for a set of variables in a domain + # will not check for the index in state_ix, and will fail if a non-scalar value is needed (like from dict_ix) + # if varkeys = False, assume that we want all the variables + # from the domain, that are predetermined ahead of time, and should save performance + j = 0 + for i in varkeys: + # var_path = f'{domain}/{i}' + var_path = domain + "/" + i + # print(var_path) + ix = state_paths[var_path] + state_ix[ix] = state_vals[j] + j += 1 + return True + + def hydr_state_vars(): return [ "DEP", @@ -221,7 +257,7 @@ def hydr_init_ix(state, domain): def sedtrn_state_vars(): - sedtrn_state = ["RSED4", "RSED5", "RSED6"] + sedtrn_state = ["RSED1", "RSED2", "RSED3", "RSED4", "RSED5", "RSED6"] return sedtrn_state