From ef260e891ac91edc4abd4833105318f0d3715d24 Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Tue, 6 Feb 2024 14:53:06 -0500 Subject: [PATCH 1/8] Add support for multiple inputs,ouputs,controllers to a flow and parse nested aux variables, assign place holder functions for expressions with functions --- mira/sources/system_dynamics/pysd.py | 52 +++++++++++++++----------- mira/sources/system_dynamics/stella.py | 24 ++++++++++-- 2 files changed, 51 insertions(+), 25 deletions(-) diff --git a/mira/sources/system_dynamics/pysd.py b/mira/sources/system_dynamics/pysd.py index c2732b07e..d910948c8 100644 --- a/mira/sources/system_dynamics/pysd.py +++ b/mira/sources/system_dynamics/pysd.py @@ -21,7 +21,7 @@ sympy.Symbol("Day"): sympy.Symbol("day"), sympy.Symbol("Days"): sympy.Symbol("day"), } - +SYMPY_FLOW_RATE_PLACEHOLDER = safe_parse_expr("XXplaceholderXX") __all__ = ["template_model_from_pysd_model"] @@ -139,7 +139,10 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # process parameters mira_parameters = {} for name, expression in processed_expression_map.items(): - if expression.replace(".", "").replace(" ", "").isdecimal(): + if ( + expression + and expression.replace(".", "").replace(" ", "").isdecimal() + ): model_parameter_info = model_doc_df[model_doc_df["Py Name"] == name] if model_parameter_info["Units"].values[0]: unit_text = ( @@ -159,6 +162,8 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: "description": model_parameter_info["Comment"].values[0], } + if name == "baseline_md_mortality_rate": + pass mira_parameters[name] = parameter_to_mira(parameter) # standardize parameter units if they exist @@ -179,18 +184,19 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: ): rate_name = aux_tuple["Py Name"] rate_expr = safe_parse_expr( - preprocess_text(processed_expression_map[rate_name]), + processed_expression_map[rate_name], symbols, ) - input_state, output_state, controller = None, None, None + + inputs, outputs, controllers = [], [], [] # If we come across a rate-law that is leaving a state, we add the state as an input # to the rate-law, vice-versa if a rate-law is going into a state. for state_name, in_out in state_rate_map.items(): if rate_name in in_out["output_rates"]: - input_state = state_name + inputs.append(state_name) if rate_name in in_out["input_rates"]: - output_state = state_name + outputs.append(state_name) # if a state isn't consumed by a flow (the flow isn't listed as an output of # the state) but affects the rate of a flow, then that state is a controller if ( @@ -198,14 +204,14 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: and rate_name not in state_rate_map[state_name]["output_rates"] ): - controller = output_state + controllers.append(state_name) transition_map[rate_name] = { "name": rate_name, "expression": rate_expr, - "input": input_state, - "output": output_state, - "controller": controller, + "inputs": inputs, + "outputs": outputs, + "controllers": controllers, } used_states = set() @@ -215,30 +221,32 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: for template_id, (transition_name, transition) in enumerate( transition_map.items() ): - input_concepts = [] - output_concepts = [] + input_concepts, input_names = [], [] + output_concepts, output_names = [], [] controller_concepts = [] - input_name, output_name, controller_name = None, None, None - if transition.get("input"): - input_name = transition.get("input") + for input_name in transition.get("inputs"): input_concepts.append(concepts[input_name].copy(deep=True)) - if transition.get("output"): - output_name = transition.get("output") + input_names.append(input_name) + + for output_name in transition.get("outputs"): output_concepts.append(concepts[output_name].copy(deep=True)) - if transition.get("controller"): - controller_name = transition.get("controller") + output_names.append(output_name) + + for controller_name in transition.get("controllers"): controller_concepts.append( concepts[controller_name].copy(deep=True) ) - used_states |= {input_name, output_name} + used_states.update(input_names, output_names) templates_.extend( transition_to_templates( input_concepts=input_concepts, output_concepts=output_concepts, controller_concepts=controller_concepts, - transition_rate=transition["expression"], + transition_rate=transition["expression"] + if transition["expression"] != SYMPY_FLOW_RATE_PLACEHOLDER + else None, transition_id=str(template_id + 1), transition_name=transition_name, ) @@ -293,6 +301,8 @@ def preprocess_text(expr_text): # strip leading and trailing white spaces # remove spaces between operators and operands # replace space between two words that makeup a variable name with "_"' + if not expr_text: + return expr_text expr_text = ( expr_text.strip() .replace(" * ", "*") diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 99a55a081..2b74bdd83 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -52,7 +52,6 @@ def template_model_from_stella_model_file(fname) -> TemplateModel: : A MIRA template model """ - pysd_model = pysd.read_xmile(fname) stella_model_file = XmileFile(fname) stella_model_file.parse() @@ -114,6 +113,11 @@ def extract_stella_variable_expressions(stella_model_file): # test to see if flow first as flow is a subclass or Aux elif isinstance(component, Flow): operands = component.components[0][1].arguments + # If the flow doesn't use operands (e.g. +,-) but actual functions (e.g. max,pulse) + # assign the flow as None + if not hasattr(component.components[0][1], "operators"): + expression_map[component.name] = "XXplaceholderXX" + continue operators = component.components[0][1].operators operand_operator_per_level_var_map[component.name] = {} extract_variables( @@ -123,12 +127,10 @@ def extract_stella_variable_expressions(stella_model_file): operand_operator_per_level_var_map, ) elif isinstance(component, Aux): - expression_map[component.name] = str(component.components[0][1]) + expression_map[component.name] = parse_aux_structure(component) elif isinstance(component, Stock): try: operand_operator_per_level_var_map[component.name] = {} - if component.name == "recovered": - pass operands = component.components[0][1].flow.arguments operators = component.components[0][1].flow.operators extract_variables( @@ -265,3 +267,17 @@ def parse_structures(operand, idx, name, operand_operator_per_level_var_map): operand_operator_per_level_var_map[name][idx + 1]["operators"].extend( operand.operators ) + + +def parse_aux_structure(structure): + val = structure.components[0][1] + while ( + not isinstance(val, float) + and not isinstance(val, str) + and not isinstance(val, int) + ): + if hasattr(val, "argument"): + val = val.argument + if isinstance(val, ReferenceStructure): + val = val.reference + return str(val) From 8bb8a5b917cb50fb8788c13297bc692e2889c2ac Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Mon, 12 Feb 2024 18:06:34 -0500 Subject: [PATCH 2/8] Add preprocessing for expressions that cannot be parsed by pyd in raw stella model files, improve creation of transition map, add support for dimensionless units, modify placeholder used --- mira/sources/system_dynamics/pysd.py | 55 ++++++++++++---- mira/sources/system_dynamics/stella.py | 88 +++++++++++++++++++++++++- 2 files changed, 127 insertions(+), 16 deletions(-) diff --git a/mira/sources/system_dynamics/pysd.py b/mira/sources/system_dynamics/pysd.py index d910948c8..34ed9fa13 100644 --- a/mira/sources/system_dynamics/pysd.py +++ b/mira/sources/system_dynamics/pysd.py @@ -1,6 +1,7 @@ """This module implements parsing of a generic pysd model irrespective of source and source type and extracting its contents to create an equivalent MIRA template model. """ +__all__ = ["template_model_from_pysd_model"] import pandas as pd import sympy @@ -14,15 +15,22 @@ get_sympy, ) -CONTROL_VARIABLE_NAMES = {"FINALTIME", "INITIALTIME", "SAVEPER", "TIMESTEP"} +CONTROL_VARIABLE_NAMES = { + "FINALTIME", + "INITIALTIME", + "SAVEPER", + "TIMESTEP", + "FINAL TIME", + "INITIAL TIME", + "TIME STEP", +} UNITS_MAPPING = { sympy.Symbol("Person"): sympy.Symbol("person"), sympy.Symbol("Persons"): sympy.Symbol("person"), sympy.Symbol("Day"): sympy.Symbol("day"), sympy.Symbol("Days"): sympy.Symbol("day"), } -SYMPY_FLOW_RATE_PLACEHOLDER = safe_parse_expr("XXplaceholderXX") -__all__ = ["template_model_from_pysd_model"] +SYMPY_FLOW_RATE_PLACEHOLDER = safe_parse_expr("xxplaceholderxx") def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: @@ -139,18 +147,27 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # process parameters mira_parameters = {} for name, expression in processed_expression_map.items(): + eval_expression = safe_parse_expr(expression).evalf() + str_eval_expression = str(eval_expression) + if ( - expression - and expression.replace(".", "").replace(" ", "").isdecimal() + eval_expression != SYMPY_FLOW_RATE_PLACEHOLDER + and str_eval_expression.replace(".", "") + .replace(" ", "") + .isdecimal() ): model_parameter_info = model_doc_df[model_doc_df["Py Name"] == name] - if model_parameter_info["Units"].values[0]: + if ( + model_parameter_info["Units"].values[0] + and model_parameter_info["Units"].values[0] != "dimensionless" + ): unit_text = ( model_parameter_info["Units"].values[0].replace(" ", "") ) + parameter = { "id": name, - "value": float(expression), + "value": float(str_eval_expression), "description": model_parameter_info["Comment"].values[0], "units": {"expression": unit_text}, } @@ -158,12 +175,10 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # if units don't exist parameter = { "id": name, - "value": float(expression), + "value": float(str_eval_expression), "description": model_parameter_info["Comment"].values[0], } - if name == "baseline_md_mortality_rate": - pass mira_parameters[name] = parameter_to_mira(parameter) # standardize parameter units if they exist @@ -176,7 +191,13 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # construct transitions mapping that determine inputs and outputs states to a rate-law transition_map = {} - auxiliaries = model_doc_df[model_doc_df["Type"] == "Auxiliary"] + auxiliaries = model_doc_df[ + (model_doc_df["Type"] == "Auxiliary") + | (model_doc_df["Type"] == "Constant") + ] + + # currently, we add every auxiliary to the map of transitions even if it is not a transition + # no set way to differentiate between auxiliaries of transitions for index, aux_tuple in auxiliaries.iterrows(): if ( aux_tuple["Subtype"] == "Normal" @@ -192,10 +213,10 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # If we come across a rate-law that is leaving a state, we add the state as an input # to the rate-law, vice-versa if a rate-law is going into a state. - for state_name, in_out in state_rate_map.items(): - if rate_name in in_out["output_rates"]: + for state_name, in_out_rate_map in state_rate_map.items(): + if rate_name in in_out_rate_map["output_rates"]: inputs.append(state_name) - if rate_name in in_out["input_rates"]: + if rate_name in in_out_rate_map["input_rates"]: outputs.append(state_name) # if a state isn't consumed by a flow (the flow isn't listed as an output of # the state) but affects the rate of a flow, then that state is a controller @@ -206,6 +227,11 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: ): controllers.append(state_name) + # if the auxiliary does not have inputs, outputs, or controllers, we know it is not + # a transition + if not inputs and not outputs and not controllers: + continue + transition_map[rate_name] = { "name": rate_name, "expression": rate_expr, @@ -224,6 +250,7 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: input_concepts, input_names = [], [] output_concepts, output_names = [], [] controller_concepts = [] + for input_name in transition.get("inputs"): input_concepts.append(concepts[input_name].copy(deep=True)) input_names.append(input_name) diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 2b74bdd83..4c54b6912 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -18,6 +18,7 @@ import tempfile from pathlib import Path +import re import pysd from pysd.translators.xmile.xmile_file import XmileFile @@ -38,6 +39,54 @@ template_model_from_pysd_model, ) +PLACE_HOLDER_EQN = "0" + + +def process_file(fname): + with open(fname) as f: + xml_str = f.read() + + # commented out for current debugging + # xml_str = xml_str.replace("\n","") + + # Conditional preprocessing + eqn_tags = re.findall(r".*?", xml_str, re.DOTALL) + for tag in eqn_tags: + if all(word in tag.lower() for word in ["if", "then", "else"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if all(word in tag.lower() for word in ["int", "mod"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "sum" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "//" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "," in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "init" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "nan" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "pi" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + + # processes spaces in eqns + eqn_space_pattern = r"(.*?)" + xml_str = re.sub( + eqn_space_pattern, + lambda m: f'{m.group(1).replace(" ", "")}', + xml_str, + ) + + # Comment preprocessing + eqn_bracket_pattern = r"(.*?)\{[^}]*\}(.*?)" + eqn_bracket_sub = r"\1" + r"\2" + xml_str = re.sub(eqn_bracket_pattern, eqn_bracket_sub, xml_str) + + with open(fname, "w") as f: + f.write(xml_str) + + return fname + def template_model_from_stella_model_file(fname) -> TemplateModel: """Return a template model from a local Stella model file. @@ -52,6 +101,9 @@ def template_model_from_stella_model_file(fname) -> TemplateModel: : A MIRA template model """ + + process_file(fname) + pysd_model = pysd.read_xmile(fname) stella_model_file = XmileFile(fname) stella_model_file.parse() @@ -83,6 +135,7 @@ def template_model_from_stella_model_url(url) -> TemplateModel: with temp_file as file: file.write(data) + process_file(temp_file.name) pysd_model = pysd.read_xmile(temp_file.name) stella_model_file = XmileFile(temp_file.name) stella_model_file.parse() @@ -112,11 +165,12 @@ def extract_stella_variable_expressions(stella_model_file): continue # test to see if flow first as flow is a subclass or Aux elif isinstance(component, Flow): + # TODO: test for call structure object operands = component.components[0][1].arguments # If the flow doesn't use operands (e.g. +,-) but actual functions (e.g. max,pulse) # assign the flow as None if not hasattr(component.components[0][1], "operators"): - expression_map[component.name] = "XXplaceholderXX" + expression_map[component.name] = "xxplaceholderxx" continue operators = component.components[0][1].operators operand_operator_per_level_var_map[component.name] = {} @@ -127,7 +181,20 @@ def extract_stella_variable_expressions(stella_model_file): operand_operator_per_level_var_map, ) elif isinstance(component, Aux): - expression_map[component.name] = parse_aux_structure(component) + # TODO: test for call structure object + if not isinstance(component.components[0][1], ArithmeticStructure): + expression_map[component.name] = parse_aux_structure(component) + continue + operands = component.components[0][1].arguments + operators = component.components[0][1].operators + operand_operator_per_level_var_map[component.name] = {} + extract_variables( + operands, + operators, + component.name, + operand_operator_per_level_var_map, + ) + elif isinstance(component, Stock): try: operand_operator_per_level_var_map[component.name] = {} @@ -268,8 +335,25 @@ def parse_structures(operand, idx, name, operand_operator_per_level_var_map): operand.operators ) + # case where operand is just a primitive + else: + operand_operator_per_level_var_map[name][idx]["operands"].append( + str(operand) + ) + def parse_aux_structure(structure): + """If an Aux object's component method is not a reference structure, deconstruct the Aux object + to retrieve the primitive value that the auxiliary represents. + + Parameters + ---------- + structure : Aux + The Aux object + Returns + ------- + + """ val = structure.components[0][1] while ( not isinstance(val, float) From a97373a6466af40cddc3434b4b95d2f34d6d8ba5 Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Mon, 12 Feb 2024 20:09:05 -0500 Subject: [PATCH 3/8] Remove removal of spaces in eqns, add support for parameter values that reference a stock --- mira/sources/system_dynamics/pysd.py | 15 +++++++++++---- mira/sources/system_dynamics/stella.py | 9 +-------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/mira/sources/system_dynamics/pysd.py b/mira/sources/system_dynamics/pysd.py index 34ed9fa13..f1d39f051 100644 --- a/mira/sources/system_dynamics/pysd.py +++ b/mira/sources/system_dynamics/pysd.py @@ -147,15 +147,22 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # process parameters mira_parameters = {} for name, expression in processed_expression_map.items(): + # sometimes parameter values reference a stock rather than being a number eval_expression = safe_parse_expr(expression).evalf() str_eval_expression = str(eval_expression) - - if ( + value = None + is_initial = False + if str_eval_expression in mira_initials: + value = float(str(mira_initials[str_eval_expression].expression)) + is_initial = True + if str_eval_expression in mira_initials or ( eval_expression != SYMPY_FLOW_RATE_PLACEHOLDER and str_eval_expression.replace(".", "") .replace(" ", "") .isdecimal() ): + if not is_initial: + value = float(str_eval_expression) model_parameter_info = model_doc_df[model_doc_df["Py Name"] == name] if ( model_parameter_info["Units"].values[0] @@ -167,7 +174,7 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: parameter = { "id": name, - "value": float(str_eval_expression), + "value": value, "description": model_parameter_info["Comment"].values[0], "units": {"expression": unit_text}, } @@ -175,7 +182,7 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # if units don't exist parameter = { "id": name, - "value": float(str_eval_expression), + "value": value, "description": model_parameter_info["Comment"].values[0], } diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 4c54b6912..31045401c 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -14,6 +14,7 @@ __all__ = [ "template_model_from_stella_model_file", "template_model_from_stella_model_url", + "process_file" ] import tempfile @@ -69,14 +70,6 @@ def process_file(fname): if "pi" in tag.lower(): xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - # processes spaces in eqns - eqn_space_pattern = r"(.*?)" - xml_str = re.sub( - eqn_space_pattern, - lambda m: f'{m.group(1).replace(" ", "")}', - xml_str, - ) - # Comment preprocessing eqn_bracket_pattern = r"(.*?)\{[^}]*\}(.*?)" eqn_bracket_sub = r"\1" + r"\2" From 06e99eda440ba6393b0135f34ab04e32d419faf2 Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Tue, 13 Feb 2024 09:51:01 -0500 Subject: [PATCH 4/8] Add handling for initial values that are arrays, and handling for different types of components for Stella expression gathering --- mira/sources/system_dynamics/pysd.py | 42 ++++++++++++++++++++------ mira/sources/system_dynamics/stella.py | 35 ++++++++++++++++----- 2 files changed, 59 insertions(+), 18 deletions(-) diff --git a/mira/sources/system_dynamics/pysd.py b/mira/sources/system_dynamics/pysd.py index f1d39f051..a10dee03b 100644 --- a/mira/sources/system_dynamics/pysd.py +++ b/mira/sources/system_dynamics/pysd.py @@ -138,27 +138,43 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: for state_initial_value, (state_name, state_concept) in zip( state_initial_values, concepts.items() ): - initial = Initial( - concept=concepts[state_name].copy(deep=True), - expression=SympyExprStr(sympy.Float(state_initial_value)), - ) + # for the case when a state's initial value is a numpy array + try: + initial = Initial( + concept=concepts[state_name].copy(deep=True), + expression=SympyExprStr(sympy.Float(state_initial_value)), + ) + except TypeError: + initial = Initial( + concept=concepts[state_name].copy(deep=True), + expression=SympyExprStr("0") + ) mira_initials[initial.concept.name] = initial # process parameters mira_parameters = {} for name, expression in processed_expression_map.items(): - # sometimes parameter values reference a stock rather than being a number - eval_expression = safe_parse_expr(expression).evalf() + # Sometimes parameter values reference a stock rather than being a number + # Current placeholder for incorrectly constructed parameter expressions + try: + eval_expression = safe_parse_expr(expression).evalf() + except TypeError: + eval_expression = safe_parse_expr("0") + str_eval_expression = str(eval_expression) value = None is_initial = False if str_eval_expression in mira_initials: value = float(str(mira_initials[str_eval_expression].expression)) is_initial = True + + # Replace negative signs for placeholder parameter values for Aux structures + # that cannot be parsed if str_eval_expression in mira_initials or ( eval_expression != SYMPY_FLOW_RATE_PLACEHOLDER and str_eval_expression.replace(".", "") .replace(" ", "") + .replace("-", "") .isdecimal() ): if not is_initial: @@ -211,10 +227,14 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: and aux_tuple["Real Name"] not in CONTROL_VARIABLE_NAMES ): rate_name = aux_tuple["Py Name"] - rate_expr = safe_parse_expr( - processed_expression_map[rate_name], - symbols, - ) + # Current placeholder for incorrectly constructed rate/parameter expressions + try: + rate_expr = safe_parse_expr( + processed_expression_map[rate_name], + symbols, + ) + except TypeError: + rate_expr = SYMPY_FLOW_RATE_PLACEHOLDER inputs, outputs, controllers = [], [], [] @@ -236,6 +256,8 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: # if the auxiliary does not have inputs, outputs, or controllers, we know it is not # a transition + # some flows also used as auxiliaries for other flows' expression rates and thus are + # not converted to templates if not inputs and not outputs and not controllers: continue diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 31045401c..91bff30c4 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -14,7 +14,7 @@ __all__ = [ "template_model_from_stella_model_file", "template_model_from_stella_model_url", - "process_file" + "process_file", ] import tempfile @@ -32,6 +32,7 @@ from pysd.translators.structures.abstract_expressions import ( ArithmeticStructure, ReferenceStructure, + InlineLookupsStructure, ) import requests @@ -48,7 +49,7 @@ def process_file(fname): xml_str = f.read() # commented out for current debugging - # xml_str = xml_str.replace("\n","") + xml_str = xml_str.replace("\n", "") # Conditional preprocessing eqn_tags = re.findall(r".*?", xml_str, re.DOTALL) @@ -159,6 +160,10 @@ def extract_stella_variable_expressions(stella_model_file): # test to see if flow first as flow is a subclass or Aux elif isinstance(component, Flow): # TODO: test for call structure object + # if primitive represents flow rate expression + if not hasattr(component.components[0][1], "arguments"): + expression_map[component.name] = str(component.components[0][1]) + continue operands = component.components[0][1].arguments # If the flow doesn't use operands (e.g. +,-) but actual functions (e.g. max,pulse) # assign the flow as None @@ -175,7 +180,15 @@ def extract_stella_variable_expressions(stella_model_file): ) elif isinstance(component, Aux): # TODO: test for call structure object + # TODO: Add more comprehensive coverage for different expression structures for Aux + # and Flows if not isinstance(component.components[0][1], ArithmeticStructure): + if isinstance( + component.components[0][1], InlineLookupsStructure + ): + # current place holder value + expression_map[component.name] = "-1" + continue expression_map[component.name] = parse_aux_structure(component) continue operands = component.components[0][1].arguments @@ -201,14 +214,21 @@ def extract_stella_variable_expressions(stella_model_file): ) # If the stock only has a reference and no operators in its expression except AttributeError: - expression_map[component.name] = component.components[0][ - 1 - ].flow.reference + stock_flow = component.components[0][1].flow + if hasattr(stock_flow, "reference"): + expression_map[component.name] = stock_flow.reference + else: + expression_map[component.name] = stock_flow + operand_operator_per_level_var_map[component.name] = {} operand_operator_per_level_var_map[component.name][0] = {} operand_operator_per_level_var_map[component.name][0][ "operands" - ] = [component.components[0][1].flow.reference] + ] = ( + [stock_flow.reference] + if hasattr(stock_flow, "reference") + else [stock_flow] + ) # construct the expression for each variable once its operators and operands are mapped for var_name, expr_level_dict in operand_operator_per_level_var_map.items(): @@ -234,7 +254,7 @@ def create_expression(expr_level_dict): operators = ops["operators"] if ops.get("operators") else [] if not operators: level_expr = operands[0] - str_expression += level_expr + str_expression += str(level_expr) continue elif len(operands) > len(operators): level_expr = "(" @@ -318,7 +338,6 @@ def parse_structures(operand, idx, name, operand_operator_per_level_var_map): operand_operator_per_level_var_map[name][idx]["operands"].append( operand.reference ) - return elif isinstance(operand, ArithmeticStructure): for struct in operand.arguments: parse_structures( From 8568c5d1366dcdf9abc4704bcad562a3a803540a Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Tue, 13 Feb 2024 17:03:06 -0500 Subject: [PATCH 5/8] Remove backslash from variable names which intefers with sympy parsing --- mira/sources/system_dynamics/pysd.py | 6 ++++-- mira/sources/system_dynamics/stella.py | 18 +++++++++++++----- 2 files changed, 17 insertions(+), 7 deletions(-) diff --git a/mira/sources/system_dynamics/pysd.py b/mira/sources/system_dynamics/pysd.py index a10dee03b..5be3f4f22 100644 --- a/mira/sources/system_dynamics/pysd.py +++ b/mira/sources/system_dynamics/pysd.py @@ -152,6 +152,9 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: mira_initials[initial.concept.name] = initial # process parameters + # Currently cannot capture all parameters as some of them cannot have their expressions parsed + # Additionally, some auxiliary variables are added as parameters due to inability to parse + # the auxiliary expression and differentiate it from a parameter mira_parameters = {} for name, expression in processed_expression_map.items(): # Sometimes parameter values reference a stock rather than being a number @@ -168,7 +171,7 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: value = float(str(mira_initials[str_eval_expression].expression)) is_initial = True - # Replace negative signs for placeholder parameter values for Aux structures + # Replace negative signs for placeholder parameter value for Aux structures # that cannot be parsed if str_eval_expression in mira_initials or ( eval_expression != SYMPY_FLOW_RATE_PLACEHOLDER @@ -211,7 +214,6 @@ def template_model_from_pysd_model(pysd_model, expression_map) -> TemplateModel: param_unit.expression = param_unit.expression.subs( old_unit_symbol, new_unit_symbol ) - # construct transitions mapping that determine inputs and outputs states to a rate-law transition_map = {} auxiliaries = model_doc_df[ diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 91bff30c4..8dbab3e11 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -56,7 +56,9 @@ def process_file(fname): for tag in eqn_tags: if all(word in tag.lower() for word in ["if", "then", "else"]): xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if all(word in tag.lower() for word in ["int", "mod"]): + if all(word in tag.lower() for word in ["int", "mod", "(", ")"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if all(word in tag.lower() for word in ["ln", "(", ")"]): xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) if "sum" in tag.lower(): xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) @@ -82,6 +84,10 @@ def process_file(fname): return fname +def replace_backslash(name): + return name.replace("\\\\", "") + + def template_model_from_stella_model_file(fname) -> TemplateModel: """Return a template model from a local Stella model file. @@ -216,7 +222,9 @@ def extract_stella_variable_expressions(stella_model_file): except AttributeError: stock_flow = component.components[0][1].flow if hasattr(stock_flow, "reference"): - expression_map[component.name] = stock_flow.reference + expression_map[component.name] = replace_backslash( + stock_flow.reference + ) else: expression_map[component.name] = stock_flow @@ -225,7 +233,7 @@ def extract_stella_variable_expressions(stella_model_file): operand_operator_per_level_var_map[component.name][0][ "operands" ] = ( - [stock_flow.reference] + [replace_backslash(stock_flow.reference)] if hasattr(stock_flow, "reference") else [stock_flow] ) @@ -336,7 +344,7 @@ def parse_structures(operand, idx, name, operand_operator_per_level_var_map): # base case if isinstance(operand, ReferenceStructure): operand_operator_per_level_var_map[name][idx]["operands"].append( - operand.reference + replace_backslash(operand.reference) ) elif isinstance(operand, ArithmeticStructure): for struct in operand.arguments: @@ -375,5 +383,5 @@ def parse_aux_structure(structure): if hasattr(val, "argument"): val = val.argument if isinstance(val, ReferenceStructure): - val = val.reference + val = replace_backslash(val.reference) return str(val) From d8b3c9df8759241b869fb5898a8808b1fa940128 Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Wed, 14 Feb 2024 14:00:18 -0500 Subject: [PATCH 6/8] Add to stella file docstrings, change methods in Stella API interface to work with xml strings --- mira/sources/system_dynamics/stella.py | 166 +++++++++++++++---------- 1 file changed, 98 insertions(+), 68 deletions(-) diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 8dbab3e11..82a4c24f5 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -1,10 +1,13 @@ """This module implements an API interface for retrieving Stella models by ISEE Systems denoted by the .xmile, .xml, or .stmx extension through a locally downloaded file or URL. We -cannot process stella models with the .itmx extension. Additionally, the pysd library depends on +cannot process stella models with the .itmx extension. Additionally, the PySD library depends on the parsimonious library which fails to parse a number of stella models with valid file -extensions. We then convert the Stella model into a generic pysd model object that will be -parsed and converted to an equivalent MIRA template model. We preprocess the stella model file to -extract variable expressions. +extensions due to incompatible symbols and characters in equations for variables. We implemented +preprocessing for stella models that fixes a number of these parsing +errors when using PySD to ingest these models; however, a number of models still fail to be +parsed by PySD's "read_xmile" method. We extract the contents of the model as a string, +perform expression preprocessing, and then convert the Stella model into a generic pysd model +object that will be converted to an equivalent MIRA template model. Landing page for Stella: https://www.iseesystems.com/store/products/stella-online.aspx @@ -14,11 +17,10 @@ __all__ = [ "template_model_from_stella_model_file", "template_model_from_stella_model_url", - "process_file", + "template_model_stella_model_string", ] import tempfile -from pathlib import Path import re import pysd @@ -44,52 +46,8 @@ PLACE_HOLDER_EQN = "0" -def process_file(fname): - with open(fname) as f: - xml_str = f.read() - - # commented out for current debugging - xml_str = xml_str.replace("\n", "") - - # Conditional preprocessing - eqn_tags = re.findall(r".*?", xml_str, re.DOTALL) - for tag in eqn_tags: - if all(word in tag.lower() for word in ["if", "then", "else"]): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if all(word in tag.lower() for word in ["int", "mod", "(", ")"]): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if all(word in tag.lower() for word in ["ln", "(", ")"]): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "sum" in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "//" in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "," in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "init" in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "nan" in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - if "pi" in tag.lower(): - xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) - - # Comment preprocessing - eqn_bracket_pattern = r"(.*?)\{[^}]*\}(.*?)" - eqn_bracket_sub = r"\1" + r"\2" - xml_str = re.sub(eqn_bracket_pattern, eqn_bracket_sub, xml_str) - - with open(fname, "w") as f: - f.write(xml_str) - - return fname - - -def replace_backslash(name): - return name.replace("\\\\", "") - - def template_model_from_stella_model_file(fname) -> TemplateModel: - """Return a template model from a local Stella model file. + """Return a template model from a local Stella model file Parameters ---------- @@ -101,24 +59,19 @@ def template_model_from_stella_model_file(fname) -> TemplateModel: : A MIRA template model """ + with open(fname) as f: + xml_str = f.read() - process_file(fname) - - pysd_model = pysd.read_xmile(fname) - stella_model_file = XmileFile(fname) - stella_model_file.parse() - - expression_map = extract_stella_variable_expressions(stella_model_file) - return template_model_from_pysd_model(pysd_model, expression_map) + return template_model_stella_model_string(xml_str) def template_model_from_stella_model_url(url) -> TemplateModel: - """Return a template model from a Stella model file provided by an url. + """Return a template model from a Stella model file provided by an url Parameters ---------- url : str - The url to the Stella model file. + The url to the Stella model file Returns ------- @@ -126,16 +79,31 @@ def template_model_from_stella_model_url(url) -> TemplateModel: A MIRA Template Model """ - file_extension = Path(url).suffix - data = requests.get(url).content + xml_str = requests.get(url).content.decode("utf-8") + return template_model_stella_model_string(xml_str) + + +def template_model_stella_model_string(xml_str) -> TemplateModel: + """Returns a semantically equivalent template model to the stella model represented by the + xml string + + Parameters + ---------- + xml_str : The contents of the Stella model + + Returns + ------- + : + A MIRA template model + """ + xml_str = process_string(xml_str) temp_file = tempfile.NamedTemporaryFile( - mode="w+b", suffix=file_extension, delete=False + mode="w+", suffix=".stmx", delete=False ) with temp_file as file: - file.write(data) + file.write(xml_str) - process_file(temp_file.name) pysd_model = pysd.read_xmile(temp_file.name) stella_model_file = XmileFile(temp_file.name) stella_model_file.parse() @@ -245,7 +213,7 @@ def extract_stella_variable_expressions(stella_model_file): def create_expression(expr_level_dict): - """When a variable's operators and operands are mapped, construct the string expression. + """When a variable's operators and operands are mapped, construct the string expression Parameters ---------- @@ -364,7 +332,7 @@ def parse_structures(operand, idx, name, operand_operator_per_level_var_map): def parse_aux_structure(structure): """If an Aux object's component method is not a reference structure, deconstruct the Aux object - to retrieve the primitive value that the auxiliary represents. + to retrieve the primitive value that the auxiliary represents Parameters ---------- @@ -385,3 +353,65 @@ def parse_aux_structure(structure): if isinstance(val, ReferenceStructure): val = replace_backslash(val.reference) return str(val) + + +def process_string(xml_str): + """Helper method that removes incompatible characters from expressions for + variables such that the stella model file can be read by PySD's "read_xmile" method + + Parameters + ---------- + xml_str : str + The xml string representing the contents of the model + + Returns + ------- + : str + The xml string represent the contents of the model after expressions have been + preprocessed + """ + xml_str = xml_str.replace("\n", "") + + eqn_tags = re.findall(r".*?", xml_str, re.DOTALL) + for tag in eqn_tags: + if all(word in tag.lower() for word in ["if", "then", "else"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if all(word in tag.lower() for word in ["int", "mod", "(", ")"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if all(word in tag.lower() for word in ["ln", "(", ")"]): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "sum" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "//" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "," in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "init" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "nan" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + if "pi" in tag.lower(): + xml_str = xml_str.replace(tag, PLACE_HOLDER_EQN) + + # Comment preprocessing + eqn_bracket_pattern = r"(.*?)\{[^}]*\}(.*?)" + eqn_bracket_sub = r"\1" + r"\2" + xml_str = re.sub(eqn_bracket_pattern, eqn_bracket_sub, xml_str) + + return xml_str + + +def replace_backslash(name): + """Helper method to remove backslashes from variable names + + Parameters + ---------- + name : str + The name of the variable + + Returns + ------- + : str + The name of the variable with backslashes removed + """ + return name.replace("\\\\", "") From ec3732a9bc5c453bede93d85a132800892298535 Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Wed, 14 Feb 2024 16:39:41 -0500 Subject: [PATCH 7/8] Add 2 tests for stella models that require preprocessing to be parsed by PySD and MIRA, test for initials and templates only --- tests/test_system_dynamics.py | 127 ++++++++++++++++++++++++++++++++++ 1 file changed, 127 insertions(+) diff --git a/tests/test_system_dynamics.py b/tests/test_system_dynamics.py index b98ab5095..a158f9152 100644 --- a/tests/test_system_dynamics.py +++ b/tests/test_system_dynamics.py @@ -19,6 +19,10 @@ XMILE_SIR_URL = "https://raw.githubusercontent.com/SDXorg/test-models/master/samples/SIR/SIR.xmile" XMILE_TEA_URL = "https://raw.githubusercontent.com/SDXorg/test-models/master/samples/teacup/teacup.xmile" +XMILE_COVID_URL = "https://exchange.iseesystems.com/model/isee/covid-19-model" +XMILE_RESOURCES_POP_URL = ( + "https://exchange.iseesystems.com/model/isee/resources-and-population" +) HERE = Path(__file__).parent MDL_SIR_PATH = HERE / "SIR.mdl" @@ -213,3 +217,126 @@ def tea_end_to_end_test(model, amr): assert amr_semantics_ode["initials"][0]["target"] == "teacup_temperature" assert float(amr_semantics_ode["initials"][0]["expression"]) == 180.0 + +def test_stella_resources_pop_model(): + tm = template_model_from_stella_model_url(XMILE_RESOURCES_POP_URL) + assert len(tm.initials) == 2 + assert "natural_resources" in tm.initials + assert "population" in tm.initials + + assert len(tm.templates) == 4 + + assert isinstance(tm.templates[0], NaturalProduction) + assert tm.templates[0].outcome.name == "population" + + assert isinstance(tm.templates[1], NaturalDegradation) + assert tm.templates[1].subject.name == "natural_resources" + + assert isinstance(tm.templates[2], NaturalDegradation) + assert tm.templates[2].subject.name == "population" + + assert isinstance(tm.templates[3], NaturalProduction) + assert tm.templates[3].outcome.name == "natural_resources" + + +def test_stella_covid19_model(): + tm = template_model_from_stella_model_url(XMILE_COVID_URL) + + assert len(tm.initials) == 15 + assert "uninfected_at_risk" in tm.initials + assert "infected_not_contagious" in tm.initials + assert "asymptomatic_contagious" in tm.initials + assert "symptomatic_contagious" in tm.initials + assert "symptomatic_not_contagious" in tm.initials + assert "recovered" in tm.initials + assert "tested_infected_not_contagious" in tm.initials + assert "tested_asymptomatic_contagious" in tm.initials + assert "tested_symptomatic_contagious" in tm.initials + assert "tested_symptomatic_not_contagious" in tm.initials + assert "tested_infected" in tm.initials + assert "presumed_infected" in tm.initials + assert "cumulative_deaths" in tm.initials + assert "hospital_bed_days" in tm.initials + assert "cumulative_testing" in tm.initials + + assert len(tm.templates) == 23 + + assert isinstance(tm.templates[0], NaturalProduction) + assert tm.templates[0].outcome.name == "cumulative_testing" + + assert isinstance(tm.templates[1], NaturalProduction) + assert tm.templates[1].outcome.name == "cumulative_deaths" + + assert isinstance(tm.templates[2], NaturalProduction) + assert tm.templates[2].outcome.name == "infected_not_contagious" + + assert isinstance(tm.templates[3], NaturalConversion) + assert tm.templates[3].outcome.name == "symptomatic_contagious" + + assert isinstance(tm.templates[4], NaturalConversion) + assert tm.templates[4].outcome.name == "asymptomatic_contagious" + + assert isinstance(tm.templates[5], NaturalConversion) + assert tm.templates[5].subject.name == "tested_infected_not_contagious" + assert tm.templates[5].outcome.name == "tested_asymptomatic_contagious" + + assert isinstance(tm.templates[6], NaturalDegradation) + assert tm.templates[6].subject.name == "symptomatic_contagious" + + assert isinstance(tm.templates[7], NaturalProduction) + assert tm.templates[7].outcome.name == "hospital_bed_days" + + assert isinstance(tm.templates[8], NaturalConversion) + assert tm.templates[8].subject.name == "symptomatic_contagious" + assert tm.templates[8].outcome.name == "symptomatic_not_contagious" + + assert isinstance(tm.templates[9], NaturalConversion) + assert tm.templates[9].subject.name == "uninfected_at_risk" + assert tm.templates[9].outcome.name == "infected_not_contagious" + + assert isinstance(tm.templates[10], NaturalDegradation) + assert tm.templates[10].subject.name == "symptomatic_not_contagious" + + assert isinstance(tm.templates[11], NaturalProduction) + assert tm.templates[11].outcome.name == "presumed_infected" + + assert isinstance(tm.templates[12], NaturalConversion) + assert tm.templates[12].subject.name == "symptomatic_not_contagious" + assert tm.templates[12].outcome.name == "recovered" + + assert isinstance(tm.templates[13], NaturalConversion) + assert tm.templates[13].subject.name == "tested_asymptomatic_contagious" + assert tm.templates[13].outcome.name == "tested_symptomatic_contagious" + + assert isinstance(tm.templates[14], NaturalDegradation) + assert tm.templates[14].subject.name == "tested_symptomatic_contagious" + + assert isinstance(tm.templates[15], NaturalConversion) + assert tm.templates[15].subject.name == "tested_symptomatic_contagious" + assert tm.templates[15].outcome.name == "tested_symptomatic_not_contagious" + + assert isinstance(tm.templates[16], NaturalDegradation) + assert tm.templates[16].subject.name == "tested_symptomatic_not_contagious" + + assert isinstance(tm.templates[17], NaturalConversion) + assert tm.templates[17].subject.name == "tested_symptomatic_not_contagious" + assert tm.templates[17].outcome.name == "recovered" + + assert isinstance(tm.templates[18], NaturalConversion) + assert tm.templates[18].subject.name == "asymptomatic_contagious" + assert tm.templates[18].outcome.name == "tested_asymptomatic_contagious" + + assert isinstance(tm.templates[19], NaturalProduction) + assert tm.templates[19].outcome.name == "tested_infected" + + assert isinstance(tm.templates[20], NaturalConversion) + assert tm.templates[20].subject.name == "infected_not_contagious" + assert tm.templates[20].outcome.name == "tested_infected_not_contagious" + + assert isinstance(tm.templates[21], NaturalConversion) + assert tm.templates[21].subject.name == "symptomatic_contagious" + assert tm.templates[21].outcome.name == "tested_symptomatic_contagious" + + assert isinstance(tm.templates[22], NaturalConversion) + assert tm.templates[22].subject.name == "symptomatic_not_contagious" + assert tm.templates[22].outcome.name == "tested_symptomatic_not_contagious" From 09a184053109ef59d8aa3bce27825e216b11280e Mon Sep 17 00:00:00 2001 From: nanglo123 Date: Thu, 15 Feb 2024 09:12:13 -0500 Subject: [PATCH 8/8] Use text property of request rather than decoding the content --- mira/sources/system_dynamics/stella.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mira/sources/system_dynamics/stella.py b/mira/sources/system_dynamics/stella.py index 82a4c24f5..e00a55d4f 100644 --- a/mira/sources/system_dynamics/stella.py +++ b/mira/sources/system_dynamics/stella.py @@ -79,7 +79,7 @@ def template_model_from_stella_model_url(url) -> TemplateModel: A MIRA Template Model """ - xml_str = requests.get(url).content.decode("utf-8") + xml_str = requests.get(url).text return template_model_stella_model_string(xml_str)