rburghol
commented
1 year ago Scatch pad of dev code
def is_float_digit(n: str) -> bool:
try:
float(n)
return True
except ValueError:
return False
import numpy as np
import time
from numba.typed import Dict
from numpy import zeros
from numba import int8, float32, njit, types # import the types
# allops is created here as a Dict, but really it represents the
# hdf5 structure, which is allowed to mix string and numbers
# this was just created because I am testing and being sloppy
# but this particular structure will NOT make it into any @njit functions
# and thus can be a more flexible type
allops = Dict.empty(key_type=types.unicode_type, value_type=types.UnicodeCharSeq(128))
allops["/OBJECTS/RCHRES_0001/discharge_mgd/equation"] = "( (1.0 - consumption - unaccounted_losses) * wd_mgd + discharge_from_gw_mgd)"
# we parse the equation during readuci/pre-processing and break it into njit'able pieces
exprStack[:] = []
results = BNF().parseString(allops["/OBJECTS/RCHRES_0001/discharge_mgd/equation"], parseAll=True)
ps = []
ep = exprStack
pre_evaluate_stack(ep[:], ps)
# need to translate the constants in the equation to hdf5 variables.
# Ex:
# ps[0] = ['-', '1.0', 'consumption']
# 1.0 is a constant. The first constant in the discarge_mgd equation
# is 1.0, so we add it to the hdf5 as '_c1'
# since numeric values cannot be in the same Dict as strnigs in numba/njit
# these go directly into state as well as the hdf5
# then swap the path out for the value 1.0
# this can be done during parse UCI time since we have all python
#state = Dict.empty(key_type=types.unicode_type, value_type=types.float64)
state = Dict.empty(key_type=types.UnicodeCharSeq(128), value_type=types.float64)
parent_path = "/OBJECTS/RCHRES_0001"
object_path = "/OBJECTS/RCHRES_0001/discharge_mgd"
parent_state_path = "/STATE/RCHRES_0001"
state_path = "/STATE/RCHRES_0001/discharge_mgd"
num_co = 0
ops_path = object_path + "/_ops"
ops_state_path = state_path + "/_ops"
for i in range(0, len(ps) - 1):
op_name = "_op" + str(i)
op_path = ops_path + "/" + op_name
op_state_path = ops_state_path + "/" + op_name
for j in range(1,3):
arg_name = "arg" + str(j)
if ps[i][j] == None:
op_str = '_result' # we assume the result from previous step is to be used
ps[i][j] = '_result' # we assume the result from previous step is to be used
if is_float_digit(ps[i][j]):
print(ps[i][j], " is numeric")
num_co += 1
co_name = "_c" + str(num_co)
# where does the object info reside
co_path = op_path + "/" + co_name
# where does the numeric value for this reside
# this is the value that gets used in the actual njit evaluation functions
op_state_path = op_state_path + "/" + co_name
# stash the value here
state[op_state_path] = float(ps[i][j])
op_str = co_name
else:
# this is a variable, so we need to get its path
# spoofed for now
op_str = ps[i][j]
op_state_path = parent_state_path + "/" + ps[i][j]
print(op_path + "/" + arg_name + "/path")
allops[op_path + "/" + arg_name + "/path"] = op_state_path
allops[op_path + "/" + arg_name] = op_str
# now set
allops[op_path + "/op"] = ps[i][0]
#allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops"]
#allops["/OBJECTS/RCHRES_0001/discharge_mgd/_c1"]
# string functions available to us
#allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops"] = ps
#allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops"]
# use dict
import numpy as np
#tpvals = Dict.empty(key_type=types.unicode_type, value_type=types.float64)
tpvals = Dict.empty(key_type=types.UnicodeCharSeq(128), value_type=types.float64)
tpops = Dict.empty(key_type=types.unicode_type, value_type=types.UnicodeCharSeq(128))
tpvals['/STATE/RCHRES_0001/discharge_mgd/_c1'] = 1.0 # constants would be populated during equation parsing, if there were 3 constants in the equation they would be _c1, _c2, and _c3
tpvals['/STATE/RCHRES_0001/consumption'] = 0.15
tpvals['/STATE/RCHRES_0001/consumption'] = 0.15
tpops['op'] = '-'
tpops['arg1'] = '/STATE/RCHRES_0001/Qin/_c1'
tpops['arg2'] = '/STATE/RCHRES_0001/consumption'
exec_eqn_pnj(tpops, tpvals)
atpops = Dict.empty(key_type=types.unicode_type, value_type=types.UnicodeCharSeq(128))
atpops['op'] = allops['/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op0/op']
atpops['arg1'] = allops['/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op0/arg1/path']
atpops['arg2'] = allops['/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op0/arg2/path']
state['/STATE/RCHRES_0001/consumption'] = 0.15
exec_eqn_pnj(atpops, state)
import time
@njit
def iterate_pnj(atpops, state, steps):
checksum = 0.0
for step in range(steps):
checksum += exec_eqn_pnjopt(atpops, state)
#checksum += exec_eqn_pnj(atpops, state)
#exec_eqn_pnj(atpops, state)
return checksum
steps = 24 * 365 * 40
#steps = 24 * 365 * 1
start = time.time()
num = iterate_pnj(atpops, state, steps)
end = time.time()
print(end - start, "seconds")
@njit
def iterate_num(atpop, atpops, state, steps):
checksum = 0.0
for step in range(steps):
checksum += exec_eqn_num(atpop, atpops, state)
return checksum
# manually set an opnum = 1 (-) for testing
allops['/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op0/opnum'] = "1"
atpop = int(allops['/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op0/opnum'])
start = time.time()
num = iterate_num(atpop, atpops, state, steps)
end = time.time()
print(end - start, "seconds")
@njit
def exec_eqn_nall(op_tokens, state_ix):
# these 2 retrievals amount to approx. 35% of exec time.
# allops["/OBJECTS/RCHRES_0001/discharge_mgd/_s_index"] = 1
# allops["/OBJECTS/RCHRES_0001/discharge_mgd/_num_ops"] = 1
# allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op1/opi"] = 1 # 1: subtract, 2: add, 3: mult, 4: div
# allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op1/arg1/s_index"] = 2 # pointer to state array
# allops["/OBJECTS/RCHRES_0001/discharge_mgd/_ops/_op1/arg2/s_index"] = 3 # pointer to state array
# state_ix[2] = 0.15 # CU fraction state
# op_tokens # Executable configs index these by the state_ix variable key which is the global ID, _s_index
# op_tokens[1] = [1, 1, 2, 3] # [ op class (1 = eqn, 2=list), op type (1-4 for math), arg1 state id, arg2 state id]
op_class = op_tokens[0] # we actually will use this in the calling function, which will decide what
# next level function to use
# todo: see if skipping this temp variable op, val1, val2 setting saves time
#op = op_tokens[1]
#val1 = state_ix[op_tokens[2]]
#val2 = state_ix[op_tokens[3]]
result = 0
#return
#print("val1 = ", val1, "val2 = ", val2, "op = ", ops[0])
#if op == 1:
if op_tokens[1] == 1:
#result = val1 - val2
result = state_ix[op_tokens[2]] - state_ix[op_tokens[3]]
#elif op == 2:
elif op_tokens[1] == 2:
#result = val1 + val2
result = state_ix[op_tokens[2]] + state_ix[op_tokens[3]]
#elif op == 3:
elif op_tokens[1] == 3:
#result = val1 * val2
result = state_ix[op_tokens[2]] * state_ix[op_tokens[3]]
return result
#
op_tokens = Dict.empty(key_type=types.int64, value_type=types.i8[:])
op_tokens[1] = np.asarray([1, 1, 2, 3], dtype="i8")
op_tokens[2] = np.asarray([1, 2, 2, 3], dtype="i8")
state_ix = Dict.empty(key_type=types.int64, value_type=types.float64)
state_ix[2] = 1.0 # value of the variable constant from the equation for discharge_mgd
state_ix[3] = 0.15 # the value of the variable "consumption"
exec_eqn_nall(op_tokens[1], state_ix)
@njit
def iterate_nall(op_tokens, state_ix, steps):
checksum = 0.0
for step in range(steps):
state_ix[1] = exec_eqn_nall(op_tokens[1], state_ix)
checksum += state_ix[1]
#val = exec_eqn_nall(op_tokens[1], state_ix)
#state_ix[1] = val
#checksum += val
return checksum
start = time.time()
num = iterate_nall(op_tokens, state_ix, steps)
end = time.time()
print(end - start, "seconds")
for i in op_tokens.keys():
print(i)
# now add simple lookup support
@njit
def specl_lookup(data_table, keyval, lutype, valcol):
if lutype == 2: #stair-step
idx = (data_table[:, 0][0:][(data_table[:, 0][0:]- keyval) <= 0]).argmax()
luval = data_table[:, valcol][0:][idx]
elif lutype == 1: # interpolate
luval = np.interp(keyval,data_table[:, 0][0:], data_table[:, valcol][0:])
# show value at tis point
return luval
@njit
def exec_tbl_eval(op, state_ix, dict_ix):
ix = op[1]
dix = op[2]
mx_type = op[3] # not used yet, what type of table? in past this was always 1-d or 2-d
key1_ix = op[4]
#print("ix, dict_ix, mx_type, key1_ix", ix, dix, mx_type, key1_ix)
lutype = op[5]
valcol = op[8]
data_table = dict_ix[dix]
keyval = state_ix[key1_ix]
#print("Key, ltype, val", keyval, lutype, valcol)
result = specl_lookup(data_table, keyval, lutype, valcol)
return result
# load source from openmi-om/py/XdataMatrix.class.py
# this is state ix var 3
# doing a lookup table (like FTABLE) where we are given an arbitrary key column in a table and an arb value
# column and we want the output to be the interpolated in the value column according to the key column
# op is form: op_type, state_ix, dict_ix, mx_type (not used yet), key1_ix, key1_lu_type, key2_ix, key2_lu_type, val_col
op_tokens[3] = np.asarray([2, 3, 4, 0, 5, 2, 0, 0, 2], dtype="i8") # need to think long and hard about these required tokens
# notes:
# - dict_ix key is NOT always the same as the state_ix key
# since we may allow a "matrix accessor" to perform a lookup on another table
# - val_col is only used if this is a matrix_accessor and the lookup is a single column
state_ix[3] = 0.0 # the state var for this matrix accessor (since it is an ftable style lookup)
state_ix[4] = 0.0
# create a psuedo state variable for the storage in the ftable
state_ix[5] = 200000.0
dict_ix = Dict.empty(key_type=types.int64, value_type=types.float32[:,:])
dict_ix[4] = XdataMatrix.parseMatrix("[ [ 0.0, 170.0, 0], [195200.0, 240.0, 8890.6], [204252.8, 241.0, 9301.5], [213736.0, 242.0, 9712.5] ]")
specl_lookup(dict_ix[4], 200000, 2, 2)
exec_tbl_eval(op_tokens[3], state_ix, dict_ix)
@njit
def iterate_all_nall(op_tokens, state_ix, dict_ix, steps):
checksum = 0.0
for step in range(steps):
for i in op_tokens.keys():
if op_tokens[i][0] == 1:
state_ix[i] = exec_eqn_nall(op_tokens[i], state_ix)
elif op_tokens[i][0] == 2:
state_ix[i] = exec_tbl_eval(op_tokens[i], state_ix, dict_ix)
checksum += state_ix[i]
return checksum
start = time.time()
num = iterate_all_nall(op_tokens, state_ix, dict_ix, steps)
end = time.time()
print(end - start, "seconds")
Example equation evaluation test.
Code: String Keyed Operators and Operands
Load code from https://github.com/HARPgroup/openmi-om/blob/master/py/fourFn.py