diff --git a/.gitignore b/.gitignore index f177d0cd..6c2ae13f 100644 --- a/.gitignore +++ b/.gitignore @@ -56,16 +56,19 @@ tests/_LargeFileStash # Temporary files _TutorialData/ + # Specific files run.py tests/GLWACSO/HSP2results/hspp007.hdf tests/GLWACSO/HSPFresults/hspf006.HBN tests/GLWACSO/HSP2results/hspp007.uci tests/test_report_conversion.html + +# Omit big files tests/land_spec/hwmA51800.h5 tests/testcbp/HSP2results/PL3_5250_0001.h5 -tests/testcbp/HSP2results/PL3_5250_specl.h5 tests/testcbp/HSP2results/*.csv +tests/test10/HSP2results/test10.h5 # R files .Rdata diff --git a/HSP2/HYDR.py b/HSP2/HYDR.py index fb8ed373..f9211e76 100644 --- a/HSP2/HYDR.py +++ b/HSP2/HYDR.py @@ -155,7 +155,6 @@ def hydr(io_manager, siminfo, uci, ts, ftables, state): # initialize the hydr paths in case they don't already reside here hydr_init_ix(state_ix, state_paths, state['domain']) op_tokens = state['op_tokens'] - print("state_ix is type", type(state_ix)) ####################################################################################### # Do the simulation with _hydr_ (ie run reaches simulation code) @@ -173,7 +172,7 @@ def hydr(io_manager, siminfo, uci, ts, ftables, state): return errors, ERRMSGS -@njit(cache=True) +@njit def _hydr_(ui, ts, COLIND, OUTDGT, rowsFT, funct, Olabels, OVOLlabels, state_info, state_paths, state_ix, dict_ix, ts_ix, state_step_hydr, op_tokens, model_exec_list): errors = zeros(int(ui['errlen'])).astype(int64) diff --git a/HSP2/om.py b/HSP2/om.py index c2d72f2c..3a7c5d26 100644 --- a/HSP2/om.py +++ b/HSP2/om.py @@ -13,7 +13,7 @@ import time from numba.typed import Dict from numpy import zeros, int32 -from numba import int8, float32, njit, types, typed # import the types +from numba import int8, float64, njit, types, typed # import the types import random # this is only used for a demo so may be deprecated from HSP2.state import * @@ -55,11 +55,11 @@ def is_float_digit(n: str) -> bool: # Import Code Classes from HSP2.om_model_object import * from HSP2.om_sim_timer import * -#from HSP2.om_equation import * +from HSP2.om_equation import * from HSP2.om_model_linkage import * from HSP2.om_special_action import * #from HSP2.om_data_matrix import * -#from HSP2.om_model_broadcast import * +from HSP2.om_model_broadcast import * #from HSP2.om_simple_channel import * #from HSP2.om_impoundment import * from HSP2.utilities import versions, get_timeseries, expand_timeseries_names, save_timeseries, get_gener_timeseries @@ -78,7 +78,6 @@ def init_om_dicts(): def state_load_om_json(state, io_manager, siminfo): # - model objects defined in file named '[model h5 base].json -- this will populate an array of object definitions that will # be loadable by "model_loader_recursive()" - model_data = state['model_data'] # JSON file would be in same path as hdf5 hdf5_path = io_manager._input.file_path (fbase, fext) = os.path.splitext(hdf5_path) @@ -90,8 +89,10 @@ def state_load_om_json(state, io_manager, siminfo): jfile = open(fjson) json_data = json.load(jfile) # dict.update() combines the arg dict with the base - model_data.update(json_data) - state['model_data'] = model_data + state['model_data'].update(json_data) + # merge in the json siminfo data + if 'siminfo' in state['model_data'].keys(): + siminfo.update(state['model_data']['siminfo']) return def state_load_om_python(state, io_manager, siminfo): @@ -133,30 +134,42 @@ def state_load_dynamics_om(state, io_manager, siminfo): state_load_om_json(state, io_manager, siminfo) return -def state_om_model_run_prep(state, io_manager, siminfo): +def state_om_model_root_object(state, siminfo): # Create the base that everything is added to. this object does nothing except host the rest. - model_root_object = ModelObject("") - # set up the timer as the first element - timer = SimTimer('timer', model_root_object, siminfo) - + if 'model_root_object' not in state.keys(): + model_root_object = ModelObject("") + state['model_root_object'] = model_root_object + # set up the timer as the first element + if '/STATE/timer' not in state['state_paths'].keys(): + timer = SimTimer('timer', model_root_object, siminfo) + + +def state_om_model_run_prep(state, io_manager, siminfo): + # insure model base is set + state_om_model_root_object(state, siminfo) # now instantiate and link objects # state['model_data'] has alread been prepopulated from json, .py files, hdf5, etc. + model_root_object = state['model_root_object'] model_loader_recursive(state['model_data'], model_root_object) print("Loaded objects & paths: insures all paths are valid, connects models as inputs") # both state['model_object_cache'] and the model_object_cache property of the ModelObject class def # will hold a global repo for this data this may be redundant? They DO point to the same datset? # since this is a function that accepts state as an argument and these were both set in state_load_dynamics_om # we can assume they are there and functioning - model_object_cache = state['model_object_cache'] + if 'model_object_cache' in state.keys(): + model_object_cache = state['model_object_cache'] + else: + model_object_cache = ModelObject.model_object_cache model_path_loader(model_object_cache) # len() will be 1 if we only have a simtimer, but > 1 if we have a river being added model_exec_list = [] # put all objects in token form for fast runtime execution and sort according to dependency order print("Tokenizing models") + if 'ops_data_type' in siminfo.keys(): + ModelObject.ops_data_type = siminfo['ops_data_type'] # allow override of dat astructure settings ModelObject.op_tokens = ModelObject.make_op_tokens(max(ModelObject.state_ix.keys()) + 1) model_tokenizer_recursive(model_root_object, model_object_cache, model_exec_list) op_tokens = ModelObject.op_tokens - print("op_tokens has", len(op_tokens),"elements") # model_exec_list is the ordered list of component operations #print("model_exec_list(", len(model_exec_list),"items):", model_exec_list) # This is used to stash the model_exec_list in the dict_ix, this might be slow, need to verify. @@ -165,10 +178,18 @@ def state_om_model_run_prep(state, io_manager, siminfo): state['model_object_cache'] = model_object_cache state['model_exec_list'] = np.asarray(model_exec_list, dtype="i8") if ModelObject.ops_data_type == 'ndarray': - state['state_ix'] = np.asarray(list(state['state_ix'].values()), dtype="float64") + state_keyvals = np.asarray(zeros(max(ModelObject.state_ix.keys()) + 1), dtype="float64") + for ix, val in ModelObject.state_ix.items(): + state_keyvals[ix] = val + state['state_ix'] = state_keyvals + else: + state['state_ix'] = ModelObject.state_ix state['op_tokens'] = op_tokens if len(op_tokens) > 0: state['state_step_om'] = 'enabled' + print("op_tokens is type", type(op_tokens)) + print("state_ix is type", type(state['state_ix'])) + print("op_tokens has", len(op_tokens),"elements, with ", len(model_exec_list),"executable elements") return # model class reader @@ -271,6 +292,11 @@ def model_class_translate(model_props, object_class): if object_class == 'hydroImpSmall': model_props['object_class'] = 'Impoundment' print("Handling hydroImpSmall as Impoundment") + # now handle disabled classes - this is temporary to prevent having to comment and uncomment + disabled_classes = {'SimpleChannel', 'Impoundment', 'DataMatrix', 'dataMatrix'} + if model_props['object_class'] in disabled_classes: + print("Disabling class", model_props['object_class'], 'rendering as ModelObject') + model_props['object_class'] = 'ModelObject' def model_loader_recursive(model_data, container): k_list = model_data.keys() @@ -447,7 +473,7 @@ def pre_step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step): for i in model_exec_list: if op_tokens[i][0] == 12: # register type data (like broadcast accumulators) - pass#pre_step_register(op_tokens[i], state_ix, dict_ix) + pre_step_register(op_tokens[i], state_ix) return @njit @@ -467,7 +493,7 @@ def step_one(op_tokens, ops, state_ix, dict_ix, ts_ix, step, debug = 0): if debug == 1: print("DEBUG: Operator ID", ops[1], "is op type", ops[0]) if ops[0] == 1: - pass #step_equation(ops, state_ix) + step_equation(ops, state_ix) elif ops[0] == 2: # todo: this should be moved into a single function, # with the conforming name step_matrix(op_tokens, ops, state_ix, dict_ix) @@ -500,40 +526,14 @@ def step_one(op_tokens, ops, state_ix, dict_ix, ts_ix, step, debug = 0): @njit def step_model_test(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step, debug_step = -1): - val = 0 for i in model_exec_list: ops = op_tokens[i] + val = 0 if (step == debug_step): print("Exec'ing step ", step, " model ID", i) # op_tokens is passed in for ops like matrices that have lookups from other # locations. All others rely only on ops - # todo: decide if all step_[class() functions should set value in state_ix instead of returning value? - val = 0 - if ops[0] == 1: - step_equation(ops, state_ix) - elif ops[0] == 2: - # todo: this should be moved into a single function, - # with the conforming name step_matrix(op_tokens, ops, state_ix, dict_ix) - if (ops[1] == ops[2]): - # this insures a matrix with variables in it is up to date - # only need to do this if the matrix data and matrix config are on same object - # otherwise, the matrix data is an input and has already been evaluated - state_ix[ops[1]] = exec_tbl_values(ops, state_ix, dict_ix) - if (ops[3] > 0): - # this evaluates a single value from a matrix if the matrix is configured to do so. - state_ix[ops[1]] = exec_tbl_eval(op_tokens, ops, state_ix, dict_ix) - elif ops[0] == 3: - step_model_link(ops, state_ix, ts_ix, step) - continue - elif ops[0] == 5: - step_sim_timer(ops, state_ix, dict_ix, ts_ix, step) - elif ops[0] == 9: - continue - elif ops[0] == 13: - step_simple_channel(ops, state_ix, dict_ix, step) - # Op 100 is Basic ACTION in Special Actions - elif ops[0] == 100: - step_special_action(ops, state_ix, dict_ix, step) + step_one(op_tokens, op_tokens[i], state_ix, dict_ix, ts_ix, step, 0) return @njit @@ -557,16 +557,6 @@ def step_model_pcode(model_exec_list, op_tokens, state_info, state_paths, state_ def post_step_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step): return -@njit -def test_model(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step): - val = 0 - for i in model_exec_list: - print(i) - print(op_tokens[i][0]) - print(op_tokens[i]) - step_one(op_tokens, op_tokens[i], state_ix, dict_ix, ts_ix, step, 0) - return - def step_object(thisobject, step): # this calls the step for a given model object and timestep # this is a workaround since the object method ModelObject.step() fails to find the step_one() function ? @@ -581,28 +571,19 @@ def pre_step_test(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step): # op = op_tokens[i] if ops[0] == 12: # register type data (like broadcast accumulators) - pre_step_register(ops, state_ix, dict_ix) - continue + pre_step_register(ops, state_ix) + #continue #elif ops[0] == 1: # # register type data (like broadcast accumulators) # continue return -@njit -def test_perf(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step): -# for i in model_exec_list: -# ops = op_tokens[i] - for ops in op_tokens: - #step_one(op_tokens, ops, state_ix, dict_ix, ts_ix, step) - continue - return - @njit def iterate_perf(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, steps, debug_step = -1): checksum = 0.0 for step in range(steps): pre_step_test(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step) - step_model_test(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step, debug_step) + #step_model_test(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, step, debug_step) #print("Steps completed", step) return checksum diff --git a/HSP2/om_equation.py b/HSP2/om_equation.py new file mode 100644 index 00000000..87c2273f --- /dev/null +++ b/HSP2/om_equation.py @@ -0,0 +1,447 @@ +""" +The class Equation is used to translate an equation in text string form into a tokenized model op code +The equation will look for variable names inside the equation string (i.e. not numeric, not math operator) +and will then search the local object inputs and the containing object inputs (if object has parent) for +the variable name in question. Ultimately, everyting becomes either an operator or a reference to a variable +in the state_ix Dict for runtime execution. +""" +from HSP2.om import * +from HSP2.state import * +from HSP2.om_model_object import * +from numba import njit +class Equation(ModelObject): + # the following are supplied by the parent class: name, log_path, attribute_path, state_path, inputs + + def __init__(self, name, container = False, model_props = {}): + super(Equation, self).__init__(name, container, model_props) + self.equation = self.handle_prop(model_props, 'equation') + self.ps = False + self.ps_names = [] # Intermediate with constants turned into variable references in state_paths + self.var_ops = [] # keep these separate since the equation functions should not have to handle overhead + self.optype = 1 # 0 - shell object, 1 - equation, 2 - datamatrix, 3 - input, 4 - broadcastChannel, 5 - ? + self.non_neg = self.handle_prop(model_props, 'nonnegative', False, 0) + min_value = self.handle_prop(model_props, 'minvalue', False, 0.0) + self.min_value_ix = ModelConstant('min_value', self, min_value).ix + self.deconstruct_eqn() + + def handle_prop(self, model_props, prop_name, strict = False, default_value = None ): + prop_val = super().handle_prop(model_props, prop_name, strict, default_value ) + if (prop_name == 'equation'): + if type(prop_val) is str: + return prop_val + elif prop_val == None: + # try for equation stored as normal propcode + prop_val = str(self.handle_prop(model_props, 'value', True)) + return prop_val + + def deconstruct_eqn(self): + exprStack = [] + exprStack[:] = [] + self.ps = deconstruct_equation(self.equation) + # if ps is empty, we may have gotten a constant, so we will check it, + # and create a set of ps [constant] + 0.0 and return + # if this does not yield ps > 0 we will throw an error + if (len(self.ps) == 0): + tps = deconstruct_equation(self.equation + " + 0.0") + if len(tps) == 1: + # seemed to have succeeded, try to use this now + self.ps = tps + self.equation = self.equation + " + 0.0" + else: + raise Exception("Error: Cannot parse equation: " + self.equation + " on object " + self.name + " Halting.") + if (len(self.ps) == 1): + # this is a single set of operands, but we need to check for a solo negative number + # which will also get returned as one op set, with the first slot a number, then 0, 0 + # the first token *should* be an operator, followed by 2 operands + if is_float_digit(self.ps[0][0]): + # if it is longer than 1 character + if (self.ps[0][1] == 0) and (self.ps[0][2] == 0): + tps = deconstruct_equation(" 0.0 " + self.equation) + if len(tps) == 1: + # seemed to have succeeded, try to use this now + self.ps = tps + self.equation = " 0.0 " + self.equation + else: + raise Exception("Error: Cannot parse equation: " + self.equation + " on object " + self.name + " Halting.") + + #print(exprStack) + + def find_paths(self): + super().find_paths() + self.paths_found = False # override parent setting until we verify everything + #return + # we now want to check to see that all variables can be found (i.e. path exists) + # and we want to create variables for any constants that we have here + # do not handle mathematical operators + self.ps_names = [] + for i in range(len(self.ps)): + name_op = ["", "", ""] + name_op[0] = self.ps[i][0] + for j in range(1,3): + # range 1,3 counts thru 1 and 2 (not 3, cause python is so array centric it knows you know) + op_value = self.ps[i][j] + if op_value == None: + # don't need to check these as they are just referring to the stack. + continue + if is_float_digit(op_value): + op_name = "_op_" + str(i) + "_" + str(j) + else: + op_name = op_value + #print("Checking op set", i, "op", j, ":", op_name) + # constant_or_path() looks at name and path, since we only have a var name, we must assume + # the path is either a sibling or child variable or explicitly added other input, so this should + # resolve correctly, but we have not tried it + var_ix = self.constant_or_path(op_name, op_value, False) + # we now return, trusting that the state_path for each operand + # is stored in self.inputs, with the varix saved in self.inputs_ix + name_op[j] = op_name + self.ps_names.append(name_op) + self.paths_found = True + return + + def tokenize_ops(self): + self.deconstruct_eqn() + self.var_ops = tokenize_ops(self.ps) + + def tokenize_vars(self): + # now stash the string vars as new state vars + for j in range(2,len(self.var_ops)): + if isinstance(self.var_ops[j], int): + continue # already has been tokenized, so skip ahead + elif is_float_digit(self.var_ops[j]): + # must add this to the state array as a constant + constant_path = self.state_path + '/_ops/_op' + str(j) + s_ix = set_state(self.state_ix, self.state_paths, constant_path, float(self.var_ops[j]) ) + self.var_ops[j] = s_ix + else: + # this is a variable, must find it's data path index + var_path = self.find_var_path(self.var_ops[j]) + s_ix = get_state_ix(self.state_ix, self.state_paths, var_path) + if s_ix == False: + print("Error: unknown variable ", self.var_ops[j], "path", var_path, "index", s_ix) + print("searched: ", self.state_paths, self.state_ix) + return + else: + self.var_ops[j] = s_ix + + def tokenize(self): + # call parent to render standard tokens + super().tokenize() + # replaces operators with integer code, + # and turns the array of 3 value opcode arrays into a single sequential array + self.tokenize_ops() + # finds the ix value for each operand + self.tokenize_vars() + # renders tokens for high speed execution + self.ops = self.ops + [self.non_neg, self.min_value_ix] + self.var_ops + + +from pyparsing import ( + Literal, + Word, + Group, + Forward, + alphas, + alphanums, + Regex, + ParseException, + CaselessKeyword, + Suppress, + delimitedList, +) +import math +import operator + +exprStack = [] + + +def push_first(toks): + exprStack.append(toks[0]) + + +def push_unary_minus(toks): + for t in toks: + if t == "-": + exprStack.append("unary -") + else: + break + +def deconstruct_equation(eqn): + """ + We should get really good at using docstrings... + + we parse the equation during readuci/pre-processing and break it into njit'able pieces + this forms the basis of our object parser code to run at import_uci step + """ + results = BNF().parseString(eqn, parseAll=True) + ps = [] + ep = exprStack + pre_evaluate_stack(ep[:], ps) + return ps + +def get_operator_token(op): + # find the numerical token for an operator + # returns integer value, or 0 if this is not a recorgnized mathematical operator + if op == '-': opn = 1 + elif op == '+': opn = 2 + elif op == '*': opn = 3 + elif op == '/': opn = 4 + elif op == '^': opn = 5 + else: opn = False + return opn + +def tokenize_ops(ps): + '''Translates a set of string operands into integer keyed tokens for faster execution.''' + tops = [len(ps)] # first token is number of ops + for i in range(len(ps)): + op = get_operator_token(ps[i][0]) + # a negative op code indicates null + # this should cause no confusion since all op codes are references and none are actual values + if ps[i][1] == None: o1 = -1 + else: o1 = ps[i][1] + if ps[i][2] == None: o2 = -1 + else: o2 = ps[i][2] + tops.append(op) + tops.append(o1) + tops.append(o2) + return tops + +bnf = None + + +def BNF(): + """ + expop :: '^' + multop :: '*' | '/' + addop :: '+' | '-' + integer :: ['+' | '-'] '0'..'9'+ + atom :: PI | E | real | fn '(' expr ')' | '(' expr ')' + factor :: atom [ expop factor ]* + term :: factor [ multop factor ]* + expr :: term [ addop term ]* + """ + global bnf + if not bnf: + # use CaselessKeyword for e and pi, to avoid accidentally matching + # functions that start with 'e' or 'pi' (such as 'exp'); Keyword + # and CaselessKeyword only match whole words + e = CaselessKeyword("E") + pi = CaselessKeyword("PI") + # fnumber = Combine(Word("+-"+nums, nums) + + # Optional("." + Optional(Word(nums))) + + # Optional(e + Word("+-"+nums, nums))) + # or use provided pyparsing_common.number, but convert back to str: + # fnumber = ppc.number().addParseAction(lambda t: str(t[0])) + fnumber = Regex(r"[+-]?\d+(?:\.\d*)?(?:[eE][+-]?\d+)?") + ident = Word(alphas, alphanums + "_$") + + plus, minus, mult, div = map(Literal, "+-*/") + lpar, rpar = map(Suppress, "()") + addop = plus | minus + multop = mult | div + expop = Literal("^") + + expr = Forward() + expr_list = delimitedList(Group(expr)) + # add parse action that replaces the function identifier with a (name, number of args) tuple + def insert_fn_argcount_tuple(t): + fn = t.pop(0) + num_args = len(t[0]) + t.insert(0, (fn, num_args)) + + fn_call = (ident + lpar - Group(expr_list) + rpar).setParseAction( + insert_fn_argcount_tuple + ) + atom = ( + addop[...] + + ( + (fn_call | pi | e | fnumber | ident).setParseAction(push_first) + | Group(lpar + expr + rpar) + ) + ).setParseAction(push_unary_minus) + + # by defining exponentiation as "atom [ ^ factor ]..." instead of "atom [ ^ atom ]...", we get right-to-left + # exponents, instead of left-to-right that is, 2^3^2 = 2^(3^2), not (2^3)^2. + factor = Forward() + factor <<= atom + (expop + factor).setParseAction(push_first)[...] + term = factor + (multop + factor).setParseAction(push_first)[...] + expr <<= term + (addop + term).setParseAction(push_first)[...] + bnf = expr + return bnf + + +# map operator symbols to corresponding arithmetic operations +epsilon = 1e-12 +opn = { + "+": operator.add, + "-": operator.sub, + "*": operator.mul, + "/": operator.truediv, + "^": operator.pow, +} + +fn = { + "sin": math.sin, + "cos": math.cos, + "tan": math.tan, + "exp": math.exp, + "abs": abs, + "trunc": int, + "round": round, + "sgn": lambda a: -1 if a < -epsilon else 1 if a > epsilon else 0, + # functionsl with multiple arguments + "multiply": lambda a, b: a * b, + "hypot": math.hypot, + # functions with a variable number of arguments + "all": lambda *a: all(a), +} + +fns = { + "sin": "math.sin", + "cos": "math.cos", + "tan": "math.tan", + "exp": "math.exp", + "abs": "abs", + "trunc": "int", + "round": "round", +} + + +def evaluate_stack(s): + op, num_args = s.pop(), 0 + if isinstance(op, tuple): + op, num_args = op + if op == "unary -": + return -evaluate_stack(s) + if op in "+-*/^": + # note: operands are pushed onto the stack in reverse order + op2 = evaluate_stack(s) + op1 = evaluate_stack(s) + return opn[op](op1, op2) + elif op == "PI": + return math.pi # 3.1415926535 + elif op == "E": + return math.e # 2.718281828 + elif op in fn: + # note: args are pushed onto the stack in reverse order + args = reversed([evaluate_stack(s) for _ in range(num_args)]) + return fn[op](*args) + elif op[0].isalpha(): + raise Exception("invalid identifier '%s'" % op) + else: + # try to evaluate as int first, then as float if int fails + try: + return int(op) + except ValueError: + return float(op) + +def pre_evaluate_stack(s, ps): + op, num_args = s.pop(), 0 + if isinstance(op, tuple): + op, num_args = op + if op == "unary -": + ps.append([-evaluate_stack(s), 0, 0]) + return + if op in "+-*/^": + # note: operands are pushed onto the stack in reverse order + op2 = pre_evaluate_stack(s, ps) + op1 = pre_evaluate_stack(s, ps) + ps.append([ op, op1, op2]) + return + elif op == "PI": + ps.append([math.pi, 0, 0]) # 3.1415926535 + return + elif op == "E": + ps.append([math.e, 0, 0]) # 2.718281828 + return + elif op in fns: + # note: args are pushed onto the stack in reverse order + #print("s:", s, "op", op) + args = [] + for x in range(num_args): + args.append(pre_evaluate_stack(s, ps)) + args.reverse() + args.insert(fns[op], 0) + ps.append(args) + return + elif op[0].isalpha(): + return op + else: + # return the operand now + return op + + +@njit(cache=True) +def evaluate_eq_ops(op, val1, val2): + if op == 1: + #print(val1, " - ", val2) + result = val1 - val2 + return result + if op == 2: + #print(val1, " + ", val2) + result = val1 + val2 + return result + if op == 3: + #print(val1, " * ", val2) + result = val1 * val2 + return result + if op == 4: + #print(val1, " / ", val2) + result = val1 / val2 + return result + if op == 5: + #print(val1, " ^ ", val2) + result = pow(val1, val2) + return result + return 0 + + +@njit +def step_equation(op_token, state_ix): + op_class = op_token[0] # we actually use this in the calling function, which will decide what + # next level function to use + result = 0 + s = np.array([0.0]) + s_ix = -1 # pointer to the top of the stack + s_len = 1 + # handle special equation settings like "non-negative", etc. + non_neg = op_token[2] + min_ix = op_token[3] + num_ops = op_token[4] # this index is equal to the number of ops common to all classes + 1. See om_model_object for base ops and adjust + op_loc = 5 # where do the operators and operands start in op_token + #print(num_ops, " operations") + # is the below faster since it avoids a brief loop and a couple ifs for 2 op equations? + if num_ops == 1: + result = evaluate_eq_ops(op_token[op_loc], state_ix[op_token[op_loc + 1]], state_ix[op_token[op_loc + 2]]) + else: + for i in range(num_ops): + # the number of ops common to all classes + 1 (the counter for math operators) is offset for this + # currently 3 (2 common ops (0,1), plus 1 to indicate number of equation operand sets(2), so this is ix 3) + op = op_token[op_loc + 3*i] + t1 = op_token[op_loc + 3*i + 1] + t2 = op_token[op_loc + 3*i + 2] + # if val1 or val2 are < 0 this means they are to come from the stack + # if token is negative, means we need to use a stack value + #print("s", s) + if t1 < 0: + val1 = s[s_ix] + s_ix -= 1 + else: + val1 = state_ix[t1] + if t2 < 0: + val2 = s[s_ix] + s_ix -= 1 + else: + val2 = state_ix[t2] + #print(s_ix, op, val1, val2) + result = evaluate_eq_ops(op, val1, val2) + s_ix += 1 + if s_ix >= s_len: + s = np.append(s, 0) + s_len += 1 + s[s_ix] = result + result = s[s_ix] + if (non_neg == 1) and (result < 0): + result = state_ix[min_ix] + state_ix[op_token[1]] = result + return True \ No newline at end of file diff --git a/HSP2/om_model_broadcast.py b/HSP2/om_model_broadcast.py new file mode 100644 index 00000000..b19cad2d --- /dev/null +++ b/HSP2/om_model_broadcast.py @@ -0,0 +1,183 @@ +""" +The class Broadcast is used to send and receive data to shared accumulator "channel" and "register". +See also Branch: an actual flow control structure that looks similar to Conditional, but changes execution +""" +from HSP2.state import * +from HSP2.om import * +from HSP2.om_model_object import * +from HSP2.om_model_linkage import ModelLinkage +from numba import njit +import warnings + +class ModelBroadcast(ModelObject): + def __init__(self, name, container = False, model_props = {}): + super(ModelBroadcast, self).__init__(name, container, model_props) + self.model_props_parsed = model_props + # broadcast_params = [ [local_name1, remote_name1], [local_name2, remote_name2], ...] + # broadcast_channel = state_path/[broadcast_channel] + # broadcast_hub = self, parent, /state/path/to/element/ + # we call handle+=_prop() because this will be OK with any format of caller data + self.linkages = {} # store of objects created by this + self.optype = 4 # 0 - shell object, 1 - equation, 2 - DataMatrix, 3 - input, 4 - broadcastChannel, 5 - ? + self.bc_type_id = 2 # assume read -- is this redundant? is it really the input type ix? + #self.parse_model_props(model_props) + self.setup_broadcast(self.broadcast_type, self.broadcast_params, self.broadcast_channel, self.broadcast_hub) + + + def parse_model_props(self, model_props, strict = False ): + # handle props array + super().parse_model_props(model_props, strict) + self.broadcast_type = self.handle_prop(model_props, 'broadcast_type') + self.broadcast_hub = self.handle_prop(model_props, 'broadcast_hub') + self.broadcast_channel = self.handle_prop(model_props, 'broadcast_channel') + self.broadcast_params = self.handle_prop(model_props, 'broadcast_params') + if self.broadcast_type == None: + self.broadcast_type = 'read' + if self.broadcast_channel == None: + self.broadcast_channel = False + if self.broadcast_hub == None: + self.broadcast_hub = 'self' + if self.broadcast_params == None: + self.broadcast_params = [] + + def handle_prop(self, model_props, prop_name, strict = False, default_value = None ): + # parent method handles most cases, but subclass handles special situations. + if ( prop_name == 'broadcast_params'): + prop_val = model_props.get(prop_name) + #print("broadcast params from model_props = ", prop_val) + if type(prop_val) == list: # this doesn't work, but nothing gets passed in like this? Except broadcast params, but they are handled in the sub-class + prop_val = prop_val + elif type(prop_val) == dict: + prop_val = prop_val.get('value') + else: + prop_val = super().handle_prop(model_props, prop_name, strict, default_value) + return prop_val + + def setup_broadcast(self, broadcast_type, broadcast_params, broadcast_channel, broadcast_hub): + if (broadcast_hub == 'parent') and (self.container.container == False): + warnings.warn("Broadcast named " + self.name + " is parent object " + self.container.name + " with path " + self.container.state_path + " does not have a grand-parent container. Broadcast to hub 'parent' guessing as a path-type global. ") + broadcast_hub = '/STATE/global' + warnings.warn("Proceeding with broadcast_type, broadcast_params, broadcast_channel, broadcast_hub = " + str(broadcast_type) + "," + str(broadcast_params) + "," + str(broadcast_channel) + "," + str(broadcast_hub)) + if ( (broadcast_hub == 'self') or (broadcast_hub == 'child') ): + hub_path = self.container.state_path # the parent of this object is the "self" in question + hub_container = self.container + elif broadcast_hub == 'parent': + hub_path = self.container.container.state_path + hub_container = self.container.container + else: + # we assume this is a valid path. but we verify and fail if it doesn't work during tokenization + # this is not really yet operational since it would be a global broadcast of sorts + print("Broadcast ", self.name, " hub Path not parent, self or child. Trying to find another hub_path = ", broadcast_hub) + hub_path = broadcast_hub + hub_exists = self.find_var_path(hub_path) + if hub_exists == False: + hub_container = False + else: + hub_container = self.model_object_cache[hub_path] + # add the channel to the hub path + channel_path = hub_path + "/" + broadcast_channel + channel = self.insure_channel(broadcast_channel, hub_container) + # now iterate through pairs of source/destination broadcast lines + i = 0 + #print("broadcast params", broadcast_params) + for b_pair in broadcast_params: + # create an object for each element in the array + if (broadcast_type == 'read'): + # a broadcast channel (object has been created above) + # + a register to hold the data (if not yet created) + # + an input on the parent to read the data + src_path = hub_path + "/" + b_pair[1] + self.bc_type_id = 2 # pull type + register_varname = b_pair[0] + # create a link object of type 2, property reader to local state + #print("Adding broadcast read as input from ", channel_path, " as local var named ",register_varname) + # create a register if it does not exist already + var_register = self.insure_register(register_varname, 0.0, channel) + # add input to parent container for this variable from the hub path + self.container.add_input(register_varname, var_register.state_path, 1, True) + # this registers a variable on the parent so the channel is not required to access it + # this is redundant for the purpose of calculation, all model components will know how to find it by tracing inputs + # but this makes it easier for tracking + puller = ModelLinkage(register_varname, self.container, {'right_path':var_register.state_path, 'link_type':2} ) + added = True + else: + # a broadcast hub (if not yet created) + # + a register to hold the data (if not yet created) + # + an input on the broadcast hub to read the data (or a push on this object?) + dest_path = hub_path + "/" + b_pair[1] + self.bc_type_id = 4 # push accumulator type + local_varname = b_pair[0] # this is where we take the data from + register_varname = b_pair[1] # this is the name that will be stored on the register + #print("Adding send from local var ", local_varname, " to ",channel.name) + # create a register as a placeholder for the data at the hub path + # in case there are no readers + var_register = self.insure_register(register_varname, 0.0, channel) + dest_path = var_register.state_path + src_path = self.find_var_path(local_varname) + # this linkage pushes from local value to the remote path + pusher = ModelLinkage(register_varname, self, {'right_path':src_path, 'link_type':self.bc_type_id, 'left_path':dest_path} ) + # try adding the linkage an input, just to see if it influences the ordering + #print("Adding broadcast source ", local_varname, " as input to register ",var_register.name) + # we do an object connection here, as this is a foolproof way to + # add already created objects as inputs + var_register.add_object_input(register_varname + str(pusher.ix), pusher, 1) + # this linkage creates a pull on the remote path, not yet ready since + # there is no bc_type_id that corresponds to an accumulator pull + #puller = ModelLinkage(register_varname, var_register, src_path, self.bc_type_id) + + def insure_channel(self, broadcast_channel, hub_container): + if hub_container == False: + # this is a global, so no container + hub_name = '/STATE/global' + channel_path = False + else: + #print("Looking for channel ", broadcast_channel, " on ", hub_container.name) + # must create absolute path, otherwise, it will seek upstream and get the parent + # we send with local_only = True so it won't go upstream + channel_path = hub_container.find_var_path(broadcast_channel, True) + hub_name = hub_container.name + if channel_path == False: + #print(self.state_path, "is Creating broadcast hub ", broadcast_channel, " on ", hub_name) + hub_object = ModelObject(broadcast_channel, hub_container) + else: + hub_object = self.model_object_cache[channel_path] + return hub_object + + def tokenize(self): + # call parent method to set basic ops common to all + super().tokenize() + # because we added each type as a ModelLinkage, this may be superfluous? + # this ModelLinkage will be handled on it's own. + # are there exec hierarchy challenges because of skipping this? + # exec hierarchy function on object.inputs[] alone. Since broadcasts + # are not treated as inputs, or are they? Do ModelLinkages create inputs? + # - should read inputs create linkages, but send/push linkages not? + # - ex: a facility controls a reservoir release with a push linkage + # the reservoir *should* execute after the facility in this case + # but perhaps that just means we *shouldn't* use a push, instead + # we should have the reservoir pull the info? + # - however, since "parent" pushes will automatically have hierarchy + # preserved, since the child object already is an input to the parent + # and therefore will execute before the parent + # - but, we must insure that ModelLinkages get executed when their container + # is executed (which should come de facto as they are inputs to their container) + #if (self.broadcast_type == 'send'): + # for i in self.linkages.keys(): + # self.ops = self.ops + [self.left_ix, cop_codes[self.cop], self.right_ix] + #else: + # this is a read, so we simply rewrite as a model linkage + # not yet complete or functioning + # self.linkage.tokenize() + + def add_op_tokens(self): + # this puts the tokens into the global simulation queue + # can be customized by subclasses to add multiple lines if needed. + super().add_op_tokens() + + +@njit +def pre_step_broadcast(op, state_ix, dict_ix): + ix = op[1] + dix = op[2] + # broadcasts do not need to act, as they are now handled as ModelLinkage + # with type = accumulate diff --git a/HSP2/om_model_linkage.py b/HSP2/om_model_linkage.py index d4b4ac61..11d5f9a6 100644 --- a/HSP2/om_model_linkage.py +++ b/HSP2/om_model_linkage.py @@ -70,7 +70,7 @@ def find_paths(self): # do we need to do this, or just trust it exists? #self.insure_path(self, self.right_path) # the left path, if this is type 4 or 5, is a push, so we must require it - if ( (self.link_type == 4) or (self.link_type == 5) ): + if ( (self.link_type == 4) or (self.link_type == 5) or (self.link_type == 6) ): self.insure_path(self.left_path) self.paths_found = True return @@ -88,7 +88,7 @@ def tokenize(self): else: print("Error: link ", self.name, "does not have a valid source path") #print(self.name,"tokenize() result", self.ops) - if (self.link_type == 4) or (self.link_type == 5): + if (self.link_type == 4) or (self.link_type == 5) or (self.link_type == 6): # we push to the remote path in this one left_ix = get_state_ix(self.state_ix, self.state_paths, self.left_path) right_ix = get_state_ix(self.state_ix, self.state_paths, self.right_path) @@ -107,6 +107,7 @@ def step_model_link(op_token, state_ix, ts_ix, step): return True elif op_token[3] == 2: state_ix[op_token[1]] = state_ix[op_token[2]] + return True elif op_token[3] == 3: # read from ts variable TBD # state_ix[op_token[1]] = ts_ix[op_token[2]][step] @@ -117,29 +118,9 @@ def step_model_link(op_token, state_ix, ts_ix, step): return True elif op_token[3] == 5: # overwrite remote variable state with value in another paths state - if step == 2: - print("Setting state_ix[", op_token[2], "] =", state_ix[op_token[4]]) state_ix[op_token[2]] = state_ix[op_token[4]] return True - - -def test_model_link(op_token, state_ix, ts_ix, step): - if op_token[3] == 1: - return True - elif op_token[3] == 2: - state_ix[op_token[1]] = state_ix[op_token[2]] - elif op_token[3] == 3: - # read from ts variable TBD - # state_ix[op_token[1]] = ts_ix[op_token[2]][step] - return True - elif op_token[3] == 4: - print("Remote Broadcast accumulator type link.") - print("Setting op ID", str(op_token[2]), "to value from ID", str(op_token[4]), "with value of ") - # add value in local state to the remote broadcast hub+register state - state_ix[op_token[2]] = state_ix[op_token[2]] + state_ix[op_token[4]] - print(str(state_ix[op_token[2]]) + " = ", str(state_ix[op_token[2]]) + "+" + str(state_ix[op_token[4]])) + elif op_token[3] == 6: + # set value in a timerseries + ts_ix[op_token[2]][step] = state_ix[op_token[4]] return True - elif op_token[3] == 5: - # push value in local state to the remote broadcast hub+register state - state_ix[op_token[2]] = state_ix[op_token[4]] - return True \ No newline at end of file diff --git a/HSP2/om_model_object.py b/HSP2/om_model_object.py index 0a13a95c..c1404b12 100644 --- a/HSP2/om_model_object.py +++ b/HSP2/om_model_object.py @@ -322,6 +322,19 @@ def find_paths(self): # and should also handle deciding if this is a constant, like a numeric value # or a variable data and should handle them accordingly return True + + def insure_register(self, var_name, default_value, register_container, register_path = False): + # we send with local_only = True so it won't go upstream + if register_path == False: + register_path = register_container.find_var_path(var_name, True) + if (register_path == False) or (register_path not in self.model_object_cache.keys()): + # create a register as a placeholder for the data at the hub path + # in case there are no senders, or in the case of a timeseries logger, we need to register it so that its path can be set to hold data + #print("Creating a register for data for hub ", register_container.name, "(", register_container.state_path, ")", " var name ",var_name) + var_register = ModelRegister(var_name, register_container, default_value, register_path) + else: + var_register = self.model_object_cache[register_path] + return var_register def tokenize(self): # renders tokens for high speed execution @@ -382,9 +395,46 @@ def required_properties(): req_props.extend(['value']) return req_props + +""" +The class ModelRegister is for storing push values. +Behavior is to zero each timestep. This could be amended later. +Maybe combined with stack behavior? Or accumulator? +""" +class ModelRegister(ModelConstant): + def __init__(self, name, container = False, value = 0.0, state_path = False): + super(ModelRegister, self).__init__(name, container, value, state_path) + self.optype = 12 # + # self.state_ix[self.ix] = self.default_value + + def required_properties(): + req_props = super(ModelConstant, ModelConstant).required_properties() + req_props.extend(['value']) + return req_props + # njit functions for runtime +@njit +def pre_step_register(op, state_ix): + ix = op[1] + #print("Resetting register", ix,"to zero") + state_ix[ix] = 0.0 + return + +# Note: ModelConstant has not runtime execution @njit def exec_model_object( op, state_ix, dict_ix): ix = op[1] - return 0.0 \ No newline at end of file + return 0.0 + + +# njit functions for end of model run +@njit +def finish_model_object(op_token, state_ix, ts_ix): + return + + +@njit +def finish_register(op_token, state_ix, ts_ix): + # todo: push the values of ts_ix back to the hdf5? or does this happen in larger simulation as it is external to OM? + return \ No newline at end of file diff --git a/tests/testcbp/HSP2results/benchmark.py b/tests/testcbp/HSP2results/benchmark.py index 79387fba..9ed1e676 100644 --- a/tests/testcbp/HSP2results/benchmark.py +++ b/tests/testcbp/HSP2results/benchmark.py @@ -67,5 +67,5 @@ start = time.time() iterate_models(model_exec_list, op_tokens, np_state_ix, dict_ix, ts_ix, siminfo['steps'], -1) end = time.time() -print(len(model_exec_list), "components iterated over np_state_ix", siminfo['steps'], "time steps took" , end - start, "seconds") +print(len(model_exec_list), "components iterated over np_state_ix(", type(np_state_ix),")", siminfo['steps'], "time steps took" , end - start, "seconds") diff --git a/tests/testcbp/HSP2results/benchmark_equation.py b/tests/testcbp/HSP2results/benchmark_equation.py new file mode 100644 index 00000000..6298c522 --- /dev/null +++ b/tests/testcbp/HSP2results/benchmark_equation.py @@ -0,0 +1,110 @@ +# bare bones tester +# tests special actions and constants. +import os +# disabled for auto testing, but may use at command prompt if needed +#os.chdir("C:/usr/local/home/git/HSPsquared") +from HSP2.main import * +from HSP2.om import * +#from HSP2.om_equation import * +import pytest + +state = init_state_dicts() +# set up info and timer +siminfo = {} +siminfo['delt'] = 60 +siminfo['tindex'] = date_range("1984-01-01", "2020-12-31", freq=Minute(siminfo['delt']))[1:] +steps = siminfo['steps'] = len(siminfo['tindex']) +# get any pre-loaded objects +model_data = state['model_data'] +( ModelObject.op_tokens, ModelObject.model_object_cache) = init_om_dicts() +ModelObject.state_paths, ModelObject.state_ix, ModelObject.dict_ix, ModelObject.ts_ix = state['state_paths'], state['state_ix'], state['dict_ix'], state['ts_ix'] +( op_tokens, state_paths, state_ix, dict_ix, model_object_cache, ts_ix) = ( ModelObject.op_tokens, ModelObject.state_paths, ModelObject.state_ix, ModelObject.dict_ix, ModelObject.model_object_cache, ModelObject.ts_ix ) +state_context_hsp2(state, 'RCHRES', 'R001', 'HYDR') +print("Init HYDR state context for domain", state['domain']) +hydr_init_ix(state['state_ix'], state['state_paths'], state['domain']) +# Now, assemble a test dataset +container = False +state_om_model_root_object(state, siminfo) +model_root_object = state['model_root_object'] + +facility = ModelObject('facility', model_root_object) +for k in range(1000): + eqn = str(25*random.random()) + " * " + c[round((2*random.random()))] + newq = Equation('eq' + str(k), facility, {'equation':eqn} ) + conval = 50.0*random.random() + newc = ModelConstant('con' + str(k), facility, conval) + +# create a register to test TS +ts1 = facility.insure_register('/TIMESERIES/facility/con1', 0.0, facility) +# do all linking and tokenizing, 2nd arg "io_manager" is False as we do not have an hdf5 here +# set ops_data_type = Dict to test the Dict performance for state_ix +# override ops_data_type for testing: +ModelObject.ops_data_type = 'Dict' +# - this forces state_om_model_run_prep() to use Dict type for op_tokens and state_ix +# - this can also be overridden by setting "siminfo" : {"ops_data_type" : "Dict" } in the json file +# now initialize model data sets (sorts, tokenizes, etc.) +state_om_model_run_prep(state, False, siminfo) +op_tokens = state['op_tokens'] + +# run 1 time to compile all if anything is changed +model_exec_list = state['model_exec_list'] +iterate_models(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, 1, -1) + +# test with np.array state_ix +#np_state_ix = np.asarray(list(state_ix.values()), dtype="float32") +np_state_ix = zeros(max(state_ix.keys()) + 1, dtype="float32") +# this insures that the keys in the nparray version of state match +for ix, iv in state_ix.items(): + np_state_ix[ix] = iv + + +# Test and time the run with Dict version of state_ix +start = time.time() +iterate_models(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, siminfo['steps'], -1) +end = time.time() +print(len(model_exec_list), "state_ix components iterated with full execution via iterate_models()", siminfo['steps'], "time steps took" , end - start, "seconds") + +start = time.time() +iterate_models(model_exec_list, op_tokens, np_state_ix, dict_ix, ts_ix, siminfo['steps'], -1) +end = time.time() +print(len(model_exec_list), "np_state_ix components iterated with full execution via iterate_models()", siminfo['steps'], "time steps took" , end - start, "seconds") + + +@njit +def iteration_test_dat(op_order, it_ops, state_ix, it_nums): + ctr = 0 + ttr = 0.0 + for n in range(it_nums): + for i in op_order: + op = it_ops[i] + getsx = state_ix[i] + state_ix[i] = getsx * 1.0 + #ttr = ttr + getsx + ctr = ctr + 1 + print("Completed ", ctr, " loops") + + +# Now test just the data structures with no actual calculation from primitives +start = time.time() +iteration_test_dat(model_exec_list, op_tokens, state_ix, siminfo['steps']) +end = time.time() +print(len(model_exec_list), "components iterated over test data setter state_ix", siminfo['steps'], "time steps took" , end - start, "seconds") + +start = time.time() +iteration_test_dat(model_exec_list, op_tokens, np_state_ix, siminfo['steps']) +end = time.time() +print(len(model_exec_list), "components iterated over test data setter np_state_ix", siminfo['steps'], "time steps took" , end - start, "seconds") + +start = time.time() +iterate_perf(model_exec_list, op_tokens, state_ix, dict_ix, ts_ix, siminfo['steps']) +end = time.time() +print(len(model_exec_list), "state_ix components iterated over iterate_perf (selective component) ", siminfo['steps'], "time steps took" , end - start, "seconds") + + +start = time.time() +iterate_perf(model_exec_list, op_tokens, np_state_ix, dict_ix, ts_ix, siminfo['steps']) +end = time.time() +print(len(model_exec_list), "np_state_ix components iterated over iterate_perf (selective component)", siminfo['steps'], "time steps took" , end - start, "seconds") + +def test_benchmark(model_exec_list, start, end): + print(len(model_exec_list), "np_state_ix components iterated over iterate_perf (selective component)", siminfo['steps'], "time steps took" , end - start, "seconds") diff --git a/tests/testcbp/HSP2results/example_manual_object.py b/tests/testcbp/HSP2results/example_manual_object.py index 021c10a8..f738f015 100644 --- a/tests/testcbp/HSP2results/example_manual_object.py +++ b/tests/testcbp/HSP2results/example_manual_object.py @@ -59,4 +59,3 @@ O1 = ModelLinkage('O1', hydr, wd_mgd.state_path, 2) O1.add_op_tokens() -return