Passing referenced OUTDGT (and others) into _hydr_()

[rburghol]

Overview

Special actions can pull and alter data from all locations in the simulation, i.e., from multiple RCHRES or PERLND entities. In order to support special actions with minimal data retrieval and read/write overhead, it is desirable to pass in a python view of all relevant data, but there are some challenges in the existing code base:

• Non-Referenced Data passed to Sub-routines: In the hydr() function the OUTDGT data should be modifiable by specl() when called inside _hydr_(). However, in the current hsp2 implementation, OUTDGT is passed from hydr() to _hydr_() as a clone, i.e., it is not a pass by reference view.
• Domain Specific Data Only to Sub-Routines: In main.py the routine assembles ui and ts as slices of the uci and based on the operation and segment.
  • Code from main.py: ts = get_timeseries(io_manager,ddext_sources[(operation,segment)],siminfo)
• What variables would we alter? See test withdrawal section in #12 Given the location of execution, when would this action take place (i.e., are we making a change for the current time-step or next time-step?)

Tasks:

• [x] check all 3 and see if we can pass in to _hydr_()
• [x] Propose final, generic specl() function call and data-prep approach (inside the calling function) for passing referenced objects into specl() without needing extra attributes to link back to local state variables in the calling function, such as OUTDGT.
• [x] Successfully alter OUTDGT and see if this changes flows?: https://github.com/HARPgroup/HSPsquared/wiki/JK-HSP2-Log#week-of-10322
  • [x] Test 1: Setting the following inside SPECL.py: ts['OUTDGT2'][step] = 99
    • Result: there was no observed change in output h5 HYDR_table
  • [x] Test 2: Setting OUTDGT2 within SPECL.py as above THEN setting the following in HYDR.py: OUTDGT[step, :] = [ts['OUTDGT1'][step], ts['OUTDGT2'][step], 0.0]
    • Result: Now in output h5 HYDR_table, O2 = 99 for many timesteps
• [x] Does changing OUTDGT in state (i.e. ts[step -1]) change the variable outdgt in the main _hydr_() calling function?
  • [x] Ans. No: We found that modifying OUTDGT2 in ts does not automatically alter outdgt: see notes on modifying OUTDGT2: https://github.com/HARPgroup/HSPsquared/wiki/JK-HSP2-Log#work-session-92322

Development

The following code snippet enables this referenced data, however we have not yet gotten it to function inside of _hydr_(), however, that may be because we hard-coded the array of OUTDGT values to only grab 2 of the 3 that are in there.

 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])

[rburghol]

@jdkleiner I want to get some proposed ways forward on this one, specifically to address the goal in the 2nd bullet.

Thoughts:

• Create and pass in a Dict called state (in addition to ui and ts) for exchanging values such as OUTDGT
• Functions that call specl() are responsible for setting things into state such as OUTDGT which have been explicitly separated from the ts Dict by the calling function (as is done by _hydr_ )
• could also have a helper function before calling specl() to automatically populate the state variable based on a list of local state variables of interest (also used by specl to limit the computations to only those variables that are used by the calling function)

[jdkleiner]

@rburghol

I think the idea of using a separate state dictionary is sensible. I had some trouble running hsp2 last week, but all seems to be running smoothly this morning! I plan on focusing on this today, so should be able to test the bullets you've proposed here and will report back on progress.

[jdkleiner]

@rburghol

• Findings & reminders
  1. ts dictionary is an input to hydr()
     • The ts dictionary has many keys including OUTDGT1 and OUTDGT2
  2. OUTDGT is a NumPy array created from ts within hydr()
     • OUTDGT is a two dimensional array of type float64 in FORTRAN array order array(float64, 2d, F)
     • Each time step has its own "row" 1-dimensional array

• Tested building a state dictionary within _hydr_() (working on branch specdev)

  • Am able to create the empty dictionary, and the add test key:value pairs of floating numbers
  • Unsuccessful in populating OUTDGT on the state dictionary
    • *I suspect this is why OUTDGT is being passed into _hydr_() as a separate input param, because there's no easy way to add a numpy array to a dictionary

  • # import numpy as np
    # from numba import types
    # from numba.typed import Dict
    
    # Dict.empty() constructs a typed dictionary
    #   The key and value types must be explicitly declared
    #   This will create a dict with keys as string and values of type float
    state = Dict.empty(key_type=types.unicode_type, value_type=types.float64)
    print(state)
    {}
    
    state['hello'] = 1.5
    print(state)
    {hello: 1.5}
    
    # tried many variations of setting an outdgt variable on the state dict
    state['OUTDGT'] = OUTDGT
    results in errors

    • Note: went the route of building the dictionary in a custom fashion rather than using make_numba_dict() which is set up specifically for uci files

• Q: why not just pass in OUTDGT directly to specl()?

  • Performed some tests of this route
    • def specl(ui, ts, OUTDGT, step, specactions):
    • def _specl_(ui, ts, OUTDGT, step, specactions):

  • 
    # before specl() call:
    # OUTDGT:
    [[73.66304779 -0.          0.        ]
    [73.66304779 -0.          0.        ]
    [73.66304779 -0.          0.        ]
    ...
    [83.55184937 -0.          0.        ]
    [83.55184937 -0.          0.        ]
    [83.55184937 -0.          0.        ]]
    
    # within _specl_():
    OUTDGT[step, 0] = 99
    
    # after specl() call:
    # OUTDGT:
    [[99. -0.  0.]
    [99. -0.  0.]
    [99. -0.  0.]
    ...
    [99. -0.  0.]
    [99. -0.  0.]
    [99. -0.  0.]]

  • Setting OUTDGT[step, 0] = 99 inside _specl_() does result in O1 = 99 in the output HYDR table.
    • see RESULTS/RCHRES_R001/HYDR/table
    • This works because specl() call happens before setting of outdgt[:] = OUTDGT[step, :]
  • Modifying a value on OUTDGT with this method does not automatically get reflected in ts
    • Is that okay? Or do we need to ensure ts is also being modified in the process?

[rburghol]

Easy Qs first:

• TS in the hdf5 doesn't need to be changed since that is static inputs, provided that the state variable OUTDGT gets changed which results in the final /RESULTS/RCHRES... array having the altered O1, O2, ...

The difficult stuff: I think that the reason that OUTDGT caused an error is the Dict declaration:

state = Dict.empty(key_type=types.unicode_type, value_type=types.float64)

Does not accept arrays, so changing to this would remove the error:

state = Dict.empty(key_type=types.unicode_type, value_type=types.float64[:])

However:

• I don't really like that since it forces every state var to be an array and that seems inefficient but maybe should be considered?
• Alternatively, since OUTDGT is stored as O1,O2,... we could put them in state in that way.
• Ultimately, it shouldn't be the responsibility of the specl() function to handle a custom variable especially one that is domain specific like OUTDGT. Thus, conversions to and from array notation like that should be handled by the calling routine

All those considered I think we choose to:

1. Deconstruct and reconstruct OUTDGT into it's component pieces in _hydr_. Which makes the steps:
   • populate state
   • call specl()
   • copy vars from state back to local vars (like OUTDGT)
2. Or, make all state vars an array even tho the majority will be 1 element long.
   • Note: I think that there is not way to create a Dict that gets values from multiple Dicts and still has things propagate via pass-by-reference, we will still need to do the populate, call specl(), harvest steps just like option 1.
3. The most important item might be performance testing. Dicts that use a string name as a key, instead of an integer can show really slow performance. Which might be a reason that developers have passed these Dicts as separate arguments. But this is not practical with specl since ANY piece of info could change and that can quickly get absurd in terms of the number of arguments needed. See the performance testing for some interesting info here:
   • All operators and operands for a sample equation is really slow: Section "Early String Keyed Variation Was Slow" in body of issue https://github.com/HARPgroup/HSPsquared/issues/38
   • Pass the equation operator (i.e. +,-,*,/) in as integer indices speeds up 40%: https://github.com/HARPgroup/HSPsquared/issues/38#issuecomment-1337963868
   • Pass equation operands in as integer keyed items and set result data to integer keyed Dict is really fast: https://github.com/HARPgroup/HSPsquared/issues/38

Also, for computational efficiency we may want to pass a var to specl() that checks to see if anything has been changed and if nothing gets changed no need to copy.

@jdkleiner thoughts on the above?

[jdkleiner]

Test implementation of state dictionary approach for specl:

See HYDR.py:

# The following occurs at the beginning of the HYDR for-loop

# set up state dictionary
state = Dict.empty(key_type=types.int64, value_type=types.float64)
state[1] = OUTDGT[step, 0]
state[2] = OUTDGT[step, 1]
state[3] = OUTDGT[step, 2]

# specl calls _specl_ which modifies state variables
errors_specl = specl(ui, ts, state, step, specactions)

# set OUTDGT arry using the modified state variables 
OUTDGT[step, :] = [state[1], state[2], state[3]]

See SPECL.py:

@njit
def _specl_(ui, ts, state, step, specactions):
    errors_specl = zeros(int(1)).astype(int64)

    state[1] += 1000
    state[2] = 99

    return errors_specl

Result console output:

# hsp2 run PL3_5250_0001.h5

step 8669 of 8760
state before specl()
1 : 83.55184936523438
2 : -0.0
3 : 0.0
state after specl()
1 : 1083.5518493652344
2 : 99.0
3 : 0.0
OUTDGT[step, :] [1083.55184937   99.            0.        ]