Equation

[rburghol]

Overview

Tasks

• [x] Data Model in hdf5
• [x] Parser (json, UCI)
• [x] tokenize into solvable pieces
• [x] Object Model: class or functional? (see Object Model below)
• [x] #38
  • Decomposition/Run-time Optimization, can we split equation into individual 2 operand
  • ~eval():~ We cannot use eval() in numba - but maybe numexpr.evaluate() See below Solver: eval options
  • Pre-rendered functions (see below Solver: Pre-rendered Functions

Object Model

Explore method of using @jitclass objects or a more functional approach.

• https://numba.pydata.org/numba-doc/latest/user/jitclass.html#jitclass

Handler Class

class Equation(modelObject):
    # the following are supplied by the parent class: name, log_path, attribute_path, state_path, inputs
    op_trios = [] # string array with sorted operand1, operand2, operator

test = Equation()
test.ncode_preStep() # generate the numba code for preStep() actions.

Execution Rendering

Parser

• json coming in (already partially atomized and good for solving)

Data Model in hdf5

• Referencing related variables.
  • Ex:
  • Entity: INTAKE_001, dhdf5 reference: /TIMESERIES/INTAKE_001/
  • equation: available_mgd = Qintake - flowby
  • arguments:
    • Qreach: entity upstream flow + runoff in containing reach, /TIMESERIES/RCHRES_001/IVOL - could also be an equation scaling to intake drainage area, i.e. Qreach = ivol * local_area / area
    • flowby: local equation reference /TIMESERIES/INTAKE_001/flowby/value (note: flowby would also be an equation)
  • resolved out:
  • ts['/TIMESERIES/INTAKE_001/available_mgd/value'] = ts['/TIMESERIES/RCHRES_001/IVOL/value'] - ts['/TIMESERIES/INTAKE_001/flowby/value']
• With the above notation can we effectively create a matrix of mathematical during run prep, then execute those relationships very rapidly at runtime?
• Proposed Process:
  1. substitution: replace equation values with hdf5 paths
  2. retrieve: replace hdf5 paths with their values in the equation
  3. evaluate: do the math expression
  4. store: put data in relevant hdf5 location

Solver

Decomposition/Optimization:

• Break into RPN type pieces value and operand
• Equation parser: https://levelup.gitconnected.com/how-to-write-a-formula-string-parser-in-python-5362210afeab

Eval

• eval() there are examples of using numexpr.evaluate() (a numba analog of eval()) with @njit however,
  • it cannot accept array refrences in the equation, so, things like ts['/TIMESERIES/RCHRES_001/IVOL/value'] - ts['/TIMESERIES/INTAKE_001/flowby/value'] cannot be evaluated.
  • We would have to do a string substitution to replace the path with an actual value (like we did in om) BUT, the numba string replace methods do not allow replacing with a numeric value, only a string, and does not have the ability to convert a number to a string inside an @njit compiled function.

Pre-rendered Functions

• At compile time, write out a function with a list of mathematical statements, like: state['/STATE/INTAKE_001/available_mgd/value'] = state['/STATE/RCHRES_001/IVOL/value'] - state['/STATE/INTAKE_001/flowby/value']
• Could even have custom functions for each entity (but this would be a lot of overhead methinks, like:
  • def fn_specl_rchres_0001_Qin(state, ui, ts, step)

[rburghol]

Add a group with hdf5 module in python

• Note: data sets are accessible via the items() method as a type of collections object class ItemsView https://docs.python.org/3/library/collections.abc.html?highlight=itemsview#collections.abc.ItemsView

  
  import h5py
  import numpy as np
  import pandas as pd
  from collections import defaultdict
  from pandas import DataFrame, read_hdf, HDFStore
  from numba.typed import Dict

fpath = 'C:/Workspace/modeling/cbp/hsp2/river/OR1_7700_7980.h5'

close the already opened file first

f.close() f = h5py.File(fpath,'a') # use mode 'a' which allows read, write, modify objs = f.create_dataset("OBJECTS", (50,), dtype='f')

made a mistake, need to create a group, not a dataset

del f["OBJECTS"] grp = f.create_group("OBJECTS") grp.create_group("RCHRES_0001") rch = grp["RCHRES_0001"] rch.create_group("Qin") Qin = rch["Qin"]

attributes can simply be added here:

Qin.attrs['equation'] = 'Qup * 1.23' Qin.attrs['default_value'] = 0.0 Qin.attrs['state_path'] = "/STATE/RCHRES_0001/Qin" Qin.attrs['attribute_path'] = "/OBJECTS/RCHRES_0001/Qin" Qin.attrs['fn_step'] = "om_equation_step" # this can be used to specify any function whatsoever that takes the correct args, but might be difficult to use in njit? Qin.attrs['data_type'] = "f"

read attributes like this:

ts = f['TIMESERIES'].items() for a in ts: print(a[0], a[1])

dset = f['TIMESERIES/TS011/table'] dset[1]

(441766800000000000, 8.82848454)

dset[2]

(441770400000000000, 8.81089973)

dset[3]

(441774000000000000, 8.79396248)

dset[4]

(441777600000000000, 8.77740383)

[rburghol]

Now add a timeseries indexed by integers - ts_ix

• You must load function in the file https://github.com/HARPgroup/openmi-om/blob/master/py/demo.token.eqn.py

  
  ts_ix = Dict.empty(key_type=types.int64, value_type=types.float32[:])
  tsq = f['TIMESERIES/']
  ts_ts = np.asarray(tsq['TS011']['table']['index']) # index of timestamps for each timestep. use known TS for this demo
  # this code replicates a piece of the function of get_timeseries, and loads the full timeseries into a Dict
  # the dict will be keyed by a simple integer, and have values at the time step only.  The actual time 
  # at that step will be contained in ts_ts 
  ts_tables = list(tsq.keys())
  for i in ts_tables:
  if i == 'SUMMARY': continue # skip this non-numeric table
  var_path = '/TIMESERIES/' + i
  ix = set_state(state_ix, state_paths, var_path, 0.0)
  ts_ix[ix] = np.asarray(tsq[i]['table']['values'], dtype="float32")

- now, the Qin variable from the WDM can be found in the ts_ix array keyed as integer
  - we can get that integer by: `Qin_ix = get_state_ix(state_ix, state_paths, '/TIMESERIES/TS011')`
Iterate through the values for a simulation

Qin_ix = get_state_ix(state_ix, state_paths, '/TIMESERIES/TS011') for idx, x in np.ndenumerate(ts_ts): state_ix[Qin_ix] = ts_ix[Qin_ix][idx] # the timeseries inputs state get set here. must a special type of object to do this

view dict key:value pairs in console by rows

[print(key,':',value) for key, value in state_ix.items()]

[rburghol]

op_tokens, state_paths, state_ix, dict_ix = init_sim_dicts()

eqn_path = "/OBJECTS/RCHRES_0001/Qin"
eqn = "Qin * 1.21"
eq_ix = set_state(state_ix, state_paths, eqn_path, 0.0) # initialize this equation operator in the state_ix Dict, returns unique key

# NR uses weird global shit so we have to do this:
exprStack = []
exprStack[:] = []
ps = deconstruct_equation(eqn)