Equation

[rburghol]

Overview

Tasks

• Data Model in hdf5
• Parser (json, UCI)
• tokenize into solvable pieces
• Object Model: class or functional? (see Object Model below)
• Equation Solver: integer keyed tokenized #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)