From 1ed3f2ca086b1a206850b8a1137d6d740dffc69e Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Mon, 28 Oct 2024 12:36:48 -0400 Subject: [PATCH 1/4] Address minor issues identified in 10/24 Code Audit (#1504) * Update code owners * Trying xfail PHE tests only on Windows * Revert change on xfailed tests for PHE * Update .github/CODEOWNERS * Update idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py --- .github/CODEOWNERS | 7 ++++--- .../column_models/tests/test_plate_heat_exchanger.py | 8 +++++--- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/.github/CODEOWNERS b/.github/CODEOWNERS index 64df5fac28..8bd70f352f 100644 --- a/.github/CODEOWNERS +++ b/.github/CODEOWNERS @@ -23,13 +23,14 @@ conftest.py @lbianchi-lbl /idaes/commands/ @dangunter # Core -/idaes/core/ @andrewlee94 @dallan-keylogic @lbianchi-lbl +/idaes/core/ @agarciadiego @dallan-keylogic @lbianchi-lbl /idaes/core/dmf/ @dangunter -/idaes/core/surrogate/ @andrewlee94 @bpaul4 @avdudchenko @rundxdi +/idaes/core/surrogate/ @bpaul4 @avdudchenko @rundxdi /idaes/core/ui/ @dangunter # Models -/idaes/models/ @andrewlee94 @bpaul4 +/idaes/models/ @agarciadiego @bpaul4 +/idaes/models/properties/ @agarciadiego @dallan-keylogic # Apps - each package needs a maintainer /idaes/apps/caprese/ @Robbybp diff --git a/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py b/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py index cd08ffbb67..4b305a1630 100644 --- a/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py +++ b/idaes/models_extra/column_models/tests/test_plate_heat_exchanger.py @@ -14,8 +14,8 @@ Tests for Plate Heat Exchnager unit model. Author: Akula Paul """ - import pytest + from pyomo.environ import ( check_optimal_termination, ConcreteModel, @@ -23,6 +23,8 @@ units as pyunits, value, ) +from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent + from idaes.core import FlowsheetBlock from idaes.models_extra.column_models.plate_heat_exchanger import ( PlateHeatExchanger as PHE, @@ -37,7 +39,6 @@ from idaes.core.util.model_statistics import degrees_of_freedom from idaes.core.util.testing import initialization_tester from idaes.core.solvers import get_solver -from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent # ----------------------------------------------------------------------------- @@ -66,8 +67,9 @@ def test_config(): workaround_for_1294 = pytest.mark.xfail( + # the failures only occur for Windows on GHA with Python <3.12, and Linux with Python 3.12 reason="These tests fail with Pyomo 6.7.0. See IDAES/idaes-pse#1294 for details", - strict=False, # the failures only occur for certain platforms, e.g. Windows on GHA + strict=False, ) From 53f756bcd8a23794f22efe25163ef795cedd5352 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Tue, 29 Oct 2024 15:46:49 -0400 Subject: [PATCH 2/4] Adding error messages for packed config dicts (#1515) --- .../base/generic_property.py | 3 ++- .../base/generic_reaction.py | 7 ++++++ .../base/tests/test_generic_property.py | 25 +++++++++++++++++-- .../base/tests/test_generic_reaction.py | 11 ++++++++ 4 files changed, 43 insertions(+), 3 deletions(-) diff --git a/idaes/models/properties/modular_properties/base/generic_property.py b/idaes/models/properties/modular_properties/base/generic_property.py index 26a31b3d57..2dcb428bea 100644 --- a/idaes/models/properties/modular_properties/base/generic_property.py +++ b/idaes/models/properties/modular_properties/base/generic_property.py @@ -304,7 +304,8 @@ def build(self): # Add Phase objects if self.config.phases is None: raise ConfigurationError( - "{} was not provided with a phases argument.".format(self.name) + f"{self.name} was not provided with a phases argument. " + "Did you forget to unpack the configurations dictionary?" ) # Add a flag indicating whether this is an electrolyte system or not diff --git a/idaes/models/properties/modular_properties/base/generic_reaction.py b/idaes/models/properties/modular_properties/base/generic_reaction.py index 531a7bd06d..6ad5644c3a 100644 --- a/idaes/models/properties/modular_properties/base/generic_reaction.py +++ b/idaes/models/properties/modular_properties/base/generic_reaction.py @@ -199,6 +199,13 @@ def build(self): # and cannot be set until the config block is created by super.build super(ReactionParameterBlock, self).build() + # Check to make sure a property block was assigned + if self.config.property_package is None: + raise ConfigurationError( + f"{self.name} was not assigned a property package. " + "Did you forget to unpack the configuration dictionary?" + ) + # Set base units of measurement self.get_metadata().add_default_units(self.config.base_units) diff --git a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py index 716c824443..39a8e84e13 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_generic_property.py +++ b/idaes/models/properties/modular_properties/base/tests/test_generic_property.py @@ -199,7 +199,7 @@ def test_no_components(self): with pytest.raises( ConfigurationError, - match="params was not provided with a components " "argument.", + match="params was not provided with a components argument.", ): m.params = DummyParameterBlock( phases={ @@ -215,12 +215,33 @@ def test_no_phases(self): with pytest.raises( ConfigurationError, - match="params was not provided with a phases " "argument.", + match="params was not provided with a phases argument. " + "Did you forget to unpack the configurations dictionary?", ): m.params = DummyParameterBlock( components={"a": {}, "b": {}, "c": {}}, base_units=base_units ) + @pytest.mark.unit + def test_packed_dict(self): + m = ConcreteModel() + + dummy_dict = { + "phases": { + "p1": {"equation_of_state": "foo"}, + "p2": {"equation_of_state": "bar"}, + }, + } + + with pytest.raises( + ConfigurationError, + match=re.escape( + "params[phases] was not provided with a phases argument. " + "Did you forget to unpack the configurations dictionary?" + ), + ): + m.params = DummyParameterBlock(dummy_dict) + @pytest.mark.unit def test_invalid_component_in_phase_component_list(self): m = ConcreteModel() diff --git a/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py b/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py index 39f81c1921..722837906d 100644 --- a/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py +++ b/idaes/models/properties/modular_properties/base/tests/test_generic_reaction.py @@ -193,6 +193,17 @@ def test_rate_build_no_stoichiometry(self, m): rate_reactions={"r1": {"heat_of_reaction": "foo", "rate_form": "foo"}}, ) + @pytest.mark.unit + def test_packed_config_dict(self, m): + with pytest.raises( + ConfigurationError, + match=re.escape( + "rxn_params[property_package] was not assigned a property package. " + "Did you forget to unpack the configuration dictionary?" + ), + ): + m.rxn_params = GenericReactionParameterBlock({"property_package": m.params}) + @pytest.mark.unit def test_rate_build_invalid_phase_stoichiometry(self, m): with pytest.raises( From 2e5a9a4245420e16445485f4a04776d0eef3c660 Mon Sep 17 00:00:00 2001 From: avdudchenko <33663878+avdudchenko@users.noreply.github.com> Date: Mon, 4 Nov 2024 10:51:36 -0800 Subject: [PATCH 3/4] Adds ONNX Surrogate support from OMLT (#1308) * fist commit * added tests * undid test issues * added direct load test * adde method tosave model to mirror keras Likely want to refactor this later * fixed acidental delet * Speling fixes and doc updates * clean up imports * Update pytest.ini * fix spelling in utility_minimi... * Fixed onnx dependancy * lint and rst fixes * fix gitignore * pylint and sphynix #2 * fixing rst * doc string fix * updates to comments * linting fixe * added onnx to distro --------- Co-authored-by: Kyle Skolfield Co-authored-by: Keith Beattie Co-authored-by: Ludovico Bianchi Co-authored-by: Andrew Lee --- idaes/core/surrogate/keras_surrogate.py | 112 +------- .../surrogate/omlt_base_surrogate_class.py | 187 +++++++++++++ idaes/core/surrogate/onnx_surrogate.py | 257 ++++++++++++++++++ ...5_tanh_1e-06_4096_tr_15481_Calcite_ST.onnx | Bin 0 -> 101704 bytes ...6_4096_tr_15481_Calcite_ST_idaes_info.json | 1 + .../surrogate/tests/test_onnx_surrogate.py | 173 ++++++++++++ idaes/core/util/utility_minimization.py | 6 +- setup.py | 2 + 8 files changed, 630 insertions(+), 108 deletions(-) create mode 100644 idaes/core/surrogate/omlt_base_surrogate_class.py create mode 100644 idaes/core/surrogate/onnx_surrogate.py create mode 100644 idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST.onnx create mode 100644 idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST_idaes_info.json create mode 100644 idaes/core/surrogate/tests/test_onnx_surrogate.py diff --git a/idaes/core/surrogate/keras_surrogate.py b/idaes/core/surrogate/keras_surrogate.py index 45a3a64a2a..1bc1feed62 100644 --- a/idaes/core/surrogate/keras_surrogate.py +++ b/idaes/core/surrogate/keras_surrogate.py @@ -17,23 +17,21 @@ # pylint: disable=missing-class-docstring # pylint: disable=missing-function-docstring -from enum import Enum import json import os.path -import numpy as np import pandas as pd from pyomo.common.dependencies import attempt_import -from idaes.core.surrogate.base.surrogate_base import SurrogateBase from idaes.core.surrogate.sampling.scaling import OffsetScaler +from idaes.core.surrogate.omlt_base_surrogate_class import OMLTSurrogate + keras, keras_available = attempt_import("tensorflow.keras") omlt, omlt_available = attempt_import("omlt") if omlt_available: - from omlt import OmltBlock, OffsetScaling from omlt.neuralnet.nn_formulation import ( FullSpaceSmoothNNFormulation, ReducedSpaceSmoothNNFormulation, @@ -45,7 +43,7 @@ from omlt.io import load_keras_sequential -class KerasSurrogate(SurrogateBase): +class KerasSurrogate(OMLTSurrogate): def __init__( self, keras_model, @@ -87,52 +85,11 @@ def __init__( input_labels=input_labels, output_labels=output_labels, input_bounds=input_bounds, + input_scaler=input_scaler, + output_scaler=output_scaler, ) - - # make sure we are using the standard scaler - if ( - input_scaler is not None - and not isinstance(input_scaler, OffsetScaler) - or output_scaler is not None - and not isinstance(output_scaler, OffsetScaler) - ): - raise NotImplementedError("KerasSurrogate only supports the OffsetScaler.") - - # check that the input labels match - if input_scaler is not None and input_scaler.expected_columns() != input_labels: - raise ValueError( - "KerasSurrogate created with input_labels that do not match" - " the expected columns in the input_scaler.\n" - "input_labels={}\n" - "input_scaler.expected_columns()={}".format( - input_labels, input_scaler.expected_columns() - ) - ) - - # check that the output labels match - if ( - output_scaler is not None - and output_scaler.expected_columns() != output_labels - ): - raise ValueError( - "KerasSurrogate created with output_labels that do not match" - " the expected columns in the output_scaler.\n" - "output_labels={}\n" - "output_scaler.expected_columns()={}".format( - output_labels, output_scaler.expected_columns() - ) - ) - - self._input_scaler = input_scaler - self._output_scaler = output_scaler self._keras_model = keras_model - class Formulation(Enum): - FULL_SPACE = 1 - REDUCED_SPACE = 2 - RELU_BIGM = 3 - RELU_COMPLEMENTARITY = 4 - def populate_block(self, block, additional_options=None): """ Method to populate a Pyomo Block with the keras model constraints. @@ -149,37 +106,7 @@ def populate_block(self, block, additional_options=None): formulation = additional_options.pop( "formulation", KerasSurrogate.Formulation.FULL_SPACE ) - offset_inputs = np.zeros(self.n_inputs()) - factor_inputs = np.ones(self.n_inputs()) - offset_outputs = np.zeros(self.n_outputs()) - factor_outputs = np.ones(self.n_outputs()) - if self._input_scaler: - offset_inputs = self._input_scaler.offset_series()[ - self.input_labels() - ].to_numpy() - factor_inputs = self._input_scaler.factor_series()[ - self.input_labels() - ].to_numpy() - if self._output_scaler: - offset_outputs = self._output_scaler.offset_series()[ - self.output_labels() - ].to_numpy() - factor_outputs = self._output_scaler.factor_series()[ - self.output_labels() - ].to_numpy() - - # build the OMLT scaler object - omlt_scaling = OffsetScaling( - offset_inputs=offset_inputs, - factor_inputs=factor_inputs, - offset_outputs=offset_outputs, - factor_outputs=factor_outputs, - ) - - # omlt takes *scaled* input bounds as a dictionary with int keys - input_bounds = dict(enumerate(self.input_bounds().values())) - scaled_input_bounds = omlt_scaling.get_scaled_input_expressions(input_bounds) - scaled_input_bounds = {i: tuple(bnd) for i, bnd in scaled_input_bounds.items()} + omlt_scaling, scaled_input_bounds = self.generate_omlt_scaling_objecets() net = load_keras_sequential( self._keras_model, @@ -201,32 +128,7 @@ def populate_block(self, block, additional_options=None): "KerasSurrogate.populate_block. Please pass a valid " "formulation.".format(formulation) ) - block.nn = OmltBlock() - block.nn.build_formulation( - formulation_object, - ) - - # input/output variables need to be constrained to be equal - # auto-created variables that come from OMLT. - input_idx_by_label = {s: i for i, s in enumerate(self._input_labels)} - input_vars_as_dict = block.input_vars_as_dict() - - @block.Constraint(self._input_labels) - def input_surrogate_ties(m, input_label): - return ( - input_vars_as_dict[input_label] - == block.nn.inputs[input_idx_by_label[input_label]] - ) - - output_idx_by_label = {s: i for i, s in enumerate(self._output_labels)} - output_vars_as_dict = block.output_vars_as_dict() - - @block.Constraint(self._output_labels) - def output_surrogate_ties(m, output_label): - return ( - output_vars_as_dict[output_label] - == block.nn.outputs[output_idx_by_label[output_label]] - ) + self.populate_block_with_net(block, formulation_object) def evaluate_surrogate(self, inputs): """ diff --git a/idaes/core/surrogate/omlt_base_surrogate_class.py b/idaes/core/surrogate/omlt_base_surrogate_class.py new file mode 100644 index 0000000000..6ae587d583 --- /dev/null +++ b/idaes/core/surrogate/omlt_base_surrogate_class.py @@ -0,0 +1,187 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Interface for importing ONNX models into IDAES +""" +# TODO: Missing docstrings +# pylint: disable=missing-class-docstring +# pylint: disable=missing-function-docstring + +from enum import Enum +import numpy as np + +from pyomo.common.dependencies import attempt_import + +from idaes.core.surrogate.base.surrogate_base import SurrogateBase +from idaes.core.surrogate.sampling.scaling import OffsetScaler + +keras, keras_available = attempt_import("tensorflow.keras") +omlt, omlt_available = attempt_import("omlt") + +if omlt_available: + from omlt import OmltBlock, OffsetScaling + + +class OMLTSurrogate(SurrogateBase): + def __init__( + self, + input_labels, + output_labels, + input_bounds, + input_scaler=None, + output_scaler=None, + ): + """ + Standard SurrogateObject for surrogates based on Keras models. + Utilizes the OMLT framework for importing Keras models to IDAES. + + Contains methods to both populate a Pyomo Block with constraints + representing the surrogate and to evaluate the surrogate a set of user + provided points. + + This constructor should only be used when first creating the surrogate within IDAES. + Once created, this object can be stored to disk using save_to_folder and loaded + with load_from_folder + + Args: + onnx_model: Onnx model file to be loaded. + input_labels: list of str + The ordered list of labels corresponding to the inputs in the keras model + output_labels: list of str + The ordered list of labels corresponding to the outputs in the keras model + input_bounds: None of dict of tuples + Keys correspond to each of the input labels and values are the tuples of + bounds (lb, ub) + input_scaler: None or OffsetScaler + The scaler to be used for the inputs. If None, then no scaler is used + output_scaler: None of OffsetScaler + The scaler to be used for the outputs. If None, then no scaler is used + """ + super().__init__( + input_labels=input_labels, + output_labels=output_labels, + input_bounds=input_bounds, + ) + + # make sure we are using the standard scaler + if ( + input_scaler is not None + and not isinstance(input_scaler, OffsetScaler) + or output_scaler is not None + and not isinstance(output_scaler, OffsetScaler) + ): + raise NotImplementedError("KerasSurrogate only supports the OffsetScaler.") + + # check that the input labels match + if input_scaler is not None and input_scaler.expected_columns() != input_labels: + raise ValueError( + "KerasSurrogate created with input_labels that do not match" + " the expected columns in the input_scaler.\n" + "input_labels={}\n" + "input_scaler.expected_columns()={}".format( + input_labels, input_scaler.expected_columns() + ) + ) + + # check that the output labels match + if ( + output_scaler is not None + and output_scaler.expected_columns() != output_labels + ): + raise ValueError( + "KerasSurrogate created with output_labels that do not match" + " the expected columns in the output_scaler.\n" + "output_labels={}\n" + "output_scaler.expected_columns()={}".format( + output_labels, output_scaler.expected_columns() + ) + ) + + self._input_scaler = input_scaler + self._output_scaler = output_scaler + + class Formulation(Enum): + FULL_SPACE = 1 + REDUCED_SPACE = 2 + RELU_BIGM = 3 + RELU_COMPLEMENTARITY = 4 + + def generate_omlt_scaling_objecets(self): + offset_inputs = np.zeros(self.n_inputs()) + factor_inputs = np.ones(self.n_inputs()) + offset_outputs = np.zeros(self.n_outputs()) + factor_outputs = np.ones(self.n_outputs()) + if self._input_scaler: + offset_inputs = self._input_scaler.offset_series()[ + self.input_labels() + ].to_numpy() + factor_inputs = self._input_scaler.factor_series()[ + self.input_labels() + ].to_numpy() + if self._output_scaler: + offset_outputs = self._output_scaler.offset_series()[ + self.output_labels() + ].to_numpy() + factor_outputs = self._output_scaler.factor_series()[ + self.output_labels() + ].to_numpy() + + omlt_scaling = OffsetScaling( + offset_inputs=offset_inputs, + factor_inputs=factor_inputs, + offset_outputs=offset_outputs, + factor_outputs=factor_outputs, + ) + + # omlt takes *scaled* input bounds as a dictionary with int keys + input_bounds = dict(enumerate(self.input_bounds().values())) + scaled_input_bounds = omlt_scaling.get_scaled_input_expressions(input_bounds) + scaled_input_bounds = {i: tuple(bnd) for i, bnd in scaled_input_bounds.items()} + return omlt_scaling, scaled_input_bounds + + def populate_block_with_net(self, block, formulation_object): + """ + Method to populate a Pyomo Block with the omlt model constraints and build its formulation. + + Args: + block: Pyomo Block component + The block to be populated with variables and/or constraints. + formulation_object: omlt loaded network formulation + """ + + block.nn = OmltBlock() + block.nn.build_formulation( + formulation_object, + ) + + # input/output variables need to be constrained to be equal + # auto-created variables that come from OMLT. + input_idx_by_label = {s: i for i, s in enumerate(self._input_labels)} + input_vars_as_dict = block.input_vars_as_dict() + + @block.Constraint(self._input_labels) + def input_surrogate_ties(m, input_label): + return ( + input_vars_as_dict[input_label] + == block.nn.inputs[input_idx_by_label[input_label]] + ) + + output_idx_by_label = {s: i for i, s in enumerate(self._output_labels)} + output_vars_as_dict = block.output_vars_as_dict() + + @block.Constraint(self._output_labels) + def output_surrogate_ties(m, output_label): + return ( + output_vars_as_dict[output_label] + == block.nn.outputs[output_idx_by_label[output_label]] + ) diff --git a/idaes/core/surrogate/onnx_surrogate.py b/idaes/core/surrogate/onnx_surrogate.py new file mode 100644 index 0000000000..52b559a89b --- /dev/null +++ b/idaes/core/surrogate/onnx_surrogate.py @@ -0,0 +1,257 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Interface for importing ONNX models into IDAES +""" +# TODO: Missing docstrings +# pylint: disable=missing-class-docstring +# pylint: disable=missing-function-docstring +# TODO: Importing protected _ACTIVATION_OP_TYPES as not exposed in distributed version +# pylint: disable=W0123 +from enum import Enum +import json +import os.path + +from pyomo.common.dependencies import attempt_import + +from idaes.core.surrogate.sampling.scaling import OffsetScaler + +from idaes.core.surrogate.omlt_base_surrogate_class import OMLTSurrogate + +onnx, onnx_available = attempt_import("onnx") +omlt, omlt_available = attempt_import("omlt") + +if omlt_available: + from omlt.neuralnet import ( + FullSpaceSmoothNNFormulation, + ReducedSpaceSmoothNNFormulation, + ReluBigMFormulation, + ReluComplementarityFormulation, + ) + import omlt.io as omltio + + if onnx_available: + from omlt.io import load_onnx_neural_network, write_onnx_model_with_bounds + + +class ONNXSurrogate(OMLTSurrogate): + def __init__( + self, + onnx_model, + input_labels, + output_labels, + input_bounds, + input_scaler=None, + output_scaler=None, + ): + """ + Standard SurrogateObject for surrogates based on ONNX models. + Utilizes the OMLT framework for importing ONNX models to IDAES. + + Contains methods to both populate a Pyomo Block with constraints + representing the surrogate and to evaluate the surrogate a set of user + provided points. + + This constructor should only be used when first creating the surrogate within IDAES. + Once created, this object can be stored to disk using save_to_folder and loaded + with load_from_folder + + Args: + onnx_model: Onnx model file to be loaded. + input_labels: list of str + The ordered list of labels corresponding to the inputs in the onnx model + output_labels: list of str + The ordered list of labels corresponding to the outputs in the onnx model + input_bounds: None of dict of tuples + Keys correspond to each of the input labels and values are the tuples of + bounds (lb, ub) + input_scaler: None or OffsetScaler + The scaler to be used for the inputs. If None, then no scaler is used + output_scaler: None of OffsetScaler + The scaler to be used for the outputs. If None, then no scaler is used + """ + super().__init__( + input_labels=input_labels, + output_labels=output_labels, + input_bounds=input_bounds, + input_scaler=input_scaler, + output_scaler=output_scaler, + ) + + self._onnx_model = onnx_model + + class Formulation(Enum): + FULL_SPACE = 1 + REDUCED_SPACE = 2 + RELU_BIGM = 3 + RELU_COMPLEMENTARITY = 4 + + def populate_block(self, block, additional_options=None): + """ + Method to populate a Pyomo Block with the onnx model constraints. + + Args: + block: Pyomo Block component + The block to be populated with variables and/or constraints. + additional_options: dict or None + If not None, then should be a dict with the following keys; + 'formulation': ONNXSurrogate.Formulation + The formulation to use with OMLT. Possible values are FULL_SPACE, + REDUCED_SPACE, RELU_BIGM, or RELU_COMPLEMENTARITY (default is FULL_SPACE) + """ + formulation = additional_options.pop( + "formulation", ONNXSurrogate.Formulation.REDUCED_SPACE + ) + omlt_scaling, scaled_input_bounds = self.generate_omlt_scaling_objecets() + + # omlt takes *scaled* input bounds as a dictionary with int keys + input_bounds = dict(enumerate(self.input_bounds().values())) + scaled_input_bounds = omlt_scaling.get_scaled_input_expressions(input_bounds) + scaled_input_bounds = {i: tuple(bnd) for i, bnd in scaled_input_bounds.items()} + + # TODO: remove this once new OMLT 1.2 is made available and includes tanh support + # overrides default available activation functions for ONNX, tanh is not listed in 1.1 but is supported + + omltio.onnx_parser._ACTIVATION_OP_TYPES = [ # pylint: disable=protected-access + "Relu", + "Sigmoid", + "LogSoftmax", + "Tanh", + ] + + net = load_onnx_neural_network( + self._onnx_model, + scaling_object=omlt_scaling, + input_bounds=scaled_input_bounds, + ) + + if formulation == ONNXSurrogate.Formulation.FULL_SPACE: + formulation_object = FullSpaceSmoothNNFormulation(net) + elif formulation == ONNXSurrogate.Formulation.REDUCED_SPACE: + formulation_object = ReducedSpaceSmoothNNFormulation(net) + elif formulation == ONNXSurrogate.Formulation.RELU_BIGM: + formulation_object = ReluBigMFormulation(net) + elif formulation == ONNXSurrogate.Formulation.RELU_COMPLEMENTARITY: + formulation_object = ReluComplementarityFormulation(net) + else: + raise ValueError( + 'An unrecognized formulation "{}" was passed to ' + "ONNXSurrogate.populate_block. Please pass a valid " + "formulation.".format(formulation) + ) + self.populate_block_with_net(block, formulation_object) + + def evaluate_surrogate(self, inputs): + """ + Method to evaluate ONNX model at a set of input values. + + Args: + inputs: numpy array of input values. First dimension of array + must match the number of input variables. + + Returns: + outputs: numpy array of values for all outputs evaluated at input + points. + """ + raise NotImplementedError + + def save_to_folder(self, save_location, save_name): + """ + Save the surrogate object to disk by providing the location to store the + model in and its name, as well as additional IDAES metadata + + Args: + save_location: str + The name of the folder to contain the ONNX model and additional + IDAES metadata + save_name: str + The name for the model + """ + + write_onnx_model_with_bounds( + os.path.join(save_location, "{}.onnx".format(save_name)), + onnx_model=self._onnx_model, + input_bounds=None, + ) + info = dict() + info["input_scaler"] = None + if self._input_scaler is not None: + info["input_scaler"] = self._input_scaler.to_dict() + info["output_scaler"] = None + if self._output_scaler is not None: + info["output_scaler"] = self._output_scaler.to_dict() + + # serialize information from the base class + info["input_labels"] = self.input_labels() + info["output_labels"] = self.output_labels() + info["input_bounds"] = self.input_bounds() + + with open( + os.path.join(save_location, "{}_idaes_info.json".format(save_name)), "w" + ) as fd: + json.dump(info, fd) + + @classmethod + def load_onnx_model(cls, onnx_model_location, model_name): + """ + Load the surrogate object from disk by providing the name of the + folder holding the onnx model and its name, including accompanying json file that includes following + structure: + + "input_scaler":{ + "expected_columns":[list of input_keys], + "offset":{"input_key":offset_value,etc.}, + "factor":{"input_key":factor_value (e.g. multiplier),etc.}} + + "output_scaler":{ + "expected_columns":[list of output_keys], + "offset":{"output_key":offset_value,etc.}, + "factor":{"output_key":factor_value (e.g. multiplier),etc.}} + + "input_bounds":{"input_key":[low_bound,high_bound],etc.} + "input_labels":[list of input_keys] + "output_labels":[list of output_keys] + + Args: + folder_name: str + The name of the folder containing the onnx model and additional + IDAES metadata + model_name: str + The name of the model to load in the folder + + Returns: an instance of ONNXSurrogate + """ + onnx_model = onnx.load( + os.path.join(onnx_model_location, "{}.onnx".format(model_name)) + ) + with open( + os.path.join(onnx_model_location, "{}_idaes_info.json".format(model_name)) + ) as fd: + info = json.load(fd) + + input_scaler = None + if info["input_scaler"] is not None: + input_scaler = OffsetScaler.from_dict(info["input_scaler"]) + + output_scaler = None + if info["output_scaler"] is not None: + output_scaler = OffsetScaler.from_dict(info["output_scaler"]) + + return ONNXSurrogate( + onnx_model=onnx_model, + input_labels=info["input_labels"], + output_labels=info["output_labels"], + input_bounds=info["input_bounds"], + input_scaler=input_scaler, + output_scaler=output_scaler, + ) diff --git a/idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST.onnx b/idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST.onnx new file mode 100644 index 0000000000000000000000000000000000000000..b0101b3ec01e2ce327467dbab5a3c3cd05ac4bf1 GIT binary patch literal 101704 zcmeFac|2F!+y8HBAQ2j*NTg^m7jN&quC0`!QW{hmBtvF0B#k0dgAhuJkd#C!8Qyzc zTN;!!r#Vtd(LA8}r*l4^`#$IVKKJK7pYQjabKl?le*gEj_j*`suf5i_uh)8B>ssBr zDR&D93-%9m@lfkxtZ%4qXi>PctK7E%a$P+A0z!i2di!|#xj6?;^>=mi>AW54=DE}( zSWfxZSC@D?2Px|s{Nmy-f9o3>Omy@0b@C4h?);BVhH8U;dH0vU|5j2<=12MAigI0? zeF8k3`$$NHvwI{Zhbzd*EO86|_rH6}bqx-5_6r)PD!E5e%J`g`L6C=YpqoLEzk6`# zj~^jHZh;2=0d9T<0b!&5_Lv4gJ*87Xn2DgX>Ls0bf<65GO!Nf?L7u^GqXL{=yq%Z2 z1sVL*CkB3g2EP8TAwHcKettUrfG}MnBU3#KO?Usm&PrV6{?tl;uKcG;P0h_aYjq0p z@eg(StLq+ZX!N6Q7iS+ICr{r1pFj7ZUHtumojZqJ;E(4vHZuEP-`O7izCRDDL?ChL6 z|C6uNQn${ob`JLR_v;+)P9edbK0n>-oLYZ>2S!Gt|JT2RKisr` zIzL{fxxnOyx%}znIyw8fc22S2pR0Cq5A^qS>a^-^e})bHbt?qsf*-AL33Tf;0ZyGg z^A|n6k&)3)kN?wCItBaxRn?9*{pn@9tqAZB4F1#K@o2*_oxS__%Kmx_{gZkA*;@ZJ z}E&9oU3k>o5xtUHr z{{8{~_5Ec0v$1mW^y{3ePJdy0{{Ff={g(b1N~3=s*lsStA)R)k( zY;UK>=IU%>;O4hH$lxarB-lA%l&jnF-+0)a#au!H1Ks?B4ctS5Ljr%e*&dzd$@Qmd zf6U&_ck<`fV*Zn93iS{4{_&apM+5O+c%Gjd|F>=v{N$+q5fAX|e`{{o*)L8@{QZN2 zI$uOU=kVxs*Z$z9oBWue|6Yfj+{poQzgy2K5Ba5rGjMhZ_FV4VDW42xIs1A1yXqV8 zPfBUY{3vUD;m=iYe`{H%qx)mLnT-DNdHk!Dlvu_6ua(UI$k>Qg+#kM9^NS^uzg}<{ ziB;Smu<9mOasRF*#VYQPUD5ruMi#5M-@67VR&l>?tyHYy#47G5ob=y7#qAy__q(*& zuPLALKcsxd`u_>#GydI_@8^i zF~4&%`VHq(|EJ}|r0l^bAS>KrdjP;HA-K_8Dl48~;W_@DTCyoaFpPjDu zH=0eH_vMesw)~q;74Pd4XYc&OGI{@+=Zf>bi1WVuu{pBEd0)hNU;ew~efg)X?>Bnx zKV*F-zmxU-Tu#jT#H>%u`u^?I1@SoyE*}5KC$0RKCKC}WIkA%azdbeR*IGlI!Sj#( zH2%7G#Y#@B-&v9`!(w`|A(y4{6AxT=D(Zu{ajMa z`oyeH%=*NvPt5xMxFbh@3&F(u`o#PC{v-+}-q$DI*C#%y?_bYHBW8VK*7yGe>-&wK z`wv;4>F;EHKbI4;zCU(B^B1W!|626%i;BfLJjFRY!*q@Qd+n|K!;gXZjkXfAJ~8VP zvp(bh^Teq};xpI9XRiOga}UHhJblGEJpb^6%l{{xxgP!>=kWZEK9lR?BPRS{vp8uA;D^zzdD!U*Cl_#eqxv>hIwL`Cx&@qm?wsL|6=g& zZ*;YI<1Pq%;gTIGyJ#?;!))@lpUAhOi+uG9J?hK!(977L0 z=IEq%Q<#m+xVBq4c%6*~lbLdCnA2-+M^`7Bc)wKibiWMld#A_m z-xN%~MatpMVb0ibvyhq%NP)4;N#s_=62*=r=Br)|x7_(*#V3OsBz1BlG1OU0bb}Yt zsjd@fcMmsMasQcRfw{g!>@GI8`VSGnH2yUWcway{W@}S*DvCPrskZ zgG}$?Y}?aR_F_yv#xF@Cd(>Gt<(_V_|J6i%zjFw$X8_Cttwo@-#}rq_Y2yX&Gn{Fe zEFN9586WP^NfFH^f2U#ex2##4P!v-dU9;y4vmtD8yY0~3^6 zTg`>_>6Y$MMD9(Fk?HHuzm`OSJ8%p`%IaNf+~_dek~U_ zErs5{Io+(<*Ktf#r9LWm^noXfytzk9s?_Fp2fb|8kU1|oG~xll)Z7mkb1~D5Vtoym_F|of~oola`mMy zDUq9ok8ZxgVjX`h^<2mMDcy%R8JhgZ|Q z@EWy=>7ZtFe8`Jr`wHO~c~Q&Lq7zeI29|nrK4ycyj|G$1q4?&01ha4ZO=hrjfhc-) zEM1b?mE1b=A^x^MH!G+NryY&r@_N04p%qbh zUOPglB&opcyQj)Np1Y64UaX-4Uw7&(rG){0jYPRYl9ylYjS4gnmm5z+w-4(X57T7s z!?AMKWA7ac%|%zqmYjS($F@J;I&C*>*;~nO%T$1UgVawXXbj>pH;agDJ_$#-XJIbs z0$X<&5Z^cxHhI)KM9*}Lb>7B1jnK4I<0UIQR_x<@1ieCL$V_4W=2U)i`7Y5Ir6{iT za!*j*M5$SnKJ!6tl`uwnH&+-i0QV7XX3GB2=>L|c9?r`{5BI&K7g^^T4JFjF@atlR3v=YR=~79d`Vg=#*x&p=o5tnyP@mW zX?W*xB=nkE%8W{sC(-%)EqY0ACxs{8z@w|jXv*`epgHF?8T&i`@Jij?9YdLvG~thvl)!qP%v1)fH ze7A_%|G0no;k9LS&H7lfwP-GeZaK@*zPrbsEm9=omnCy$r>An8%2!oL9GABocOf53 z?UhJ;%MNxZ@5`<{6c0^@4A|&1e$?G5ojo@+6Ta#7fOE}N<%PK~*^}dU;JNEUQs$^k zkIBBLja`!AP>?O=ZY7mB3`^*oUYd-n;|@`4S8bu{NDWrK-HK7%9gD|GvT@rzZB9vV zGj$xv@!B3{wDR11Qra*KeK(CI_mA-v8Xj7x-F5{c7hR;!e2V#TRxIB)tOu%$?9b1!lZZz=v;RXl8fMih;H{l&xI`hRyTo=w&?0#4bf)vKMaGNTXd&%*5*%MiBX+ zFE_DL5dsD#3jGV4VdvTs`sKb0lm0=5k%%(n_AkpHYx2;Z&a!?_xJ*g#*^tc`o!rKJ zuI_E|9+RMAt|ApIm*;yd-bNmX^r&Z)9oaiXn;0fv;rb80M-Syaqn`^ZiCT%VXiDA( zQgpnVnR9bmX-j+yv9&r!&3b4s?3oD0ebF@0s4yKm`fL{CsFx)CGJIP3fH(HXmo=Pc znmo3PE>RgK;I$^X@~sBiq%>*RV!5%pD_78GPFl>UFXi0SuTL!WgVxa%b7Muulj69W ztFBOKx)iD`B`lqGzT%qJ>OfQ20isr;!@WIyue`csR4?PCn&*8ijXfro&46F6{8s17xe?X)4=VNnY%kO4csk z16zj<#$nd#IC-2MiCDjwZoqhvU5+|iHr5`;$fYq_{!h80hP{mX#nw!(07r79%F^jyGa(VzSpLeSBt_z_t%qNzRruGF!FA!*U~OVfVM+@ojDn z^}Ch`Gkg|tS^efiptCJj-^u6Qc5K6-F$(;&SusRINr?RT7$#!F9VX1FiG)@kA+M~r z5YNO6M%BF=4N9wl_Iz!8tLDJAcKtxdDQt%&<35`|)$R#G3HLJ1mhWZ`Nrt4!yC3!6 z>LuLNgA*p~++C5lAdFj3A;sKZcEuvmMMgL#hA}T_D#4AyMda$eNVKf-Mf0{;QcRX3R- zsi8P!b0MbCab!nS2wq=p1KrJn>3ENNT&6di^@>lzZ&TGd>j;+0FRB%l8pdNwt0Rsm zTEXQGXF>C=A5NTn?8F4MMRbM!_l7(L?c#CwcN1&nY&zDtC=lbY|;w7`OFmhoI z*Ksn1*5qEL4W+ZVblagcv33TRe{YA-w|E}AePukpa2Q1AUfu|gZZ64%LW%8{7g2+rO??Gz1gL!E_0RpH-Yai z4R-z9CP-CWh!&f=vlkp?%Z6?424fu3h_h%C#7@~qs};`Sx_HNuaq=f<-HKU6(r-7` z2Tz5(L-`d3LnWa4@^Q3pGi5lZi&($f7o$#}=FW`F6e*d;kWp`HpogU@T)izrj`Vm) zFIC*?M!%c5W_oKSm$>C^9{Qp-?SadZz9WQLQumybjA+QIC?Z_`NlSzW#$A`yo-1a#_+ zAN}zwy;eFKFeLE2PUd zI}7o_HozvG(>=5VFrY^*|3g>%tD^Ggn2qM=E@B;Vz<} zv!P=5%5YFy;e##nrc<*%%_QypQa*9mYi9tyd!9f%UFVpDu+yooLRR2{+<|xu$HMSpjF)t5vL-F$it17Qw3H z8^JwNnx`)((f6y&(Tzns-P#R4W{v^V$0NYx-c1aX4WJd3(JZ`whoVcZaB@f&oLEsp z7AE#&UyRrc6E61P6GlFUUdCP)r}ljy-S2g!%GSF4g^`x9z|Y*WXpuohl;=PioH31m zVx7i?4*AFlXKW^tJDhNtR6b5WRcSGI&Jw)-tVtB0m(EEyWz!lanlumE#drPo!DPg^x*C=e!KOx1Kkc zU_{q@yk>l7-bF?56R74jokAC0sK4rX4Kg|8D|GQ#wQYA1#-hUVxPB%e!nWq zJ4
qRY$#lxEN7 zhHG1&;CnA+{#jZo8B3m#swMJ`Jo**Z`@B3e1-6!XLr6fFcgQM z?a8j1{uoayE)(RHNQ2ahm2jLU!|=JO5R`Pca!`K1ipKMqBzSc_>>jn9_^k+}Z)|Bc?Z6T8{H=gpJX{V-C)nWE0aqE7${JGE zXwTX%Fd(l^8xi*8dR99BF0co3xd)Z~*whh4CY;Ve;X+F(WUfyGsRlhLZ+Rs;CTV1@ z_&}CR@7onJUYIfy3#*xkYZ-KY$w27Se<@!(N)Im12_=Q8%Ea_~ICpW!5n?fSy-=#9 z7kIQQk=popTG#gwt#`{5+J$9uZ{{80Y-IbH3v)-}(TNcltJO{i_Sfbo^&d|QMi0ZX zu^CV)y9y8Pmx8GqB3`!+~g4d)r&8 zbtIloJW_>u##=zeUqiTOy$8H{8bG&i*aRo+3&^9nR;>Jn?!bKQ3LD;>#AnA(n!o+- zhM#>6Kp;>iJ=Qm2Y2sPPli5)EVE%i~FSUrQd=gE*Zipa_A)kqfzcuK8iNXyxqOoAO z3i~#z6b2tw$AwD|klHKv!MO{?_k&}h`(UKTU*8kSU22f2I|J9H7D3w4+hwKKE!Zv6 zYw*^gBIuB+#*tkz@Uh(}w45;%cICgrT(1SNw6QyTHhr0(dFB-UuEAT7x$un~XbhlP zMMkK+M~X~ey@E;8_Jfq;8W4BRp8fE;3-u1)h`rvl;^?P}EW~BQiClHASCj?~Q8!{Y zJ1Rk|-bP|u{0s-DUlug%aE66*jo8Bz>H$Tv5L^D53?3_4VP|!YJeau>7RB78i+b** zW~GTt$meAA>N_6_JBEO%c5ikLPUaIsY-o4-*ka?_p&0u05=33h0o94Gpy}y){=pyv zHsf3nvoxiZrnM}Z+sVK!Or9K*!&n`ZJLgu^YV#f_Y%6T zxGNe-y#v>ldg6I)Ja(THkIm<{gQw~yycA-GFHcV~@2Hl53F`Y$>h3BWbyopGZ;gk9 zBNK?4sX9#Ex&zZnreJeTJRXYA8Y8#g4g|*UM31&?gc&0jg3^sZyq>HNeWU8YFG-6( z+>(S7x2Yo6a2hQ*JR`N$=49SFRg71A%3ZX3O0qV^Vzl5IRViHz!Vq;hlx$G( z{m==r4^*ht%X5`PGnPItI*<1v(=1vAU8(D~0J68obQq>`o$z zmT|`5C>cX#_h(U^o1eI?>P43Rbxt&VY6hM!>EMQ14#kQ98S1KCC$wpkr8XUV%a>Kp zWAzGGbJ8hIXm@N2tdtF*lVt{>>}hYDD!K%RjaNbA_yRb%O=QvfULNBjFOlYB<}l2K zfwu2mXsk-TXu?8cyr`j%!S=l?KA$pX){T=OhECB`@W_q5a>p2=XN$0r|46?tvt_$+ zwj9lpBG3BW!Eu%*#P(o#x!uFA{Edcme86lY&TDl@(UeUz|133sAM*lc+d5(1Q8}82 zB0el8lpgChkgjagX0MNz0qfM;?Bdh|{5xeHZQFOk)b01l3B@AzlbSOcPWA!G0Y^n& zyQ}e9-7XU4*)!m1j+|iX^Ic@X8ICSlJe{3UIE0UIT?p9Mflc?1;gtr3!|8dCVAe=& zc6HZw^Sh~8eA3P>{5%a~rX?8hft3>9MNt{*va-qQnHgBuuNy3^p2BaoS%?ExeI^68 z_TfffccT^yq*zyn(P+Bd4Hor1$x7_8XRq#QrcdK-*y+=^nfpKQi*HxR;=^Y}FrX}z znRap>ZPg!ae$vVwp9M&<-ILwOQY|AocSZtDKfi;z*p1?i-uA=3Uc>nH*Q=<|S(|)1 zz76qF488ZY8Tyu;0WHV-#4qtNCSMLD^E0R6rwnPDzC>8obMZ4u&fg)DrmMl|^myKl z7Q(~#&Gb0im))-1M7PT)ll!SbcyMKaz=zs1J4IKvxd-ulel zr3=hIxsSzj@>#en-VF~vBP}Fq`v&fvkJ6NhzId&+X?SSs?Y>qbHU0+FS6#Q zBKGHovSW1CaMgzrsaE=2ZteULY|lnh(p&NjO^e>lORmqRp68=jPTGO+d*8!~1yuCu z9ml<>Dd%zXXWVndg|4|q;EYB<`JAja+R*O=X^me5oA5c=-q;28y~5z=lM1wK3c<{T z6dX47GklPrWid`NpKiOTMjjLy!u&ZBEKbNJ9urIG^}~_SX51lMJN3LUrFsvX9cqQ= zc23~deanQVfwr{uU>cR0nus3hnlM3s1_m8c0nzB2m~{9hGrR4&`L)q2@q_sp2D_fH zJe+Gm>W1`S?@qM0FjTk4OxN>V(Za{f{nOt?hZK)79o_ff4Bcnc>f%M}W;dTKaq&H`+bBwGFJVTm4cyN7UCj#%X zlWJPtr_~qy@cYo`L{)VStQcrZ>W{962T5v{dNRl0klqo#$2=+2D9B;6KTZ*v_L0Tp zp=)uYXcIftM#>^lrU$>WHjE$PIF3YR*mIL3j&h+|y(|2^uUINVF3VUSA;Gbcl>^)+ z5tW|(Vc6Fx?BFN&n0J;QmSXV}G;)GW`KZR3?vH&ZdrJ@zVt2IoO_B7g^Gt>+6IEo-3h2=X`Gw8shVh+o(?T zZg@EM2yD^U6ZzX+q#f@fL3*D8dD+7O-3z74lg)Rd|Djb_|4or3?<;{Pn+>UBR}Xqa zP7@aC&*5Tax5CVM^U0ByQi9LDmSL#VLOOSd8kjBV4H^rTaIW77Xj>ge)@feDtghO4 zt>!8!r;o&SLwAsMmLuTgx~}{cT*JO8GQ(j`4!F)smX$D8Vs|nUtl8W|kY(PMe`!C1 zzJ~MZ0j&=>=*3pf@@8Lt#`#`?UZ=CL=D{%hAis|2n_TMjkp$#j?p7LkC=yj(t>jK& zZ{9j{1K!p6Kq%|w&U zJtV~*miHxVwr|EtLo@b6W;MAmIg{@rbjLj__u}Ta*69573Es|05xLEh5r$0Ig(i%E z+oSi0=5|bEL&t4KxA0LGuX{eGYrJm|H3L)2Z?pST2N#b2-u*#&ziJ9Va`lNq~FX5jrhj8H{ER=D!8RlW$$)*wM2BvHau^ z(Zj0yczRzkv8a0w!EnsHS3wobC*`Ec;Uc5mK;ri{xfsQi9ON$y$(la z>*5n^YRZiAr~TK#9P)^f(eWwI4xU4%N<`8oo7=SS3k~wNsDh4pd?fc)t^+6W&Si2RC=FNq-5a?WLKp zeq9+@9gRTu2kmrM$9LwWlsEWZsAl%+Sqa?dxD!471nN|12H#(%WAGbgR@EW}mJh$i zzqt@jjxLF#CkG8fOJORW^Lqx5(`A@VEbUtUoV?n%8(GU5@HLS^uP!p|@gg@mP+JoWM;(WbXTYyq z@Cf5XeMpWI($y)}_+`};OwWE>S~JBL21V?H1AUvoT>2EPZMub#bv*^Ko|aU0coGE3 zM3CGZYp!^S1haFgF8{<*274@*LpPfNyuameewxi;82D`xj?Ed$f7VT8maQ*A3zuS$ z%9+G_J(|Y$C|V5bYChrnrNL!d4(YHU#0$Th2C%_jZRoCw0M6`tFFNs5CaBeJq~q5P zgm&LgMXaSKJi^Je)~&UKd&FZCsu)+j5MJ>FFvC_XUyx$ z?!k&J>I~z0n$+}{2j?YyNx%UM=zWC(e^eKzJ5f5Tx`hs2IZvRy(;7A%oI*WMo+S_0 z3UL~h!|4(mN~?Muq3xF}{{Xwq8#~h+NeH6RugKX|#2XHzpSDswz24kEL;gLs!VZufg ze6w-_xSY->mkv0mu4a{$BQ33S{_S`vS8Hf zN(`MmM6k5q8#-s*CX825qDrDC^vQ#(Q2TOInTujAL@4&5GYzeUBU45(r?<2T^%6Fa z)D0nU;9ytmzD5#s-cNxAMWM`%gpXuE<^s_GEp0k3l&9k&roq<(RnXF}3Y+Z5kjvXu z@vzHTv-SzLWK;Z9h*dU#mxZ6mz&&Z?se3iEEbbE>v9~`C`h1CcHh!mPlgr9_Z??xN z_R-|sgCTG^pnyc>*E5RpmuN4!(|D48q@y;qRWy&?!M@0NLKnTKM$UE*@^Wg zpUJI)tR&N99!2^} zN^)DaxRE)fPgv(LHB7dR1OG1#C9#gXJ+l`a zeV_@IZkfz`917jPd&1}M?{QOJZ!+^rjU){iC5@@iwd(c}FO()T>h)_`#Ua*JS@!4ZMD-luDS((JyzkVEp7& zSgsO4>#iAzpid>(x_L12_wy{jR4;+&tDRy0_jZf?FPC7&A|vn%nS;BJx(f@uLTCy* zg!ih~0r`D7WQ~gj=JwS_X^Rk&y2BHsj*WPreUIDc`l|-C|^x?L= zp2>|n6j6S7^L#$?h9#LOR3$?iuag$-zR;afLCZP)aQ@SDG%dZz9Mxzg31iiSvQxG) z>I=`2blHB8bj*z;`{$9ya8){I{A95Fx`li8S_VIwhe2!YI2w~tNvvi1^L-obuwaG{ zss@~-UG7Xq>+?gonY|`~@>Dsz>MDV%8v`tAwkcuLgR!hCjl-q(fyf@84|>vt<;tH^ zxPE80A^qM36E4r=ROVZeqFMpp$NCdB$lk@ZJ`ct31A|cW=>g_V4>|tG{RV-$gE@X3 zf3?#`ZY0H{_Tl#e3AVWR4Q8%uGG=P_s<`cXgzF|spetQu@cl|%RC0>KJ6lu0RU#Na zDF?z9YaYiOQ{#{9m9Su5>BBeo!F1#819bR08)z9M#rM!uB2F)=u(bRtk|U2vPHGg) zHm+a;E*$1=Up$GkKWl@m6b~DxgrV6;X^86?jVh&MiN}03wsO!|#%G6g`R4v{;Bej? zY-%<$3pFl+&aENrMeCh{^NMl6jd{$CHA*7eB^D5$QE#Y=yd~?dd>+i26S>7+XHag2 z1Exqy5}(Ae6zq<};=(d$(W=3Ok}8OpB*oh=ab!oISp;7duj9e<502X&7Q)dqTRgq( zAkIvQ0;$MY&{`jXR;#+L3ZO_CHb0zpw zYG*)^F~tRA-D%h!StL^zp~cf#e629KbItQAZjn1c2ELuimPF{IYWZ#=TH2RwJ5&c| z10o>3P9Ge$D2T|Rbzr;XBfWEEU**ic+2oSC8s%Oeve=!{i}jLq2aOZY=%P9U8hUgp zU6Vh7I-Ze2+mc6AWLJt8yT%bSt?j~qBOAc3`lMx6dpYO0zyxBG_36>aFG;9t7#ZZ! zldhEeOcEDKl*xaqWS?KxfQHvwVf2;WXnA%MxVj54p-vfY%vRUn@Z{?v7&*I{eA|>wj(wjf2%35uQjCtmCdC-^9e$jC3X^9YSI>hfi7!yA zSRcZD)4^XY9$Kotl`s7`o2HJN1+*y~4{}=cdVe|4>;I6HgjLajqnFa?2N#Krq&YE| zCt?&*lBiC5Grg48#9VmM3JbUR61g1$(XQ~p*sWhX8g{M;7Cl^!*(&|Z3-JjEYxZNV z_F(p8N;GjO*1>_Aiv$Kk1L5(YgJkuPktF!yMkY&P3$r6n4nChzV0Bee7~RAa^f6$m zS#Sbp}mvNt-c$kXJO)XvghYLJP-!LtHbI^5c6k)AiF|n6wNqvSJ zdX_fgngJPHp&gI!3qDoM>beAmz17ATW(0q3iXF3{_j+jTYlSNxWpQ$ydM?7V+G3#l z4bHUXmid_8RkYRD2{dM?Q{}DNyv@X1cxc8r|2KfP7pPcW>MYI$2==tg?+SQ<2%j zXhvOtXKgW@^vX1*&@7e)&DJM9du~GSpfUL9-U^y-nF^}FlAP+d`1IN*dZb?pv3tH1 zkLJ5Vg|-e{n^X-VW;#YbY=GdM%gT4(G3bm_y+rP|R>J+bK$Q1;fzo;!{HK_4WVm!6 z$UUKtjnZ~}wIHW_rd1xC)JdUlR|`8o5o0k{&k?SFk6_ZwlR$c%77ag}Lc-fJA#r*( zc)9Fk%xfySEcqQorl$e;cQ+$<59i{J(9XE5t5vv1*^oMX6Op59tIO4Hgk#slsl++U zqtdN2=CeYzh;wSp6`gs^%eTPI| z8_2ZpZzPFVB`mwO#nJT-kAp_xIhy^goOwoa@XN%9%-|~<$g!`|<@uV4;Qnzgc{qDJ zJ#8|EpLpX8w%*(g>O4bZ>ykOU>t>L~59<%-3}ncLxYlmQl=M#g9FIS8HCmj)>; zC4MDw;1e6g70)_GJWe0t>MuLOE62-tT?S0-`U`Ril?MSqZa%-?ZpDg zQ*dXv6Kp(mp`7+;Vs3>fLh;&T;Cb#7OdxT*UihWb6)U5mr)^hIHArF9)blNx*RLZR zzw~E|wnp#=w`7torTOSlDMYjBEy8IP4`{NJA#s0I4)Y2WNO0eQkfWo8Q@VG>aDh5R zWXRB$7nKA%t+tS(i7)AzCohB&?~7?@of8@0Jp}7X7O~j443DUp(CLoz@bc2$6*Kn_ z27jwGI<4dhAMx@j>Atuh*;Ew{L5w}~?S41QofonxSE$zMtQ+&{Zdt@I?f@JNH)ID6 z?86VSl!I+&yj^&x+!<1WdK24^iFjqe=F*DeBFN4kL5a!)I;hftJ$Gs^nQ-$cKD6%7 zhgBs%8YZ_cfIKb>FQoP}B|J!xG` z4NSNGiaX=BaeF&|zX})(_cxw`%JrY9=cx~*@8@lB=C}zgJ2aSIxVxOJ*u7ft&3rzM zzL7;HO}krJ6dBItT7D!Rsb=PbLv_i6b@#apDHSeaK(i?0{5;Z*%lSK2duhu3W^Uoj zBV2Djml;<#i%GTmP7kfwZCUKv2fa7%q6bCsRMAd{NJ_pH&32eV_b(X+5AIZycO*<< zUmY$aGj=OO&4VpidOMg)8ZF5^HCjN%_8rJv%p1+ML`Tv$?kg<%nCFpLBQ36{0*eRl z9^%|5R6@U3yJ`9&X}II80uPnK$QZX~48E&PBO*7V_FYe2)>#SmInUa3tY@6CUzT^Abq0vON9#)mTfQ~lPqZ)OVwT8RkZGCCxP6-nukO`~ zBzL=TZnYsaru;l>+sP*U`=^luBa6v!eX=Aw^C*Y{c82^2)Bh91WH-nYj6Nce>%(9(45l zh;PbJP!YNhN=xHm-iL6?e=6ddW4cjkof=x`aDhzkp@?nXh0L?!jqudKl*6$~uqo3O zHg8;miAPt{2E2?lRtFdZ_3K>O$qrI(x|1oD%PBuva8-B{B4M2J0dRPb1}ZZy!SKLz zCbF>_Myhd)$%k{)I_)&4y?&uZrT%0b?vVuXL3YBq*~yg4oJD57T)_=Bl)*Qbo#3|o zSV+HT${KD=#xd#)9-fx~O>Y$0fT&3RW9J%FV{kiFeO1VB$eG5kl(8Z4kCySnS1DHX zi1LPEThicQlMOCf(+Bz*y@dU*iUk`s^~YvSV^SUwRk5UO5qG(LCw;FzjBa&|gk0;s zH0yvAyFM@p75dqdh(mhxTYnalPE=8eVS2o3SWo!gFo-I3`dlM*qDkbGV`N3}jItGF zgCHVF4>sPvgqt*`v+imdSRa%C<{t*}&yt_g=BOe%+1U*wDr*Ee$CUUXC!WJ0hfVA@ ztwM+oPRC<0)u8I#lXp=*&iiz;hp^_ZmKyIA*m+4ln7%rmVv-l%eBuraI=_}YyjO-c zv%Dbq_#DXQBZp^Cc0mfL2aH)>e97S{!5_EP5A{vz&H z^G)u}G;@xbH-j5jT*W^ec!Qj?=}T|sv$&^y1YO~=n^8}%A-1a~!&%XGkyniyj_q7y zFWhU&SJiicdbdJ&lUQrEI&v63u%P(admO$Qx(CN;C(zneN>IDz5J;)?LyO5fut@DW zvFJj}F71wnyveso`&0eO7pethfp;E_zSjyzhH3%mm$3sB-_xTV>v`S4PJLmMh2yNB za6=WG=xHAxno~aq*G;fRX`A&T663>pI(#6O)5qh>k_Lg>I1{|0d5)=hy^DtJFF?g1 zyUCTIdgTkcHjxLZKA82yh;mn($&MKvOyMqNdVlOZeAiF0^75@i^xTZz1TU;7yQJiB zaMgBFxit=3yib$-rS+V@NuB8P@iA0l%RtL>1KyFG=?BV*XgJY3zMfvOo{Tv{T}Tc7 zM$J~oFj~b%M1F|_^h_T`hq%a-;^PY}%K8srCcR3c1rt=rH#=1r{eB?EN>8dx8zP6N zdM40$UXxhI(oz~Wc@J|csShuGVj#`xs!5(+k*JjS%BH0O`%o)oR%I`Z{d`SiARV** zDt|CvwOr{{Jyy9sVO}h}OwY#6rus#$s5|Sr`8w-(ASMI&v7x49=JgEH%U+j6E6;IY%FL;%_l1TN5bXL3V3XG=f0PVdvV6KA#5L~_nh0~t>8Q(5@n6$;PCQl z^ePK1+bvl_+rt`Zj8g-m_kLDMDhb=x-DHT>O5QO`fz5w2j6Xqi$!?r$zN6E=KNynA zj+aqku~ibj@PT~KtKU1K<^0Po#8iN;gcEN>7qq)s@gylv{u7YftHS0g}DeVE>G4e_oxm?8vnn()h%+P>D zmsa{bHU^6PvgzVGGBkX?8s2=c4OPxfp$|?sLaAmKcKlffkPQ65T?rw)tgAB#`Xa|3 z=w`$oJ*j~uaWUxFPnQ2K(^w`!f%+m|>_ps+y}%-jvw3wwiL#I(V*b#yW9*Jm-N4K*ZVPw0?G4uM#Pcj>xS66iE) zr|4sLGd|TTLeo!!*vK_oLtJ3l18}2ItQ*x3%T;9M= zunU1cBjhnjw88vEUI8g}I?K!{o+fmaP@-c$g>wlbhO*n6I*8$>D$>4j8(KEx(3v2N$C_^nl*;?_CHMf3yP>&(Pd)kT0)v#Wgu@FQ~q?90rv9LM&ZkB;gqHn zvb$grqkc;U+RDF!ZQLDM zW1}_}fZLlwE}1Fs)cF1RX%nKcfm=4XlfqlI4)iE8%4IXk2LdiJ{WDJdWzl9#~|Ibsxs^PLI0rJaR2&BuJ-f0;3Y7U?>tDm2uS%GmHkvdp)nly+Y_R>Djg%|B7 z@;A2SLdt>^cx07EeEJ@vJ2pHcT_=X4TsK{E)*urXA5DXZ!fc4mILBQ{oIrbw?gi_g zZ6e$Ic){w9MeuHc64|2ph1&LrfcDzjvXO7bL(rt*=_skO}8E3*g@Tjcc7Dr~^y)QI;?IhaPF^;<|8i5{DtsoBf zLw1jWG)G&38x?zpgbAW71$)oYyr^sDd_G4eP;XLatOL*aHJr|%PPSCghrW*8Nj}c$ zWOKP&Nd4&o5`Fx3`PG~;)cnFm;p((l+!ged`YN?h!{)TgPl`)0xqb}ovO0;Kv^Eb% z^gYFGohrrWO|O7qL5}F~MS&IE+k%NkiM;L7Fj!_8hs(ZO(f0I4UR1HUT)N+OD9Aem zv)3mP>y4d$?%T7ZBz%&E?WaIq#q~TrKri!+wJ|)lwo~)Dmuc$v5Pl6>pj?(WGNYzg zb{spwf8IWrFIKaGf)QKrk^Bi^FS*CK$>IsSH$57rYo;*yc{8wz3`*!)^&rj4pUz3-2`HFD38h+hto@P#kefx zB|YA0pC>EaBB>_FiD*y{D4TSJnV|T*BHv~j$$EL*JofGom~XX?Su?HJlAcHLFIS{>d|0B7L|*L*Winc$FV20Zm5oh)~jXwhRi(T;k}F*5UEYmS6n3p zRdIOV{1#E&u@h-wH&|Ni$|ddUBMPvb!Zf?q(!iq$L?-DD36;>sa(+GyQhY$3CzvqW z{nj(yB?}0nh~$%v61?A?P1enIq6Ws2q-N9!{CKXN?wk0A`kbFn$A7d(FSEm(-AyM8 zO}R~Ak<|w>&V~~>Y|k$B(Z-keZ7Lq_`$CP2W1!pWqZRaXAxYW27NUpuq}>)IL-hQf zmOjZdQ7>gMJ>eS&{R;Z?mTNX*cjJXbU)z8=EIS_!jw$2wilG1jBPyFhOhBvcXt{5; zG~?30h3T(30bWhqj7NJNC+6FRlI5?a!`?ai=m-lvYG7`&Z zLEs`HA*qTR-tM6)4oiq{_fm7?!HQsK-AYX!Mvy%VLeXqrS0K}$nN8@smtIR6PKVvR zDN@k!VA`gYa1)M5(_1}Hz{hnK?D`j9c$@`39<8$tI-nqkeiY^KV{nZBF4jD$_OV|y=ZkaQJP((ck) z*tD{bTGzyKr)(auRX+0}=fiw@zd{3N6%Nstu7Bu%{X;n2iGg`JpTkCK7IrIowBP|bsok_5nS&x6_W>uv$Z^Vikj+0cm zMidnu(+R7Ok=2n$$w|#2w0>g*uPW5pO*igf==t+pzTP*waH|Gxd1WZ@!|taq+-^cp zLJ__m)Ww@BrZ9175N4Ii**d;2B4g6*;mk<`_EgbJs^(n?zqef^g?n;Hn|}-xb`{dV z=GCCunPjVQWCV{5sMB51H3*^-7&J5&!ksMPJ1D?~`wP%sZw%<@dhn)OufSrfAR1O4 z!)xu|fnHfa#TE^deS&BA^>R#QYVtT`K~{GH^yHLr=6+(CXPEAoyuado%r#&CDm z=Guw+H_vo(DMqQ&8-Tk}%f%%KL&&VNbVuWF=oie&g#SfxwOP5iJ$5e4c=3j4)n?Ho z10{zyA`+1AxtbL@afodGB8-Kug|I#sQuU9N4|C~ z5cK@NNsB`ay%hO_t_nE@f3=#Z)*}b#?^Hw4b1HQYu8E@72NibX3voX9i64=?A&%d4 z)ipDHtWpf`&D6z;{}YJSF+bS`Bv{Q|jL-3dNw<8gH5ikgZ?E9l91fhM=6 z8B9|`L9y~IY}Y@f4_=z>BL@%b zLPfb0o_WVG3Xw8=qHriwG~Jom|1c7^m0cuLnx{iv^d&OnP)(N#`Bu9Zk7Y-0T*bg^ z%P`==W;&&30jybAPsbUUu=}UqpaqkK*>@X7*?{gYsF&YGn5Y_hV9q^m?(Pte{w zf*PnQ7PsTxgkj3nFL3l{A`L88=cCWQqswzUm`82iL_B;PR{fZdnVWL)s$~}*?K6jC zDo63A-D+Yt<_C%$vLfSFAt~7$gl^kx@R;@~XbF@?-MZIR2@RR}c7-O0YS+{4OR`Z^ zb{lLrkDG2bkV-0NAE8HO7BjBT)yO+{6MExxDOyYrW`j5B;p;X>*z)NzbSbEl-5(rj zjU)^50*y}9W-*C4b&PJ%IRu(AhTKZ+aQbSmGihEbHB+_nAi1}^1R7qSq?TA^AD`uo zvT3>0y(Na6;dQY|R}$0o^?=za-GY! z@7xS~?dGz#)Ct%=8luZr+tc64`7qx@pPi-s8T#FGNp{L6w(wv&1hn3OE9ELMGx&w= zk_Xpl-?}W+^>BccpB?03nH9g|Y_+XO-#1LUf0rD-T}Q1n-NDFxD%lcYM*W6_`6&}N zL+O)9D0_L0QG0WVXv?M$7n##IU1m8lJH~)-e>&!;1+#88Z(zO`OJ1)(32Ccbsze1G zuT$(%()Y)lZ!FJ2Q6+s|i;BXNchh0_3010>exEaM?7*6BH8ahMmct4uE8N^On-AZ0 z09RY&!KN5%?!uI5th~bumNjRfao`C5(X^1hv65i*w_T@?Jx_sGtQjvqb^;jhT?a?{ zPGD8TO3b)b#LqRnOdV=hS|1mR1H%RPYu8<$%3mLOLZdz>vm%S~Y`b=q$5ksUhWu4o}9ou-g_R8o+*I5 zn_aaFOL~DiucrHdc%w=Fv8sk?jl?j5pjv_*{#v%hQedB6bc&s}{kE zz#JG;@xr<*TZo@A@g`XveVbMLYm3h>@5k~znbfjz2Nbn+@ME+sF#e7VtC(5{8+4Du z$n8{e+{v7MWWSVmA9IgS?$w0Wa+5bk(xp${w8;dq-A z9FrJB0t>2bns<1kebpn_G<`40XKyh!-C9(_>=2oKt%~hxpT*J8L&UWymrnK_=F~TM zlB-o|D5F$GU$n?#K=KOq?}jL>*}jB*VwzdwdhZX0O54~hJZB8Q1YXabCnXq9Sj~=` z;9jxu^V!-|wu7x#i15AGVRdf#b1F_30=$l&xhx94p0bq5yd&p=5SRIblx4!l$`K;K46Ki(_8YYQtQt z(@=tUV`ot3w8HB3qYVNr@G>S&GU3gBIWVP$PkFzHZP@oD7>+#6ArG|YgL2_fWLI*; zfAAQcAu>X4`|by&oLNwtIgvZ;B!L5&@*r+A6Z<-!!I{lTP;w{<{y)6s^#8|O2HjT! z-FY|cUm7;jq04W{Mc*Qf8>hfl%i2Kd#Usp|8PeF+upZM_`opy&6*MMY5gxAXC9H-B zKWXYgD1RLe=YMS=={Kt|EXIv`7pGywqA#%i$1$i-e1cW-PT;hjpxNXv=yr7lZXUZA z#&5d-BQxAx0|JiSxrkp7bfF ziBn=fgMM-l)TA`Q)^7zMbu|z-XT_2^E8fuE+qZLK6R)A1eh~I_o?{jsu7?5PNV0!# z7d$n2MPtbYPVrGQc9}> zxvH%|3{Kk9`$@@V3QJOSDOKNqJ2O0(6(eZ`5onFO%UfbFFC`yD=$GttPmWu z!@xX#BO1^9O|6%lfV^|3S=up~>2puR8&~yMS5YDUzu7?jOMAL4mMv zQ42F~zc_v{EhFbD_JO4MKMYbmf{!O^pk9_CoH;f?122cdSnX5jyNRI5n;Hfzi}6rT zJ}6{PX9q9F(>)_4>>7ae6dZTW=H_k81*x+|Ev@qK-^?=?Xk--F&*u5Jm??>TS=7XNn!M0ZZ(r;Y@%LkEoTPZS zqu-~()lKtQ4bxsM6-&YCVzWS6r5c9xPg3*a_VDAj8au(moc|)Xh;2TZ!TmgA$S(4n z#Jg_q$9-pxaLH0f;f(zmdU>uidAubWS51n?ne+W1KRgyZJ6>{HE_0Z>&zi}QH;=|G z#ptXc|yL59ZDpMqT|db z0>$qFSuMW?#vKraO@TZdo+-scAWs99KO`p&_km(XAB37m;uJRxP<*(U{@DK%>XTMc z1s!YnZEeIT|G7VJ;Y2d>wCGF@%!w_YedONz_0io_lOg7oVDeM6a17`x#)r zT`Uw8$l=1l9CBAp6jTg+5ZWSe?88CMX6<_B{9=KhJHLrO4VXhTa&2Jexp6SQ`V|GSw3J`z;yNsbR73|ypfrUsRN%#ai|imw|{AW-!==}(vA`NiGk!t zcYd{$rvuo}lwpG&RpC;pfAq7$GN$^DEc&q^lx@2 zeNS>3anp04De6MJ*f6TEbQu(a(&2RQPg>ji3%X@QsYk&=c2c7d+hh&Q{t{KYP5u$6 zGT9ExAsZ)ZorC=?2Ox@yvJJ~B=*mOq@yH(?dM&7tiBVR9c^0mi_--DqkvdNri@QL* z*cpzUH)m=$UV`+`nV|S}55N{D>p#7CNTKNP#ekJ0WfJv~b;UO83{DW)v zGH`H~KU_Xg1d2-#`0Hut6c=E3`Jn~ucVCGL-;{WZ3FEP0Of+8TC?HGz=D~nZ8MX-* z9j*c1Cnp1d+DFrBXoU_Bvo7T9=5t^Gb`tDcsIHLXO`*G zE*)nqw~HamKRtz?`$G8?X>YQ->i{ecyGOQtNXA2Ie{gT7H5Q)mM5k#p;obWVa^&%9 zJZ93(+2-B^6Tej0efuaGHAT{<{hFD+=OrxEYJi*JLRi3?6TM+$m|tp$-g+*?Tbrei zcg+KLtNPl(tsglv^I)=l-}%~i#z$~jo;K)Ir_-0K4^ow*-8MY|pa-WOU4?(#Vy1aBi95~OK4Q?K5>08h7u(U51ZfTpr z4;+V|j8xJ2QYsy?s6@y1i!{C760kcDH5D`=ATADyzFQJpQb1GRzo$V5B-r11a;%7* z5AMmo4jb!|VDEDo{HHgEw+SnR$L~DIg~vMVHE&D4#yT5&;sS8eYF9RRU5AFCFtZ7!hcZD#xpAuaco>Unrn#B?-t%{jiVY|vv|c0hS);tHT7k!0|oJin)?igX*RzyYBwPJ)=QgFzB>#a~taOQB{>ZAB=AzSQ9_h#i8L$@hsy zK{qEI>R*#)+(gHf3f`rU7vS~12yXAI(90_m<#$rF)G#65M-IRgon~w|3xk6AQ7n-a zr@={DBsxl*ol;d#e>H3J&#E2Jyyr8%FOP!lHA@L&yNJD{HjzANlH*po{e-hYg|suH z5biZj7v#Ay{MVS4+Alw5Q}@NXyzQiPI_b(uXcQWQf#YVfXgi*77U&w!42OyOn<)72 z++0E5Qo-E3^BZfP`SAVUuYk}qVfNAO6?j$$(PTiBPmg#)xv@vVqfg3cNF(7}qnKkx)Mc68C&ymeve%C^42; zoH;~xT$f;9n5TnwhXwzt*^iGJ8LB<}C=B1P{>(IC*lQI0?gVk{TL>l{9Glwpn^$T?--NksLBlh=BX~F)~rw0$p4Uethb?zR8p#>8d z=@s9ZnHiT!CN^WmUJ2f{dNCVna0Ke^9w!_BW((fGnebHa7zUr73K7`pVV*PT3&m&jQnVf(eKVi`_pODww>%41rD?M}trXBF_yBv= z!ilS9#aSztY#OW@4fa|-%&U>Lbm0RTcJ7H`_?;xiwu&DhCnx38qTYL;r2UV&yt+aI z{65iH>0`;C$$K3 z|8$%thVvW9x5xWY3)AS_ywAj4SitI+zZGQJ_|M$UJmvnz3~-h|%s@DHfP4;j#CI)e zpr2!3BU5|<#1fy=ufD%Y&d>#jN$0Tiiv}8Ox2DIXI%}7%IgA-i^SR?ColO3RpI|F( zL{(2+gBSBk?0LH%+$`T)uqCJl`VBte=^8EEaA_1uwNDd;GiDH??gQR2ndIcHha_CT z9?D+N2FY)(+=9I!IIty_>+RZr{rY`mOwe^=>XQf>c`wkbOq-Zw^3-emBkJcGjEDe|1mk_lXK&mowqYEr9v&lF=?=V4oH9DK+Z;f2{$G)k7m z^jZlNw_Hw?NDD*@TfnQEvDhdQMAqitr7itNtcbxpYCKk$&S*1%pm#BJW$8FRgNoKR zJbp+ki`&Tc!ha-ZOfVlNHl4RQ7z%RtM4)bY6cRBVYWm_FL@zGHu=7hXGd0z|qiyVfTvC3>!g)S<+VnO}tD<>FB{F(&R7cgom*8IO+vKaDj0zY>{C%H8C3n!BiLcd=456o7sXP2J*hY}gV z@O$JZSlzKB_k*tqJi04ziPs};t&j%1bN|iex@DoPmI|`W!I~K}8W`;>K6N51&3L~O z2kO~Aht}HFptajs{H~h?o4*KQ{f#-`s46VzfctU%cTe_;l}+7`JG*CEKW+lmkOdC8 z$*tsxbr+cvUPBrc?sMAdv)S1ned*k!6Wrk`c@FVh3V(i2B76vNfvt6KaOQ?=xHr{; z`m}!{iEs=uWPcHt4hNX!>gX3R(LUA@I907H?W4mGx<<2N{QT*cT{<17A$wr=L_VF=mdou z?2n6+(3Mw$bt^1c7dOH`{_u>14g0e$&zG=YXXikurwcyV<^s?3`pF4pc~c2>Z5Yq%4D{_YZx}F2wwFl*_I|fQ&zDmU01(Xix?k8_mPfyT93h}{x{_vZ-mR)0T(x1P^w^*>#hzNvsTRGy-LKga+K z3sTV_HjN<4rFn;DI75@qY*bPNZ;>fD!=C%xqtmK&b~4gMZf+8|I9Jq zAu${OT)zZdmpe%s=JBbx99tTo%vzg8LEGdr)FXT{|7vRutcd%Mo8!!~1^yqw-QoeP zD`|sRoqSNT5p{T$XUCqXk;8gvmaSJ&1$*y%IBd6%Wc@si$$tgDt4UsHo$t)}+}wc+ zG*ZCknm@a2zc9Y);Lz5>fHoah#0&pyNwJe1D}8!9`*-6xviOe)-#;SAir85~EX_%0lo&-7{SKP#$(zHj}kK4&Y8+kdJls zB)SQ`CR?MY@#+~L=o$Kr-0(R9!`~OfNRl_deYX)mD?ttwziXgb(w4Nzp(L;GHw@qGdqpS zLD4jB!V0P?`JIZMvIS51wZzHPj<@C}V$avD^lqdwt{X8$gGF1|J$k)h-c*efm@o9A zrUW{tF9f|cnqX2Z#_zPy#|v7Q*a?d$vHoSmPt;0;`@d~4$$KJZ&tC-f-ucM?$i~p! zi{Sd)Pted32d7)E5MfJSYC889k=&XNzt$%};6NN+99Kl`9frxEqBJaczaY^1w5s$7-q4?+o(sSRtPBX{KuL7i6W}G(@4^>{r1_r z$}l`El`2ZG7!dZCsCY}WOStuRMj=0$1#eZsdVeOi7c=Yv-+ZP))sBB}*2Z|=GK5=w zf+jjM3p2Lo!8+B=P~|>_%6?vefvv{?UjmFdWD3nSv&opy4pLkH4m<4skgX5masI-UKA!MZN-K$cw_W@`Lzq z#shqI#~MbA{t->59B_R3juWl@PC7jVo^AK>tdyD-ukIkmzn?ytFXMK?S~VN^XqZUm zUR%pYB}uSCOHJvU_y!OgSO8UfOYJT#H)k)Z2t((WJ+#35I#@{eaNAmPASa{~_6pRD z;;+{!FLfN(3wWcX(l&>S;0c zSH1|V)Y_Q-g%a%jko`D5ErCq?90(4ZvvEVuGv=GmKl(@NC%9aZ2cvDT61y@~p8iq%fn_r~1+)Iobnr?% z>AZ578mx;WftG5-^OGI_RBSt$v8$fzzi2PuC2v9ShJ0%AAssAxJ_tN%4Om-qpFH1u zkf|NzNsI4TXjKa4|2(#(t;>7ikHjQgcq#{WaYxw8f@^bnt`=W#HjE3lihF<_wZp;m)h4TCBh~_e$A%D z@qACz9cb-Ey7i|Et7eu#F7}jSlgJUqhD*a~9~a`sI|7#ei4q!C>IgWZ?Wh~}h-=zj z2Y<&rB^g2Y7_ImFNN9usJRCX+ReMg;(D5lm@s<^9Tc3x)FWT{!K{R|;`bO`Sit#3s z|HC_ngJAWXBn-_cgRcRUcu&%S9eYLCEhFme4}%t5uy_(#WjPgA2e!aHFJ1PItQyY0 zQIA~$57EDaub_GLPSj{Wg{lI5K5(5Ecj}aY)8P||L}eOrsqzJ`*b<~R`IDluOZZYa z4u|(!$41Qi(NL2ctlPfB&{A&|s5L31V)Lno=BfCg_gcBw{5-`c-UBGS^Cg)EGJ}1|P zi2?5TExHD>!zyayR>nbKQxE<&`2m+>_Sf!vRYdn~eM}_e=VAKBDeTVdAiQmz3b&76 zhrYCT)K}b}9sGI-uO%#{O|#7KL}xx-vF|L6tn$Yf8;h_cz5>)NobmkqYxdS6%A}}8 zm}i&e5vvy+_{De@t1OprIZ6FMc>Y_2-P4{3rK3Wu%?o|DQy~W3^8{X<8830d zFJF99Zb(yOf7*w4OvHInp6okMQU0d4BTiVHhn8cTiPyPNIO01NpWhq=iQD>YSLGq{ zPBWLT*&~93{_{|_B$$Zm#gJ#?u0wwwg|6;$aw=p1KAae$rBdm*<&hpwPETgzvmc{% z#aXJ;*Khx6TM*qXTZ(as_NYC*0L672*e!nR;r6Ld+`>z~lymK^Rj*3rcXBt-Bw!u$ z>iu7wF~$QQ?VAL0AKc*d9vvJ`*@5;y%dzn5Z1|R6hvUACph|=qsMZKW!G9vGVx>Km z|FfD*J-LrNL;f&%x(nGk$_?~q#C;OlGae&cv-!%UH_6HEHaM^CCRzU4942r4%dGWk zARYTAk_7EY^!;0aG%^y(*J=^+Up!x{qD#vIW4PJ@ZPZUOfyqAx;nw?$G^SJp{pZfd zEh9rjZqqIFnf08+c~sK$WK+~Xz5}8oCc=SN${3y+OP;Tiqpu!ZWxDU>!TvuQaM|Dz zDrv7~3f+$|2d~s~#)fyuqx^XAi3kO~*~dxD-ffJ+m;+?x22;*rg(}QHt&E`?N{IU5 zAYxFi13P!=K|tgwSS9rsyyS1-t2iIjgDJ4u?G3l^@BT6FV-C@8`A73}FhXa-C zz-3UDH?T8;I2#QVdM#kcpLBrWcxQ#UG&p!X4&OQ6tS#a9bJoK)I9qu#140R)ZTN`k z+g$=gu(#v1(FC@(pO84gy?SKQHhB4x#fIsY>^JQm$T^}&bz6manI*f)aIT;~7S&=O zny_$v@&I<*-oa?r7d&q2!f6{@tf?}AO~QsarDz6!{2Y&MIX4AoWjS3DV`n$^pD2lu zNkkv55IoiwMPsH;fF<*@$+cbPs5SXGv+h4D%u$wP^JyHg#dh?+nB`Dpy9@_={P~j~ zD($3SuI62tG}`^vhmHGa!>@=_q1*WxY|PVEbgjCDkxPyFz<_d^Q0XG5A-~dRPet(j z+eYTa-UBr5)mCh;JVK46-r!?{0cO*~0^;}43qM=SgXPOsnDlxU-riRM({9N^Exk{I zmno90S6k_BwIu7pylzAQ&Af~A z(-O(Uco7g5@X!|WiNvZ*kf%AEX|J1kfc_ZhrR}2&1=&7YJl>VytwI4<{}+tCni#Z1 zeG=s7FJTjp--EcOMZBLbz8~K^Vx zG&R2llNQHP_M{_?bgZMVC;P&_W#g&8l?mA*!a`B;Td2Qj2I^15_{xVm#Q(tuI_%O3 zuadh-V^cjEh70;(absxNH6A@KE79~Fl`wHI13zX`^0Rd_J^uOz5mOt6J)&pf$TX{;AtqKDL8{h zCt}OiSh~RI4|$z2mK{Bt2wj^;&^pT$yt<_Mlq7q$#Ns@9KiUL*xe!}doW!;YGP>8C z`;WOtGx5jFVz{H*fU+j((6PP>H=kO~{<*UUV@%_q;QAKa<;CFjADVFd@MSJJ?FG*E zE1_;%he>|sUU=~Q3l2#ygZl3(G^|ya?-rINUjM~`{0C3oq-72-)pVV+KPk%x^PZ?Y zs}xg@N%P6)6u4&L-?(S+FLnN~8K>1hqgOspff3Do@EY8MiJ7zM@Fxx0K2{SOALrJ- zJE{d!pG*aj%B!I9a4LUMTbew`a%FD4*?^aRU&Pj--}v7mIqH(C!rtf|$H+WghKnyB zf!kLN!A>w+xKvlq>d(=F!oxD?@79a00!>4~umj7I!^rZ8U8pd!9DZVKeIHD@>1W)Z)#2 zw&08PiL8Z)JfFN&9)6dP!&BAikVHM%7u+^*zOb9wAzegw%VjYNBR83Us+&;r{WRSD z-xBy4a}k3*7r?DJRd9y=V6UG_=LW50)dsd;_b&}1U1J1B2{WnabwzL*e8?>GT2AE3 zCeyOR6Y0~9vV#72m`JqFAbg1!R-P;(%Z2n&A-k5DXZM+@<#oxPE4B1Px<8GYx|C`& z>o^ldCA4dZC9-@4W7=g1TW21`+wv@F4Ejr6e%MHs1v=3=PbSkJygS!D9LCw-ybqJ- z|7Kg4kZImX}O^X5tft(muFkhey1Q`ID9hHByn*3TsX)NyP?YW)3D)U2;K_5 zOy8v2VEU4FXgzcRK3h+J7`LyH z^!`N}GiwqCu0KcLzf7f7hf~1A`3BdrxeUkeJ_Ns+0iv+&6x{WiNA%~tC+o8ciP%v} ze<$3;iWKfyHr5W*EJwCR-R^<=ToWor1X7007CZg@H#*^~er7ZF2E z|CF;Hy(Ns@c4KH;tw%r2y@a#s#_$>s&(}Ij&xS|FOPDoMT3mg7A$lb3rlpc~F!^yO zZJB6}uME1-wcLiRF5QL1TZQQ|sSzY0SJ2Urwac(l?kzn3*8ArrAvG?)fcH@1Na5+7&V!@_lcuvGLUdihR3 zKfOc@*|(HFG};HFsw4Q(tBiJZB=PlvETVW|AZN6ID`>oirwo0W*(+W%bLXey8$A#7 zx%(Uozlh;LLM)vwYf04OSyr|FAW0Ha zM0}kI`On4!O-CHb?W);u#PYddU(9OQk?cv;&m_=%!>N4mO$SKq8DZuem$x@CcE#fD zt4Vn2YS`lE0X?~n7^OMNWc&-SZ4f_=CnoKH(zoYfs4{`ZEyyL0G*4pK{RDc;%$Bu0 zCyH}l%;Igv>VfvA2N34B8BUwe#MO%bpsw$ZAm=oaP889>j`vRRWZ?sF44%OGwRqtV zE5VH6n+R@RUc}XQE~J&sBK)J28XV8GfrgbkYqPG3Q}k^>(D2QP=&7nWY8 zPvSbc88^gup{H8~u9I@?=vKm>#(s9~R4d;2${Mm{qi)^Fb$X<7GxKKEZv5t^N(Q>y$bSJMxca69iUlU{YR18=N9S%ZZ;OBoH!oUavx#)B(_-3` z53v$Nz<|q$uqhD<5b7q-V_W91vs*5Mb-WV)?(i47V5pI3YY}!*+f+98=>ynn7Kf7u zx5A#anOs?hHg@@bpdt5Nd9E%Mk{ zdmJCUk)sKkgA0ZnDzR_DM)!cO^941#K(WzG}1X`;d*4#de#<|;Rv6UbrRPHtu zr3Ar}rbzb3!5hHn>j-H0ZtCF^PS=ef?OGIqU!N_*ZKI*M>t8i~J$MSAwEsn>&Rf7h zz97wya%tQFU$7#DpcyO3q^*=^Pi9YHo<-KtWg5*)#(FbW=R!B5m0^m_R7TKCx(k@b zvt&(e0(sJTrZ)Xp1N|>=kY<_hrQbKM<4y^i!Mlcw@H}!pZflf+E#m_quSyJ8ztvz5 zI-etZ)|>@{CTZN|E7-Rr>duD!mjg{%BlOYrlVrlQ)Y@(KgY;9y5?bDK4|842&}UCQ z)gH>nqq%Tzd>Of#Dv2gvHshGh&ak*@92|eP2n_^L zP8%#!aa>v=_32s!##+ni0e3l!(msUSVy&QX(g74-_!$ozw1IQR;;`V57iBz!VB|&{ zgl{k;pZ`|VfZTDsAjFVM?DIrxof>AD06khb?HbrTg!|zG=iN%Ch_{d4bA^e#uldr zf6n0ywS26Iw4;-S!|F0->>^k16~Vci6Cg0xA02i*W!~B=fRL#z))};yDBzZ|AedjQjy7R;5vRdg7uoGbjli0_$=lMy9KUi zy9-0~!5tfRcSDi=?@htTUv|Le*46CNSQVDtauw&+2(l1rlCdMT-roASJnz)8n_knk zz$UW*+$HdT{ki*+4qBB!hoTNK`RR!l&ufAyxy@7uPDN>nWbS5(6?jY+tkv>Z3`yY! zp{-RJ?g;9L6mtt=pRka8b4r6Pf?XUY^GEEBoX3NMdmpwmDc~jTb+EeW1kb!ZNJ3wn zAw=DScJ!TKN>*&fXM=`Bhp4cYW>UnnDv8;?DUaCRmZILTMA@?giqvIy6qc+I!~Vy= zG4$j(Z2T*@r}ZmvrO*fa{XvJAS4U;Yy{1?)cwN}O@$hLfX8UPw*MeKrQSLbDs#Js1 z+bYOtp%%)EXprggX*6lS8L$eE$fU?fN+LyxoX=$PpPv#IBq|e3 zE>tyX3FDfz!`^hgz)v;i1eJ9%gaM&fLKkm{)MV$3Hlb-T83omy? zL;bn8NnSjbI8JuoZ zfLpTPf`+^gGi_WR$nTKi7ypZ*epG|6nQFwJjts^PVrFn{^cH{N!9MVx-oQPGE#;^7 zu&_zQAN{Qz1+y(D9J!eYAwg@vGs+#O;v<+eoWO4CE8@JSJM!}@-RV?a(Yk9@U3lGT z27A@A9oL9Uf&RBRR(IboGDD|}21e-eWx=uhZjDd;3s*NjvO$>-{S`%5bx-5TF+Wt> zw4b}uHUad_CW9`255s* z6X|usE~w~=BMVo)CO2!TBJ5G|T6kVcgqzK+>Pl3c9$6ZhTFeET>Ycs?EM))i2vY#ZXi z8Sy7|tnkddRwlO75>?&2SPyx7K4e=zTKu--Fa7+^cYm7&k>$TwoqI=!{|r%nQV)k| zE`G4%R5HEuz!7o;Gj8Rf10c2L78j^o1w#KW!LX+SGg73&yJyw0`6(P9J@E#(8=8Y{ zrXizbSdTiZykOstYIcF(yFY|k^6xO7KVSO`cDpWwug|3U?5qmj$D@djmOruG{nm@s zopRXWyH7i``zO%vJ8rQ#!&6xGDXXA(wIhbV{Y5v%q@iX1931Q7Ey!*tLigGa;o0VDfj_#f(M~qJ! z!Hx=1_R8z!eB>)FzW3#PZkt;&nG+(4na4-zK)D=mry+@d(yHJ|y9U2E{2OEoc%ST? z?|8>O7m9Oa9OU){klL6MYKKpRX6WCsS1HUy4P}Qf$;e1vlnj1=zyJs)6>WjwH z%f&++NiN4-6+6Ldr#+u+cbSu!Y)q|J+-2U1x=@{ySCJTm^V4^W3w)zic>8Gw3}+Rg zsDm`TO|zm44)oAfwV60M}Ov$90jRKT_yA^(n9Px{tafL_y}9$D}K|fV8{H z^81P;u;Tbh;Lcgsj`3GSv&d%dcg{hyTwX)uF4v%TZw#03WQ%VU)6pqh9JdH`t_Le; zkk5CGc|&I@fye$g_0Bs-I$u3U%jcGOA#WdFC!GbRzhhuemKiLo9Ob%KET9d{4q&xa z*{R*PvB1TKRaV!5vmvKJSNI`Kl*t7_o*x;wbrQF1xQ=1J4!~gHYRq&v%5a6}F*Zzu ze;}AKTPQ4t<{%fC9q+}eMrg8T;gj*=g)HdbRR{~p578T<%@p+&*@%P1=(qC@Zf`h@ zM6;MCz1xLG`!3<9T_;J6pl%koD1>_V6v0leCp1)VJ5~CT;(>DPv>6?-@vB=mP6@l1-Tz13gSKHyn2rc?D6Nw$80zBJ+=X= z3`Ed3K7h6|li5XK2}G;#HaVttkM6It!Is=$uozg2NufQo>XacDpiyZrJwFBHBh2kW z3`$|?c~_qNY{V_*`#};9V%d-o>ib@xlk_iO?D$-IK>9XltjZvc@fr-LbPm#mrn9Cm z^KitikNG@S2rAUth?Hg)mV4-+vFJs+i@KT6a8(%%-roVy{8GApdnW7}iAR|_K^Kd^2kmJ|bYSr@%5Q#MTknCSW_L0csiPUttHYrg|HCVtYX*NQ`FIX@Jp^R=knE z4%`2!@YVfZU_0#{?z4Q0@6IoW?*lW~o5s?>?!FJdoCW{g86{}DdI3}iqd`gM8P*AA zbHRc=$m-sLnWM=^NHhGxLPsKy=<8(uuUfqF!Jf`j)x(h*N9w4s8s1Jk%0`%r^O0c_ z*}BQk;DC2CzvPK2zbM2N{yTq@3=KPg5RoA}9g5j`RZU>wJ(b^~Jc8zOX4rkb9sgSq zhpoA<=nD5sP!SPNoZgh;hyBv*#{gA!SKuA&Y4^j%SD$fam>j<%$pMFJrm)*=`l*Lt z2kYStO6ckmhjQ-+U}9uEl%(9E>~0x!HarApB{z|2nrm>ImL>jM<-!;#yg{RnN?>Oq z*cahm0a~07$i9^X#o_E)$;A8I)7#?@cYIgUMK_;A&&wllPwg;#92Y?$PK;)a2IIfh`Mma#_Wzx?oc@1#OVP$Q zlvv}9Wj}oR)E((mA$c{5*tp@L_h$e?4#VKxJAw*MOi* zFlElL_r0!LDnw_jta%&v&03d#}CLb^U(l`I#!NImtyHNCeXX zPx4)FtFZ0W7g$u1Mp8A5nPa5TjatIE~Dw__{rm6?lMmfeP=z$sAg z`3v1&t8luz$6-)!2Mszmk{RdPP96pEF9felkQH;_6ETg6{CyRs7EfU4a(CP+F@{~G zGewprYZ#k8OA_R{h8b8F57C)R(d>vC6Pz9b$9D_R>*Hln$*#BPv+@c3U;37||I@cT zF-id1`%Q87%=2(h=PcYGRe|2_L%9u7S>UWzNMcV7(l+`SZx6_G+fS!sz$HE(2<7B; z^=(n+KQBDA??1S8(;3Hadx%HfG{Nj+GR`rKB!^F*rOMmea4fB&6ZR>Afn^F+@H>ZP z^B)QA9OE#(+Z;bMtU;aCJ!D_v85({<9_vhQkc5ZPqBP$MI68Y3@A`TOxA8rl7BmxG ztizzGAW0;@cNN^-AB;nCMl%W%r+|)6Cp+HwAk_UV#I^CWS@ZBdHrS>BFRWb%`3sg| zSIIN7Bis;G>W72pC11!Gu>jmR?m`h;i`x9VPYS-#L#1`7GyE(V?ad}z_ZMTE_cC%d zUWHrKdzWhQ4ELeMBSG)URyeU^AF&dR;d_%W@O@Geww*KLe(ciWW=w9ScVeShr%AFX zW9R_~Hy4oCsv`XMYZ8nJ_JecVk3;U?EZ|05hK<&t&=3|&<3|~DuLL}Uq)v$(%Fy9{ z{5XbNx6h~ZXZVY*zdr{?-WMS7v^8e+y{E68YiUO3B=~Xg0*z^w!9|nYapb-zGFIuA9!rQxmW54e3=4(8=trs{i~ z$up!qk3>76?<)gCsBVvB%8J@Z@ zjhy9Ojs`K`KxN|-Fj%_+3XE5Br%gn7B95Qo8FN&+{sZ;UumvlHMexfGk!2qQFOy@K-Jzb+5gGIR!po^;;7{dV+Diyb&JW_!N~un&!MqAZ2+M(LjF=l%FDa zZt5I3AQgat)qLL0{uvqgVTf@*9@0$-_u=e`Akg@H598IEY0jUWaOcD(NPf=e$r4f2 zU$hl-Re5&MxY0O%);UNqp2hCD(1y{svW5FWg{wKgknG&11Re8E(-X%lsDm}-J=6|@ z@wdF;pZq~?wO1^;DW?v^MHWAJ2&u}UDm*wj878$}W=H$-9FXTX1^fM+|D$0ogJ}w0Mmp?@zX*@(-=xm3BBg^p7R6QHm7F z`O33dkAI4CHhm*m6HgG2h$krY2qU>aGEn7}45166=#6hD$^OA>@aFC)d{!1G`g+I; zeTQtsWuD)uV5cQs{v(T#4I#u3W?`acDWm;Ik&8YUfq#{M(A^)GLNdQ=TDxHc2;NJ` zwD`5`BSsI-X{>`O=y#fq<;ACP>}VM**t?F{j2udITTh+=7xF z0y=M6IDfBRNEYuM$=ENghj#ZYkk+e!Q_VsqNx=dOvaK+3b_3i!>rH+~^4Zu|kLi-g zmn89a6P$7J2IJ5vWW8A?mNTiQyu%DC4A@1$WhKuv*TT zOK=aOnQHA=ds~C}MGdfaGE-rp&kY>*X%qYo4S|x5Xma6-fICHsFg#d`^ZNA(J-=6i zo}w(BrmTb+XY!a!athqHsliZt(uk?bO2?s-vms9t@&4*n+y?z;TuDtbDz!~v3hKte zLc^)}%+C;H+D3CP$!$DhHjbIDY{{91{)BuPKF9E65Pcqv$58E8km#zw-8$bww-rc1 zh3jFi)AkB+`B%!U)8Bzr=6!TSmoh(dtp&ZUAHi+vc`VIo$AYzK>=J%2W9@Jh=b0V= zSH&2R?RVqUcOHa;=`u9(k}~4vU^-VZ06re{#|6uxu;HsbcuiZ3e=`*5j(kO)`_L)z z+@cKbv2&O>=BsGMZ)47JWeUCb@D8f_oMpWXC!^B7D%`B~O;nTc8*cj>F}fDF&?~GS zCv4gYB`LpXo;u*{7Bj|n;c{}c>=`Y{2?KMLA#{wG@&sa$*>bdz232H&eDKV;nPef*oh`?h2>pQiz@}`>}A^P8u>8g;N$a zz~WnpP~lQZ9wbeUf5X*Cc@2+CZb7Mt+ zm4opllf`#qSAlhIEVW%Civ~4vOzNJoutH}8=&rj2_GWJ}Ht_5NERe^%b}{Ycwn@y$&5Z zcj2{?8Z>R5NjyrIVWP)ECgO$^TM@?Zr>$>-qG1RfY#N4Qp$z9-nF>4i`|y0EBJy$5 zFlbFRVlK@$!4(^g(0pMb$fO2<|KwmO@)1+(=4^35`ds{W$VvQYrWP*$or=A)B(Q?@ zhayIT;mbNmr+6W^TbqyzbM3)jWh`?)RG+JhTZ9R#Tu9Q5P3U|=58A?Ja5w8dlMpZk z-&I+xZunQKA|rvEO|xLE+f3&8)b;SpV+7;*GL$^t)<)J1S7$2uoQdamWhQ;pU!2l0 z2HPKYvo29z;OO`7n7q4zZm*~U{T6?SJAWOj`TWJFQ-7c|L;vAH z2J&SB8q z(n{|Qaxm4u2Ildd$NhCb$?Ijo7L?X=Q<8`okO=C zawe+wx4=rP9d;Zsq+C;psIJTho*4Y2OjQBtGkK1dKb>xMlJv4?sd-0Bredi82l>=mJjWu_!Ga1!_2+46#!;sAf$hr~v znAk81XKs(Baj`KNS$>5aIk=eSpD;z))+_uDZY-Jf&>U{LC^FOaJgK5}0?BClL2b5n zksHHhne!WsF~YWs7%xvD#R^v-PP&{;_GzHQ<gzTXn=OH;-MG{2VpAHp1px`{1k9 zO8W7Z7R+j}0fvc0RnuzP_$h^CPuzwvI-O{;{4KO4{eqc$r2(5KGrz|VB~R@SlQL;( zu)LZ}%Tql-x*?bt=>H?jmBNUR*=)2gzlJ}PYSFps1G~0zBp!e9j&!-2fqDO8xT5_C z3idcL!RkQ}*Rq2wb)AF}Gqys2S|It{eT%G0L^!9Uk2l7A;4|z>cb`2M!6(w@is6mUVy7=T~ zD9q|xi&Z1ulcv*g5V5HOlHA(Km9TiC&hqEjZ)?~b-o}ntSR^@qbsp|we?ZzH341d) zljMzAhUt!};QH$(G5D&&X;v0sd--`Bt{^Lt4K%|gUYCfgv;r!W+R<$h88Ef~E$Iwt z20Q&yP*#;i_X}Eh`C>WRopgsC_pI@Ehdt;_nnujmC!>GCMVjOJ3v{Ft0AkB&xvDnU z_B2s?a~tIF-@|mCx3ckZAQEpc^l@pz)HfH2_cB{pSvLoNjvIpZsRDGn*h}>WDCx}1 zBhjzBXe+{-ZU4CPnlD@rrjtTzn2O64X~QihO;S+@TRr{UrlVnh1b5pq%LzDy-SN4U2tZ8 zw8g;ot?3}3KgqulZcux-i)1}vVMoJIIJo*I$fQT3{&EKEGc-YF+e#epOvj&#H1PPw zBe=Um8t-km2u?m_lBByF-2b*4_g@PnMs~}XtlsxnH&96ayhx-4+vehf%`=&|PZU8v zYznDW7z&c_H|U;^>3I2GHha5mC9YEP7Dnx{#OBv8u*g9h@<*hQA{%+8q<1O$Y*!&x zo2>Bl>^D>@)E+*+{RGDE#MDMn8QcEkK#F`6x;%2`21XYOe>3XDGx$1@F4&C?Khxl! z+BJ;0wHbas{Rs=zj^TloMOb{=5L@+Jz$4Z~7;SPM_skS>)f1k=YO5NO*d)z2f}QC} z3o$uhIFal-d>S%mej)cw%$V*Ir$9E~B%SghUzEBf2PS??p`VK;L#@w5xH|ML@HH9A zeyxVyLM_;^%T;K;XEzz&6UOKAdxfHpr9{Tlgqc}vM5ey^1Vc@2xnYm*vG+%>!I(m4 z_O__s8mx1W6d%;FHiJTMxu z8ufysakyzJJ}fvSnBaL6w>N*s)d3-(dwC;dU2x@;2c%%cj5KoeJc~{v7gC>h9&}QR zDf7Ww5srRn!jDA;^tEdrRJ*sZHXFC%i`z{sSE)sZ-j;`IXE}s*vk7Y-0%1J@j`qcY z*d&^KEsVglpfjL%W;7#geS#NP%Hs>F$4mn?Y$z5WO){3S}mM=T0(QmWSHeT z-2{8oz|J?5ijI!PB^O(vBEJ(?@{A;njmC^iiaopLE1x+qGNR`-3H;1ern3P(K}wYa6HWnUs||B8LW!%0T(5pa7ko#>rhg{SY! zF^#7Fpk%4f`(j)W1fCs>){mav5BWApQ$k4e`G+7rWpVSaqad;FCAF>G^pUoBO@Jyhu|M<-Rv{_8{+7_}PBZ&HH z9OYTEDRkV#Q(z`)V_RqzNf&H}O((aZ(z45>i@&>Tt{#SKMkGK^Rut47?PrzxzLOXw z1&mhUdpdVJ>AC?^82#Cdc|SRlj<_Jr^9(0bZr49lxmAwG`uuS9l>78-p&PWm+d)0v zgh9E&4f@RQBJ0;U9HSPG7JdI11_yHUN#KVBw){j4_?$Wi7R%FsEm#O+XPcwnIB#@3 zSciF8(wtMc9waB{k=W-AR5aO?t1?c&)uZo_y}9!_`;Lrh7Hn9uQIT4 zb{SEyNTZ^%Kw*lsGN-pK0CRd(nc!m&NW8NmhG|`;$CXa9KOW46^m}8W#e6%RlxtU@c!)Le z-7$>!7#?RjuGW)z$pPf8=r~n9_zmL{9#e&X7Ie0K9o#hXBTjD}nMa;~;7$|I3+0Xx zoxfVlDviHHJ#7U_4a|@nESbPeR!E1?H3|&-Sb;QbW^mJ~w?x6Z0(NwU;+L#$QN{WZ z%y5G;YMCg@loWUG*kJl#jtNGlyWDLm1kNHJGeu4{W(Q4x1W+MVG?=keVj` z_jmbm{IcXTl_{{mXQrb#_ggYhv-z~>*O4V;;GYy`ihil}Q-dk{yi8OGfk z;=qimegyt|XE}~9ILe(~9nM(VPGRWgXhy&O31_uohU2VRQ%TQQ)hdJYCXNp`C^?!+ z*U`yt+nBgHe;C{UDDyCbO#dzX$~2Bnz(>YUnD#vfk-{-pw>8u8$<{A)=Dk>^e9AU% z$tx@7)07CtqpcMC@4shys@FmQzb9ThRhMzDZ{rj*j&MCOfy~k|&zU-9!>WD5>Z`J5 z2668`+Brtqlu|Wye>C1x#XOj@zUpgA9OG*!aC|f4gyWpNC0v*CWV-s3dX@3C*^UFB z2bl#{vE1X`?Ocb6cGUtsqpHydQ#sd-ml&_2R_GeKf?eTZhL6|2qB}Ql<#YKDsRKoxyT=Tyug_dd&_ zZpn8#p>`5%9-$1YL>aVruOi)}@Pa*dyh|h%@RXWvw})w06kyl$)mR-NgpXyHc)zU{ zb=((DR^>z!sZ1xwHS79e-==psHvKcrB!d!nn=hh=-!|Yd{w#~zy&hcm9wqmzC$T4d z%Hd7Hds^9}g7H>Sbh2bP(H$B_^*}6QH+N^fRrM0GbZe%mhxE_O1n~XWngBwwicQaRR2nN-Z`?!Pe zerOs-qxZ=*>_oNcymM#^%rmzX4k76fd3_W_3*?zEC4MltvKT5-)iLRf9x-}pMI*L#W?WO|E@o{$2sr=V>-5ZnEPNHD&s@pb1hHZR?q|$f zIOy62Ps^3ib#E>XELI>zP43ue8VhTf0HV<23USl(vA=c=_+%W1<$48V?a(~7!v7N7 zOgcsXJQHhf<#~u1jdYfLRwNEksaBHHl;pXXNeYi8gvj7Wk39# zmdHx63Ai%1fuujQ1JmOdM5ESUgI(f4?nBsKe3A8!SUjB#h94S8;hod)sXr4q^&r#< zmc>0jrS$zvbJVQ9Og44gzzZX{< zEoWeoOa3FHU_ddFD|NYuH72`2Fl$icFu922etUy9QTxy^R)_nRJenMA4#p1={Jo`r zt|%yIHdFA=1X37;t7AEs_v#visWs5PO?F^wu>+LsjA5StBQUPpX)b3UX=zof7;0S2yHk%83)#p{t!QNAy-8XVd)%2jCB783qJP|-RrB6QGJDR zzw%))HxK`eDd(Mqia3x^i4h6|P=ym2CC6H##Z1JaI8CEIID_R?EeC3IX^9zZA`b+qG+>bh1 zT*m6sqaYGU(3#J0Uo5!F$_K9ilSNPAjCnlUrkQB&-3mNemG)j!Wj2KKv#Bsm z$LJJiGG;|P@wLfhpPyUDsTej=zbtoL6+QuM&iy7*(qGAU-5``)$q`9^NTJ)k9hsWJ zT=3nm&n>>M3@2OUIEUMEuqC=5%5RV32Cvl6EZa^Y2?}KvuZX~Chm}?P`<|oa4Ock3 zC5N>-+a|*qJsZZ1iqFQ%ZKKiP z&NhjGGSBY{4+HIg;rQV{W$u^Z2bdl2&jn>H!Uwe!W0Z}#^LnAUeqsTqIeM65$!ZaW z%X9Hi?l|tAdNvI4bB6Ne|DdlZ0t>6WVNYl>-kx|8mZ$83_V7EnFvXSWtrZUs+lDWtSdg5V%dk8k}zD&&@&u6S41a?l0 zA|YOOY$)!<)QzERn^!M(FJ)14Mixefs51U;)gp8ID4N>WNp$fL!Ns$5olTkK-1InXm{Qhw<^3O`y?09ByBU7xOhCHceDfMulUV(V>LPZP9J=-E)Gsi zOu#u`B7`9)Z-|n(aJK8F8&Nv51k$a%soJUaAY<5ruk+>bV}~=7z2Q8D`(4LnzPjAA zKgX!z;A`CWb04huas--`FJRehIkZ$&!|&}=QP%%Cn!R7b^htEMaX;+gXX8j}TsoSO z@`;1(HW#3*C5--CHG+Fs;fS45e}t^A6O*?499XXTN-Tzdq`y7$$ZzA(?C-~B_yi_$ zH+I+4mZt_>v+{53&h{nWx>a!O;B}B)5Kq&Tp1{*b28^qPHq&!v6}K^cErw1~M6*x- zL{nC$VLc7TgPJ#BMTrEyS(?FKsktcocovp5DKkf=@r3CImaxI%B>Kucr58q>B&s{p zsY!P%G|ixJ$-)oRRwT0p25s!G{2HQ`>Ww$YT!JOG-pr@PcZiEiBYAk&o^Cm^6~^`& zVq$SVXqVOCf06oJVUssq<7xR1cOnjG4b_zHt4-Z%LNQaq@JZIa)nX zgbx-^iQe+Z#6Lv`j$a!CAM--tKuHm#-AyAKRCV!O0K)`h2W{cc?%pML$jMm#-fmon zf7i=1_o^l`3;GY>XWnbOVxAi0dV0fzmI>U9pm#7oK#6%a{UQw8HkGL{$s&76DE9Q$ zi@g6+Cb(A%gN7|*a;&7`{Kz!X%nhBQMN9hNR&ypl>mSX8A6Gz;%y%*)DhV3$&f|uE ze+gs2dqb{t6T=7n82)J#=KeWF9zK3gledh(7~w&5o@&E9P@2GJYgyLz&jL`)=oZZ? zt;aF6p15h3&U`!F06GYlt&pMYiOr^8^^Y2I5LN;AcUZ2#MNqT7q!<6NH?F!4?? ztoErPr`I2czS~Fek1xf$X*`x$x-WP{u$f9fx^_ zg4}~Pwn0yaNc~w4#Wh-B|Jsf2Je5Mv_g1k-Utd8NKDTeGd!IHHZss$`cSK5=5h6c> z>rk^UpVlT#2h(GOC9v6dX*`=x&6i*SSGNaCg$f%i)d@uk#t=&Za4AvdhB`D_K8^Sa9YiS>ZGX|>?s zJDeLhagf_^EFAi^&GAvbA*{YlL2hOwd$6H`wq|wU3WWfOe}0OdN`4Qr7YRgio>&>B zip>*`!~3O|fh|1&YfcW(X|1MY5B#9AUsBLVYC9Zx%rihYucjYAi}<|oBWAeHMtG|e zikpnGarou&)U_{~s$Y%A?VY|j;!go?n^;L2o#tX+$VKSw*9X_QBuG0Mj5${}Fiv=r zaePvQP-{yoHQU%^nNzU8TZcGpkK>_7Joej?zK@WWQ>$!YJ5eo+IRtso!<`AUKZhor%C9TXUV+&8i(88C|AXWyoE0^ zf$SUG6gK1DVP?4X2>hZJO6!^%x#w@rF&Qs@GjH^-FgG;LGJ!W$u=0={F-;f+CK}i2 z@=7CqRkW9NKiCJV&W6k_r6llqe3%CLt`>d35!}j^vds46Y|-!khL8zEbMe5Ph0Lp= z1q?gqFz)?TPZCxba8Iq*fbaPaxXr5thgeNv)Y5ywcL={LQ(MWw!_)B5MT?A1dxqvC z6=9a*0M4I9FwM6&R8NDg*c z*)oCab~1gt#&aKsoaETePeJ4pxF7}fO%%}ja>wPL#w+>^ZVu|> zW!zo-G0_>eUN?iErkOPEun&FQC4(=~ibfWW6U{Hvr5h5zlGZn6(EE&GWX9c~&7Ukl z%its(GG-%w-hLbRyeY?zW(K&$YcCP3vZQmvt6{#%LsT2|CmM;juq^u<{k~um?3P-E zZ7ru+Gp``@+3}a%!1siP?6IT%Ta$3{;as@%UIQK+3gq7)k$Ci#6dv2}NaSYk#qx{5 zrPXNQqKC)G76WZKQir&4(jCH#j>QyS_o?9>!#3t6fN< z|7MZ<8jY}8W*G{08=?1iS%?}~3NO!Rf$LE%94xdz%%hMgNhj55OK|Z)f-et9(YC@} z@LblI(|RYxU1+GqJ|8(|fteInYwpKI>S%D5rhAz5Y4^D<_AdNapU!NwkL9$MZs(Fb z72xU5zj*#a7dLg>NBnxV;=0EE_1vnMR1&>jhBI$-W+IfPb31QpRdw#43{z?jbGJ)# zK)d7>4n#|FzOSM|v+52 z7_;^9YR)N61zW5Vm_y56$ya`tx#t^1;IFicjmYLgJXp0(4mS)DY{Gn+nS z&k-YKbxw2lCVI>;fu6YuxOCJ3ynpXDpPvoHMaOpFOMNvy11zR{dU`0adW3Rm&+*jw zA)+-Q@$hQwbZ{?iC$c8VRQ^OTT~#yzs=q5RubL?>v-Krf8jkppsQ|YZ74T0=m%5d8 zvjgYipw^~ea>eKXR%NFGnGyj{`JCAx&l;VV6HDKz2lC&<49waQM|R}r^8AiaviMFi z>AK@iPUl>KTNQHP@%boiO3cDVDuX}>6u9dQ$+e`*r>?7e6o^d9epvP5oPAN^;1 zn=bkHPjdQ`EwkMAG__t+P9#(CDtO25(X5y^N2lQYFHzjf ztajr4Cma>a@;LjB@8rfGSLQ#%1lXe94e5@rFn5g$ZQlL~TE{rxuE!D__VhEVU3>xN zi5Bc!H0&P-`O5il z|3fi{y_aBI`!TGXB*)G9oCM|1)?w?Q3GB5l&|MKx&^FCOXGA5aN^$v@*}x~&yS3x zO&84QluX`HyC)9D>YQcw+uGpJtx_;$P8^?^ctr0d^E+83YplwVg>2@7x-8*O>9l@>_;*sN0CoTxnzR`bWfRQY4*pW*uGZrU;V{e&ju;fDOUJ zxXe%?I)%?A(NS~B-T7}smx9E!z%>~s-U|VxM0c>6>rD%FdLg>?BC$SZjec`^FT&|} z!q)XmS+~m>_&)ateOQ}81~-i7Sz-e)W~vzuoA>}?GG9{T&S#{k&lmPzpUliYqXKD* z9eDTKd{Mvaatu6CNa8)Bz~kXG=+gW_)05}onD9o_RnCFl`8lLXpiN)%UDQV*QD9>) zpyxjfkfO-(aE6(ViJkh)1=A4t`!WL`Z~Q5WH!#B5gVyj-uz~d3>9RXl?}f`|W_Wxk zlFpJ8Vky-FFS5Ppms2ZY(qqIFht*Iu)ec6*%cIuY8kimYj3%!=3Nr#OSGIrhWzEE) zF!Xw_@cflT>}dTDEJyRMy8C6cW_b;+&FH0V2~E&!{2gSab7@)0OnTnzD~Ud=$_>qo zMBCyn9RB$|hSW(huh*Hg%LFrF@n~(X+VvNi8Dfq7>txxX7f)bE*;`D?RHxBLM$)DH zerMD7VT_k)4;!0X!TS%Z=$OVRyyg@q?)JaKzb(_a+qs`vgSOLf>krQ(9 z6$3?nnqad{6RZYG@I-wwX3DLnU%C=us>(DRRI$u<*;hyM6zn;HjyFQK+C0jh@tda6dY}!%Ocm|Zv%ZCHD7~?P1=Ju*S&@GYZ;(l zKM8F|TB6K&4Q}*z0lpa+3KcuM@%~Ot+MRtA2G%xZZ3dN~jJ|NDu$08hi*jiGxRi*e8aC zX*GwZj0{cj@G}3)f&7n%J5H`5ogaETN zY?%^30>@f#j#oQqR{3w*{W*dvOyj3TYsbNlR~g`3H(zqx=OSdy7V$aaJCNpkgKe5$ zfR&-j%;N($p+2>hn(`gy_cw15->K_}_Ku~jVQmq7Ok4uNW1nNW)*HIF_7aTU_erFY zH3afK)^KAxG`RY9V~F%<#Z`(^sM>{0!Y>|)k`3VfTi4;c%tB~*t;ig9?}eU2TR<^K zje96<2QmToP^zp4mFo?0<9e+Q z@_#p(vt5`E9?S`N#b*G$GC+TjuUD3$NP;F=i$!E(*kpDGk4R)ZxzQ)WJOr?s$Z<0_| zKEFOD6Q*w~0GGi~SXFm{CVz3I^Rv?6tF07lOEjh4o$p0+4oO4y(lGk`-dp$w5yG+= zk8n-aH5l5nfzC`Q#JJFL_@w?F?=x|Q78wBP{}OP|d3C;%a*fWcFK1;xH;CqI^8Lz9 znmGT&S=5W}q(k*HJyy_t(&^Nu;Ueal zY4Ba`PE2mzfPJPfsA|70vIDBzL%U9{f_(Cy3B0an0Ez6PuGNSqg;CP#70IrAr1d#KA^T&)vz>dCABR&0f+fHV*P{R zT>qkKvR|+fe|jZBsZbHjAJ&6p%MVe=n+y0r%K;Bh;@yX9qUgSThjDp#3Sv|>-R9cB zrJq#i>I3c({S;mLVA6fsaC;D7Ln015{tu0BO=Rkueu1&-5Y!a~!BqLHFzUsBoat*f z{C3ixYf|R>{Y_mw+rSN?TwkH}X;-ZJD#hH`I-XJcb_7%}H<4*EELy6ZAc+@`fgPjC zY)q3yx7Q~y{aFbZt- z5uf|gWM|qt{FHeKVSF)tF|(MBoZ$&>XSxv{9Do!29AU?4D_FtjNpFw8L;|)dU_-zV zT;+C`lrMQq9irl3LY+GbRVTrj%k3ZoeRyKOChjUYi*lPXh_h1^SvqVm5LJZhevri&iAVg$ui(FgTw$ zjnKh`(Q{G#mw=9ZS5G%hZ^1i%WWZ`N&#fYKMBnu$ecgMWbk914S^CTI z*E1Q$>6Is1cO+q9rzY3@iVwOYbauiT@UkVM?~9GWs3%Sij*W zd&K$yTs~Ng{xF4E(G<_me=P9o*l#FU@RrV-I-bTjI&z-(|D&Q&e1~GV8WY%-g(cw% zxCnfpM^ORBmQ2FD$>s2C`)D$4(MiOM&1AQkGyZb$!;ZDb;qB)_s0dxeUHf{SYG0IK zZ-gSvtmLROwc$=DD>4y2jr79AnQ(2l1bs5sg5Z%a{U|q*Tj9BoTs1g|r?ui?@Yhwk zuf3m!KD+`^wc}wG=MGu;0d0)V{uV0+d|aVacEg`lR2+GZuf?C@r38eA$YNiCO~6PtKy@fxFBXRX?Cm3tsv0D3#0Qs$*Lwh5blVh zG7R+{AnmbwAg42i=@o_1cF8FEQ|2=1bPXnE772p&OLgI3T8VhTw?@=nH3^IN)#5{? zR_g8ZkcHF;{PK4a`M$k}4RcvW_5b`Or|Vpp66G{7%!?!K+rz}=pU#1YYNY6C{0kV^ zehw5byVD7JBWSAM8PHzu0;hgt5Vcb!Bt_d1_at;eacK&jH_HsJX7{l&Iip~zvIX4g z6JcZd33!;51DOk(@Lv8&oNi{obvZ1i<4>fbX?7a?J!;PcSo{;J*)AfJ`_ky>)?<=D zw?*9e5BEgYy&v#^YKOwXuP~_gZMr3}}naHMY#vSd|@VqcvwCwX8c$t46dc`+j z?tKwnF&M*W1en9;mKfNk`+*qcUc%sg-EhNrg@~`HLXd?#Tw7EClSbJxYcIcn0~&F( z(kd6P$JOG2X?qx*ZHw^#>$kN1e|^i_Ni&F&V=kNX@&@jnAV=3_oDxjZmu4aqv|+*h z1rY4q0vaP9f|8~wtlnlq6{qm2rO*!OY&}a7{_!d6T{mf3wI=cV(IOQnAnf?1mG9K>aPYHi4IjV;r$83eS4fP~o zQy?&kYhbUMI+X74!BIBWO#F&qTyr#${i(YS7r&~eFYWyee=HN8k13XXo9M&Nd(b0ln7WhIEsheF zZ8T!TvZe?u>f(iiE96BBKUIkFm>MgNZI(PDLx{oiY$bm)j7>Gjc4QbP!#1O~gZLOzfrn^u*tnbUK99Ubnwk5F%KVKb|V) zC_218>m>f^Z&G>wn6tRP_*>=mQ}zPlD-h`ItE`+7cvv7Dd@7Jxv(Dj5NwIyzYYhjN z*9R&y0#zltlV6D}bIK*deykJMTB+Cx{JZSr{|y!2f0-`+XM0F^`tE)Q_4zj|??2ZQ zdx@`CI&IVx*E|dn#F| z60nu-4rzjO3*E)0{ucy6qEGhgdX`kwn#l;d9n}P13*QOO<{c1hoHs<^p1-T&>(K;p zVa`wzc)U{-@ia!*G!Vji9Nogo|1~F{e@hd#Rhd6OLPf%Td)bvAq{QnM=CDOiX41N4 zse(f@&IvQ`NwY=Y>MEb~3?>0qbbD3#R2JFq7oNzq?xXYF?c5A(?xFX zAU}8%Qd;N2pSlF785v1F8>N!(N5`Yd<0xuB`v>*xQ-z2TxuVRYNyO|)07-IC;VgK6 z7B9U);}Qoj_{4L^Mm`dykFS7*N=9_qrgjX{nT|^)bd!~(nlQ7Y3VzDnhKBA&+@0=( zF0-T=^NaFW@JJo?oGr+*nfGC4p%Cgmo`!;JmL%6RkaUGNkc-!XsiS`_DUCWKDhlH8 zt%WOE-n67Enpa}OMK{*>SQr?DwTb3x7Qx!`v6ZRQodg2|$Hk`?oEOYl8^?p2qy#5x zW{AQMWQcdTU9K2b>g7@Att`fg0kFl%OeP1b`9_LVDTwD39%~5i|J4F0UCbDAp zxs#QnkG-y(?&n=OZ+C>)rZQNvVY(rk{rj+hiT4nj@%ZSK!Ri8Wv4TLiYlP&yorlCT z)>dr%JX|cj)kx485+GJ~+U3x5SV=JK+(d^_haQP1Sk4juE&gOb#7b3s%*(p6N_V+r zmh;2PXFtUcN65Y;tIQn$9&A@|x>LGg4wo599pV;?eASYB3Kzcmo97YbMT(l;Sj)iILvL|B{2I^Sh?OKq;g+s zxcJ}j@XC;?W^u*a727K(kXoy6#{a<>M)B7qbf5Q!T9m5eZRv?%d9RaB*S$x} z^ha>_txMqT(p9K8=?E^dF=soJ(#Sr~0vKHD!*<*FgZ{2w+Izv4e)v2T+C8Rl@8@o& zUHtuaw}jtqnfgHV;TTwDGzGIllz@&tK`I4b>86USXu5w1_P6tl!rOyH-T&fojG&nx~}v2yg${J zr0&fW{`Bi2kM)ZM+~Zl0({5q8{v|A1y4pwhQBG93#QV{P??x!)ypzoSXWq2?egol@ z>QkIa*EBx*QqA^1bz%=|-(swdi&>eMF6^Z}KJ@mWGZa`QPcLZRLzTj5tmc(O?ncCP zu32`z@QyFI@DT)!s`^7A6FC^Z^Or+syz}J zsLA5p-6h}w+{wfYRG^*@A@JTuCFB=S1#&K5CjKwF<J*;+3gfF*L!s{J+tWDD- zvT6JxY#WOKYA-KGtk@#)t=xWOJ^;yew2C0T^=>=N`% zryCutEhIzihsg&yCxOokYqa#&8M1s*De>t|H2;dV392>~C->^l22+ePh=lv=K)%*d zXy5J)wr&nUIfXNU-RTWre{d?Apo55UCV>#72_awaJX(>R4*Gtc%G_m4A+*%f$yeJ&f0GL0hLztM{^v4>L}vJ8q3Dj5&&+MmLvbL z4~fmf*U+R`4;ZT0p&17TVDADF8R}axyB75k=X%Nr%aSqZ^K}p(49y1dUo@%WHMh|= z^O^YGSQ?VgK0!R_{Nx~gVjFMei*#UUV+)JlE=GNGeL%|R9nezBfz(}}jz-$DVcYg{ zH2I$sIp=B@-@s!YXrJ{K{IIYAliuWmy0B{WcB>59EG%L*dJlpwS_c{5GQ(O?Q<1UR zD5^h2K$m1qWS`lO{+q-H>EF_Ut4lQ~ntqfZijDBV-v+2B+AC^b7cl2{7VsY=vdD8G z4@Czw2-s&4@NwOKaHGZ+b*nv~pd!PKNeVkl>j36KUN(Mg<93l&8PrR+{6Ih;n54>n<1=|8r z$?Vs&@WhWo^2P8;#{1LgcefNzBRfau0luxZaFyu9TxD9i<9?GknJZ_-)(s$?N+xqAjVHLDJ; zkCwpmDP_T~Rz+mBy8_0n*JEG!^`UapPAKss8ks%x!57EokRGmC$U~zNyt!F~j0^XKJ9FBk#crIVb)|-$I$o8YLYEeM|;ajS zDoyaz$5>b${f zq=#F`5MoQ@FJ!Ph80^2S55BEFf%jXEqb%#6%(m@0g!@Doc=p`~oOsel4CodKN*`f# zae61e!07}rBdZPSoGwJ2TaTlz(L~T7Oebw?FFWj#{SJSTy(o09433I5hWBdtpnmoz zl$yr}kN(ak3+}Hb)i&zz-!8}oax>)6gM}Ah!;usa;8kJowvt2dUfye*twkVysuEh7 zCx!e=5`l3;3sHY`GaqexMs{_o@XRGck=tYe8Majc3F4NM|5@AyC8a;%^KL6NJunwO zP;&&1UnU6a*k0&z&`n^yYi(WS@fcus@f0X8=J2zzAfmG3k|@J)ka(Myg8#Fe%qD!T z06GRW==|bK4z8a|VA&U4pg*1ul;?TzP0j5{r#*2*>b~bdX*f<$ISycBqc1FdEXBK1 zdjfrpkRf)q27qi`Rlu&&MJIJmqaEThyr1f`$b=X5On0pulX5u#=AIb@KL_>_?;Bk~ zk-rgfWwQ(!X;n`QFHs>&1|(sJU^Ot(j76Q#RdCXsb9hwh46`-f8GartK?m}Wqr30g zp!KCjMrd#kHhHgs@-1l~Q5OJTkL75JPX|(om`0kTK)AK3Td?7+1S$9LDyn*T3eYX% z1kCXyeAj+Q2id7)SEL``*n5wl`OG@#`Y!;)?e0Qlc1)gF-Y5gQVzi0t%5I2LzE33httD#mEP3mnXtPq* zE)3)`z?{$lI~0x!Yr!jcB%u<{`avMxwmeeZNuF%HQvfy<)dG9g3)RO;@@c zX|TR{DKI(E0^(gq8z+sP2FG7r;_bP1flTrJ0ft+akf}+>z|{pCU_G+}Z8#H*_@Nk9 z9&K?@9X-OI`S}MiIVYZQv46`eez6QKQE3EEH|pcQ5CNDd`wI8|QUbz{EYiHN1r2Qf z&HGc33vY;7g7_=;q>Gm?|K)HKK7Jve=n&ol8)tFEw%56E>#kz(_^tuE@w^C?&0Z?{ zk2?@eo6V6FFAP*}J&68IDg-*7-e_#}G}=Q1gH1Ju$@w=4hd=i>^2OBah`FyD$z_V( zpkksGCEXVnga|5uw^<&Nds7LzBHEF1`6FoT#yaeL;{l$OhM_NFJoxF;Ah=( zMVuEh3f}*+2D@9!z**OoU|53_tXA?v!bsl!;-yrM@B>Iu(& z_$bn@KjUxw9nE3$5>8V~pUcbD=B(bl!lzR+X|UuVC4Xvy^ZamyvMG8%(Z-o4DO}bVaS0Q0`AF=ettWr&UD=$+bb;`0^7puGe!sD`FIWe)f@`9coBv z)U+@PcZ^9zt#EE>eIYsXqbhDHI?g5gTQW=h*HIo@W(c{*uer>xs~lH+7xLUDyP|~I z4?vYz49phSM>2a)qJtv#fOWuJGQr6eEKL#RpidnFV^ikfY0>uBjK%hytHQyC|As-l z`e7g&ZwRu~`;dRfbTGZN39j$9!;@{xnH?93_<4bLDC|K5TIQ+&_GHPT{KePMCB;YN zicjjiG#LY&vS~fI(jkjFHdza98H&y(>HyQ{yADUbO9C6JEkOU-B}|lVKG^ZHnmKG^ zgFGVm*lr08THehdA59d%EbDyKAn+h#56lPi6atXK^JU};Stqow|0XFhwHOu6nhPag z+K>kle8~$EuGA~|3!LP~6SLo`p$Vx+u-S+Nu_LKaX`dk*H_j)B$DU}>?nCIXmJoHF zHb+J}5vXBLJBU6}f$W^Ep~1&sLT_RYDvY@V??|)7{VX3Pmv=Q-SFtUgJ9k98iCiXe9#n-g(z4FZx&2PJ62VL zYu0-}f7vZ4p6tmpBjb^*dmxrc^%waWGRY$Ehk_Cr9jKoEPLQ{EH!}KL11b)N5`Ldg zv--w>%>7RboF04!o-aCyZvJ=(76s~%6D|wMU)LXj-oRAQ+0}?1tQ{a)J{F=m?nglP z-L)WD=fB1b@o0d9tH7J!V_;e0Jaqe*ABghPfL2ckxcZMP>8JmWs9rKe{CvESd_8p@ zaamK6Ja_jSe0Fg*Gi|jvyF>mmFcLe7RERBDQff8v?&DJQXthe?^wBONX0j5&Z=4Qq zMtv5@hx&p7)u)1VeGD2x`M5@1k8CbY1Nu!L!KO?Coi9p+9^O*?SleHM`1>R^{rLfxwzFg_8lx<;*1VN9wWBi??Cw~b!g5fZ_-6N znTVs@VV`ytIMct8j4V0MU$lEKweg%CZ5TSvx`a5;O(tU8x4TiCho>05ZPh!vYvLt4 zMN*1$NXcQH6ML!bf)s9N&l_%Ico{z-08{^3WrY*7kJ3vf)lpHK?{L$!rgFY>2e>!2 z)r>J0MZ3iZ^L0hsitl!txWMN!v>>^eyUM5tx85J1M$XOWu)eE^|CG%+Jp9eN6(|VB z--~j?e61Yaueb`G9ZI=m)9v(BsdQS?1(QP#63o9eL$0mZSg2uRg`@xN;(GUJa9_qY zQAv|m@Mo)(QG?W6)?bO^dsggcOU;bvuWC}<@QSU%TTxD||CcISOKt%C?S6}z)C=*_ zYYeh-(;#y_{YX=1E4Da$GC!vM5IA=(9>mghn2P#^=Bz!+WV@=6t;QQbs`r3Ho~Rz! z(KDA!oE(UitmVnh#5z>HNgJ1)muBv?U0|+y-hj`RClNgO7R(IzfIf6b!y~*+VBlaH zx$N#_%=r0|(SrGSvC4d86cLI?-}f+clqa*Fd{fC<6H$NRKa;t8)SkTj;w1jQ@EWem z9!323&B)k7nv}fEgMZFtW6vUEv|PFieGpxWu8cP#gQWe)wX=5c8CrpCQtkrO4P7$f z-3Yn#K_=BwBTaA5%%$C4^wY1miP4gknsnK+>$G$oqIxI9Xv;xu>Y=qR=eKGe9diB) zztdHX^F8g(n$(N*i`$jyAX18ZQ@WBCOhj}3GZclV&Ml|ErxsJ+mMs$g{FKC%1o7yd zirwsv{Zh1jxes?VC64SopUc^Iey1<;O}STDH))OI*J&?_G_G0Y09QJ15`FFca_+?4 zMAm0Y5mljfhboFwpbu`ca$Kx9MsWA$bIU|pVoe`)?&z6++@Re#IwaGZ(K%Pkef!=b zVvHI%epr5oPC0dq>(VG^Uq|xjnZ7!lO8)}7MwBD^;Fu9Iw~|2@x~GG{0U>X`mKWK7 zXfDL(_E4|Rdtf))B5>eqD6k!TB5-y(2TC`{l5N>11#^THv82cVMcnR#?ddM4UTqTS z>#8Dl`KXbvmuTRJfr{w9-))?;_a%rprV36h80Iba$Ixuh1xeLrkz2N1C3g74fNMd@ z;A;noZd96s;b8}~u|kGfF`mW9HF=_jfW^o`rweQUItyofQ^IQ^G;lXx3wPvpqDN0v+R1`lBS_DUd_3FDkUGq1@pI5p_`r>r3oSX^ffF|N@%`x(+_71#u(`~-= zEjg5$KTY&4972-;g-3(UkgM_$bhAQ@MWab%g~10ww=x4PPXB}P*9(c^r_HF$bQqd- z$pU)L8Dw|*He7nC8nkokd1rg3;7R8wxk2_#-KgE3AzJw?P2u#jBDA2Tt6z$_}=!3;Ip7Z|Wjnfn}(X17{$T(+Aa4;@` ztmJ(ILsIo%_Vx~93|D~govPH`r-i6y?Mcw-Qv^QTF9P(`gJ?Wq6Hv?ap&o~iAl}ls zAkVN9t@uKsT@8k$O4wRb4%|m_LrcKEiWp?;sf%7}2l3onyNHeU(uBNME3p3Ui+B}M zRDJIe5`26KmW~9VU)OFz!}b8d*ZQl13kH68s>EDyWAZxk+I=lln$tiOy=ou=4w#YO z2CVqj#n#X_w*p8sXcIefx07f3>QR)R4CS}VkTjmR1W%1!gX3H)1a1jmkZ`9aSiW=x zZ*9yO!PG|rNCBec%g;ME+LZy37md% z2Z-+u!^Nlnq9f)f+25)5xI||SRQc2)co`Ij%M4|~mE*B^-j-wNufsYpl=THnjXz0V zeYk=2H7x|w<_3VShGs7j$JCz@4q~pl6jBj+Oy1%%zSj<$Te% zS#r!QC0P_Yb1ph1n2P7q!KkKUD>H1gjB>LV1L99UgRzsvjJu~2sF5{co)}KTQs+xR z(cJ~$(MwZo+kApFJi^1_1qZ>x>UGF1eLebX{144gE&#pD4ES4W{^ND0-Un*(d2q|r z?;xG^po*XP!IwWbfzMM1f%W96pgm6&v{FxyRr)ZQF!2%?XWRnE0`rjhG&{tF8xdjC z3do>XKak)iNm#YqLIpbrJZW-1bhXPxYN}RvTd+2HKy;sSAj**ceeHR^m%#?OLv=mD znQetf`xC+7+Hi19LJj1br6Dvn|d3}yK0Ahx;lx+-E z-G(db+OW%ya-^GLiLN^^g4S#5cumM5Cg@TP_a~%Tl&5YYr1EmOrJsLr+O_qZ(Gw|Q zvA+~|xj2W)pWI02#yN93uddO%tYoQC+UUA&w3bWzwkT zBy+_DS=6cnMvgDOP;5xsHIUr-gl-A9V{LCL)B4^bC-tjoj+6Wzav*|7T~`|BHYbg9 zedacf>71IQq4Z|B*23(}Q`FzXmh_kD-|4Qvz3ipA z4vtbgDWS%faTN1K$FVsO2n+Hpgu#}Rgt2~$9nG*2Th0Y=x&PTa4$f)fG-FQ*(>zCc zyE9(E2c$D-!yovbt(}r-U4CK;SR@Ge2ZpzP zBMetn3cAjWpzK@a;KNo^K6mpNsJ&_d`WE*(+>$*4B*{HUd14`IB!|d-QYuJ)35Bw3 z*PvV1Pa>IPGsx&=CqaFJ2b74DLzy4lVVToaG*?R>{o1z}#E(x#M=gHCr-l&t|IlZC zTQ26$`BFyosGmYT|9PXYod7IwwFFx~y5h6#n?P%l1^(pMh(hQ4!ZRDn9PFy|1$W${ z_+MHysCWK7@S_Y6nBE*CzZc(YxOmxy`NCHu;v&wI?7>OsM7g~nV;F-u_F|-T%N1n& zR25mdlStOd4CZWXhDwKAvDQgBq*UgPj-;O=yzX|Q@D?RDEHR&d#pyElt1^+xtPW&l z27l7e6?)hR+KJvge?8atWvPh0{FcuZm^j`YRuIDd(`mkTJpINdoqj7T;55eKs7FVw z>9Nf-=zL8l`r;!E?$WB;w8T?4e)*+e?1-fqhtXcHtWTPI|23YSSXa%A2S(Fc509|} zsm1Kf>618<5HZfx>l;5NIh!j8TT6MKvS-}HPEa-(JbG|kRp^$wlCED9$6sUfm>LS2 zpmh%~r)%5fX#Ut1O4I5WWw-bc#k=22zdqE!R!Sxc#ZIVm?}jGWg9@S?#6`2Y=><8Q z`sh5GXz8Kj15)6%VI7d|WeX-2%Ypa1R0MjnXms9FmCAg05GU-vN^E$Uh|03cki|n? z=|2=Ln=i!5Zc$lNJEz>znO$T;~jT+tpiYWRtx7Sj7O2^{@Ugx+iFfLlfP1Op$mft;TQ z*x>Ys2#(#3`i&ikNZC(Ba!?i&CB2~RMahiQoyllAS1LFWoP@3o_X++e-9i#S{=jF^ zr9k?kpS6~xT3%$hOpk&iBvTko9Il3@_?D05+yhd`sov;7UL7aysNw)$O`9W}8 z<_!4Kk%%WBNfwN zaL-?r^o(#r!#k^h>|P0=K7{}sPD-R((r;8^eI7nOxx2x@#*my8mx4F#C;=ZTrSPX< z5AgkaBS<~q10Em9=jn~ai0YP|{6E4_Y zbLJW1M$V(~vq`+%Glm%S7y;VK7TL}ZKpU^}_{n7l`44%f@NuywSa9z=7z@usi@pBw zjeocCEF?+1=J#!s(6JbO&Nu-6?fD9Rc&9W@T_1`h)B-?Eb~1#whM<0D6WTDz6=XZ_ zBeXZzV3&~{puIztJ;u}`wnd((XqJML-(Dd6^~?ydt4E2);(Yui`xlxy>JPep-Guwh zz3?#ynsiwvLk?~~{HU<`Ul?iq?eQf^aO_1hyv!dZ1PP(mHc{XC7dRUZp1Cjy9HsotpCJQGAd zngr@p)`O69>%hi{#W+W7F>rdiituzwMrIiYkud%QYd|S2(0}U`WWQ(w+JlmK3vNe&$bIGL$g&|InHCFwCzKHnkDtKfQTtG{ z&oZ_;b{e#I3By+Rn^D36h{RTTIjA>fA)UZdVt&9&;)k3G_^W#WmAUuBfl0*;y@Mu5 zuUZ|r7K@SR2UlV3m4W2m-~E*8Z$s|Nf-h{r+A_9Kv6xECJiy@xzxb0Ak~pW$T5Q*e z2=;MB71yVKl{v}5tBqoQNT3|}O(RAZh4Os|Xzcs%m(| zxdj~^av~fr{&L`6Y-FviW5~g%Xpvq`0nG1zf!=SILv!|hL0W>nXyj4_kt}VF+$0J? zon$W0&_oNl=Y2+{8A7NSluiciso~p<Z5s>C5#J?6FQpWxY$gL z5MP!~zWcon>?sfuebjAa@G%35En3CIPBLW;x?^EfK`0S&q7VRy5a3EZ;VZb*pu{EZ zU`n47XkW1#6|P!A>?~VBxVK$HUNH*jQs@$V*KiJ&ed~$NYpx(er<=l!#+L}4J!2qV zzJ%->r#uyaq#=j-a>v1Te-%*UTSBGg z>7ecLYj{?b{e;uQAmAHmicqr$(x0b;7Vi@yO%Lx!|BYB6=Sn4_hKYkSud9&DGB<+x zjv{XTLndtxEbBdR? z?kU{bpn?nECNU|^r@-Az4Kz=4IZij0f^PGdFy5CJLL3-?jkZ)E&%e$f@w_Pxzat4! zzqi4va|NKn*b2=OCy;NE4ytTu1i9wj0=uY)Mx*`fz{kJ`gn?8yDpAa$0#2)PLeEv) zDd}7`?0E|Hvp9vjyV{ga)$`{v#lNsahn(5X_sh819TA*`1;JL`zrZn61~q(BmRtHs zNZC8pbBn!A=;5`;XcwUk<2Qbu^-Y>Q#Cn@YZ!mW<{C-RVwNWSo!uwc z-JNDaHLr4d&%icLE%GpT(Z`4NU2%=xyGxlqdt04e+qjl?3$o_k%}wHz_WM&Mu|OC$ z$%8w+;0vYTww*4w$)_)N2m5OBF6LHTs-GTBs;z<24S^p ztVOf3(9ph`-uGSCaplKIa{5>ioZ%#asV)1V!V!G|w|FnAzBiXmjR?b&`*)-Ao~uZ4 zx`605Uk0kR?h&egDe|||6g=Bko6PYy!NDIr;c3O^;G~5vM$*q9KjInrs7Jsq`=-TY za<|dOhcQ6rjYy|5-vV^XzXMJ=OUc;rGH_z(G71lu#3l7NQASh)q3oBB*B77XgGJ@| zxOfqs8>j&J=>-4qBR&8d7&fO^hV1{knH-rkD##dmiAuL`BhAuE$(W|YIN+HZlu4Au zb%F=zk8?MOUlxLID7u4@%U2z$L(PHXk%!>WC380X%GzsiltO*FmFd%+(GL}(} z!O;q7R0L(k2#0@xG-Wp+(Q*>(vAPa6U9sgW|3n~gMFFuA&j*Sy39IB!!Qs(LyicSBW1wUXfskf zM4=hqXOfaT3V~?w!C&>*oSOKdi*5I%L5;khWQ*-wOs9S0^U`mS;`%#CEn_WgBKHTK zkK6>d>9hbPqn)5uRu@$Fz6EuMRpFwKOA4S$M(f zVxd9@LCN+nVYEc$>-|lT&dgXyzrA~h*t_O2FdocAKNC-a8l`TcnK1_qdgF|v%?)yC z28sNFQv?Aa`A{8{;MzKU{8{iFdCFrj-OUB)pu_sQ7+AFuGvj{p_-6*8{Ap=sPSH{n^>B`8RzCr@%U1E$%vV6xtCoPG zUTybP=Gq=oc2T8_p}t;L5_i$S`GJKH}080s8I)#`iIpgb~&6$ z+pU?+?F>xde7+v%1R4W$YQp3k{M({w4BGXiQt`!^Lc zJ%iqOTg0*u2)GZ6WI4|(!(4T#kh#;I&&Ag~b{W&H)}+m%681+oFp3A*dJ0t?|c*CqA@ix=JyiezL)}T;fpx>;`BtitMLjMmR&;W z-Z$l*?iVqLvzQLUM1wb$qu?_v=MWzmk^V7-6G!p zy9~Dp8hHyOo57anwIHzT4|@MT8Q&KV2b!mDK~R;3PHa3uD7~*>o|{YuZ#U0`$0k;Q zWrj&0=meO zG5Kf16qL661QGZGoA8Cstj7t9lGjLGe`ZD#9BaJNq2$ zf2#50_#wZX(%U*NFkGt%raWV5#ZL&|BPL_vrX@%rBb z80nXTUW(m=JD&@{#uuHu4H+xRas6?E{3!`SzMdtV10%_ejyLc<^XH80ifW?mE`@fTUwdxHA1bWqDwqw7~v!Go$1Xxg-fsLj6t-;rXxuA@@m zhv|9rEX9!^Pxc6I{=175XBvh$a#96SJX&TpP%_HVwV@j|RFHjsOOg2-H+2h(9yV zp|MrZnWN2sKd)Q|oikX7&QI8Z9o?=CqRRqC5d(cdj6=hBv^rz<=;&+C6C8 zo&}p$^#i;3>!ebm6;L#ZH902>6^SD+BNt5+&VADD$ zyz{j!b|}op=g#Jn@BVeLip2_yw(ual^gt2u8dmZbMb|OybO=8@k;F?r97GY%EJ^p$ za&o-c4pl65!uz)-l8<*ZBRTg|{0LeDyNE^!13gmA@wC$-jqw~ZGT$0A#$)6>&kN-0 z4Hh`(Y!*&T*8yEU4@qW9o?ut>V*GeWgL*s~!8$d(L@!MYu%5Uc>w@cXK*12UmOsp| zp0|$eF4t#@4DO*m-F9#!?;3Ibq#{!@vKK#=Ge--r1TiM}whE^B7>j%zYjA)~0k~sj z!^r96kzVspvt!zRHOu0d> z(y2odX)GA+_ysE4((sg0T{3^gY{o(KQP z(M9a_HcR$NT_@9HeS&|ssGocGB%Oa`gE&>V)C@&tFszqy6?b{deD*;6ero!Gm+YoZ z9~s_9$b^Xqm734aF*zGEu+Fxd%<*jrlzD$A9YXnIeo_G$HnD;--4%x;ocYu>P{WkI zW#|t^dzs#{DCz(W*_SDK+==BB6=tK!#8q{1dn@z!6LC+e%5gbI1z5tqN^Hl?ijve9 zhdW@s?@DTBwuGFx9A8!3Noz~^wvhEZ-L`KG*{{SAY?YuS!f>hyoBJsanZ4>W%LAj#wT{A*Yl{Jv37fo)7TI3mc~ znD5}wx5_~xW>}yodAKpP%NyMmLHS-4nGtfQXE(CClX?AHJR4r+Z)NWkh~tjT!)!`H zAUK?>3JcCfIlPCNAjRoF``U{~80?U1-?w_P!`{Nhc*-7Y2Z;t-@LzkPeO7x({f2sd z^mL%d;mE1gbZ6}=fiqRxIBC}u#tky_L%7Z>ltw;300N<2j%}h4*&oF literal 0 HcmV?d00001 diff --git a/idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST_idaes_info.json b/idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST_idaes_info.json new file mode 100644 index 0000000000..d2970958d1 --- /dev/null +++ b/idaes/core/surrogate/tests/data/onnx_models/net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST_idaes_info.json @@ -0,0 +1 @@ +{"input_scaler": {"expected_columns": ["feed_pH", "pressure_bar_feed", "Na", "Cl", "Ca", "Mg", "HCO3", "SO4", "K", "Sr", "Ba", "HCl"], "offset": {"feed_pH": 4.0, "pressure_bar_feed": 1.0, "Na": 0.0, "Cl": 0.0, "Ca": 0.0, "Mg": 0.0, "HCO3": 0.0, "SO4": 0.0, "K": 0.0, "Sr": 0.0, "Ba": 0.0, "HCl": 0.0}, "factor": {"feed_pH": 8.0, "pressure_bar_feed": 405.32, "Na": 135.99, "Cl": 180.0, "Ca": 10.0, "Mg": 10.0, "HCO3": 10.0, "SO4": 100.0, "K": 40.0, "Sr": 10.0, "Ba": 0.1, "HCl": 2000.0}}, "input_bounds": {"feed_pH": [4.0, 12.0], "pressure_bar_feed": [1.0, 406.32], "Na": [0.0, 135.99], "Cl": [0.0, 180.0], "Ca": [0.0, 10.0], "Mg": [0.0, 10.0], "HCO3": [0.0, 10.0], "SO4": [0.0, 100.0], "K": [0.0, 40.0], "Sr": [0.0, 10.0], "Ba": [0.0, 0.1], "HCl": [0.0, 2000.0]}, "output_scaler": {"expected_columns": ["Calcite_ST"], "offset": {"Calcite_ST": 0.0}, "factor": {"Calcite_ST": 98.19}}, "input_labels": ["feed_pH", "pressure_bar_feed", "Na", "Cl", "Ca", "Mg", "HCO3", "SO4", "K", "Sr", "Ba", "HCl"], "output_labels": ["Calcite_ST"]} \ No newline at end of file diff --git a/idaes/core/surrogate/tests/test_onnx_surrogate.py b/idaes/core/surrogate/tests/test_onnx_surrogate.py new file mode 100644 index 0000000000..623254c285 --- /dev/null +++ b/idaes/core/surrogate/tests/test_onnx_surrogate.py @@ -0,0 +1,173 @@ +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2023 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# +""" +Tests for ONNXSurrogate +""" +import pytest + +pytest.importorskip("onnx", reason="onnx not available") +pytest.importorskip("omlt", reason="omlt not available") + +import os.path +from pyomo.common.fileutils import this_file_dir +from pyomo.common.tempfiles import TempfileManager +from pyomo.environ import ( + ConcreteModel, + Var, + SolverFactory, + assert_optimal_termination, + value, +) +from idaes.core.surrogate.onnx_surrogate import ONNXSurrogate +from idaes.core.surrogate.surrogate_block import SurrogateBlock +from idaes.core.surrogate.sampling.scaling import OffsetScaler +import json +import onnx + +# onnx, onnx_available = attempt_import("onnx") +rtol = 1e-4 +atol = 1e-4 + + +def load_onnx_model_data( + name="net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST", +): + onnx_folder_name = os.path.join(this_file_dir(), "data", "onnx_models") + onnx_model = onnx.load(os.path.join(onnx_folder_name, "{}.onnx".format(name))) + with open(os.path.join(onnx_folder_name, "{}_idaes_info.json".format(name))) as fd: + scaler_info = json.load(fd) + + test_inputs = { + "vars": [ + "feed_pH", + "pressure_bar_feed", + "Na", + "Cl", + "Ca", + "Mg", + "HCO3", + "SO4", + "K", + "Sr", + "Ba", + "HCl", + ], + "feed_pH": 9.5, + "pressure_bar_feed": 1.01325, + "Na": 0.230858556000748, + "Cl": 0.106701648328453, + "Ca": 0.0245274696499109, + "Mg": 0.0311348703689873, + "HCO3": 0.430482141673564, + "SO4": 0.182204065, + "K": 0.000500561, + "Sr": 0.000761853, + "Ba": 2.50e-05, + "HCl": 10, + } + test_outputs = {"vars": ["Calcite_ST"], "Calcite_ST": 40.4270772546529} + return onnx_model, scaler_info, test_inputs, test_outputs + + +@pytest.mark.unit +@pytest.mark.skipif(not SolverFactory("ipopt").available(False), reason="no Ipopt") +def test_onnx_surrogate_manual_creation(): + ### + # Test 1->2 sigmoid + ### + onnx_model, scaler_info, test_inputs, test_outputs = load_onnx_model_data() + input_scaler = None + for key, items in scaler_info.items(): + print(key, items) + if scaler_info["input_scaler"] is not None: + input_scaler = OffsetScaler.from_dict(scaler_info["input_scaler"]) + + output_scaler = None + if scaler_info["output_scaler"] is not None: + output_scaler = OffsetScaler.from_dict(scaler_info["output_scaler"]) + onnx_surrogate = ONNXSurrogate( + onnx_model, + input_labels=scaler_info["input_labels"], + output_labels=scaler_info["output_labels"], + input_bounds=scaler_info["input_bounds"], + input_scaler=input_scaler, + output_scaler=output_scaler, + ) + + m = ConcreteModel() + + output_vars = ["Calcite_ST"] + m.inputs = Var(test_inputs["vars"]) + + m.outputs = Var(test_outputs["vars"]) + m.surrogate = SurrogateBlock() + m.surrogate.build_model( + surrogate_object=onnx_surrogate, + input_vars=[m.inputs[input_var] for input_var in test_inputs["vars"]], + output_vars=[m.outputs[output_var] for output_var in test_outputs["vars"]], + formulation=ONNXSurrogate.Formulation.REDUCED_SPACE, + ) + for key in test_inputs["vars"]: + m.inputs[key].fix(test_inputs[key]) + + solver = SolverFactory("ipopt") + status = solver.solve(m, tee=True) + assert_optimal_termination(status) + assert pytest.approx(test_outputs["Calcite_ST"], rel=1e-3) == value( + m.outputs["Calcite_ST"] + ) + + +@pytest.mark.unit +@pytest.mark.skipif(not SolverFactory("ipopt").available(False), reason="no Ipopt") +def test_onnx_surrogate_load_and_save_from_file(): + ### + # Test 1->2 sigmoid + ### + _, _, test_inputs, test_outputs = load_onnx_model_data() + + onnx_surrogate = ONNXSurrogate.load_onnx_model( + onnx_model_location=os.path.join(this_file_dir(), "data", "onnx_models"), + model_name="net_st_net_5000_STM_100_s_2000000_60_5_tanh_1e-06_4096_tr_15481_Calcite_ST", + ) + with TempfileManager.new_context() as tf: + dname = tf.mkdtemp() + onnx_surrogate.save_to_folder(dname, "temp_model") + + loaded_onnx_surrogate = ONNXSurrogate.load_onnx_model( + onnx_model_location=dname, + model_name="temp_model", + ) + assert not os.path.isdir(dname) + m = ConcreteModel() + + output_vars = ["Calcite_ST"] + m.inputs = Var(test_inputs["vars"]) + + m.outputs = Var(test_outputs["vars"]) + m.surrogate = SurrogateBlock() + m.surrogate.build_model( + surrogate_object=loaded_onnx_surrogate, + input_vars=[m.inputs[input_var] for input_var in test_inputs["vars"]], + output_vars=[m.outputs[output_var] for output_var in test_outputs["vars"]], + formulation=ONNXSurrogate.Formulation.REDUCED_SPACE, + ) + for key in test_inputs["vars"]: + m.inputs[key].fix(test_inputs[key]) + + solver = SolverFactory("ipopt") + status = solver.solve(m, tee=True) + assert_optimal_termination(status) + assert pytest.approx(test_outputs["Calcite_ST"], rel=1e-3) == value( + m.outputs["Calcite_ST"] + ) diff --git a/idaes/core/util/utility_minimization.py b/idaes/core/util/utility_minimization.py index 8ea9e35bd8..34ce574b27 100644 --- a/idaes/core/util/utility_minimization.py +++ b/idaes/core/util/utility_minimization.py @@ -28,7 +28,7 @@ __author__ = "Alejandro Garciadiego, Alexander Dowling" -# Temperature Epsilon to add or subsract 1 degree to avoid division over +# Temperature Epsilon to add or subtract 1 degree to avoid division over # 0 in equipment not changing temperature EpsT = 1 @@ -301,7 +301,7 @@ def heat_data(blk, heating, cooling, DG_units=pyunits.Mwatt): Q[i] = value(pyunits.convert(v.heat_duty[0], to_units=DG_units)) FCp_[i] = Q[i] / (T_out[i] - T_in[i]) - # Generate a large dictioary containing all the data obtained + # Generate a large dictionary containing all the data obtained # from the equipment exchangeData = {} for i in pinch_streamsdict: @@ -413,7 +413,7 @@ def pinch_calc(heating, cooling, exchangeData, DTmin, eps): initQw = -sum(value(exchangeData[i]["Q"]) for i in exchangerdict) + initQs initQw = max([initQw, 0.0]) - # Fill Class with all the data to initialioze Duran-Grossmann variables + # Fill Class with all the data to initialize Duran-Grossmann variables PD = PinchDataClass(initQs, initQw) PD.initQAh = initQAh PD.initQAc = initQAc diff --git a/setup.py b/setup.py index 407a027159..d5303b756a 100644 --- a/setup.py +++ b/setup.py @@ -44,6 +44,7 @@ class ExtraDependencies: omlt = [ "omlt==1.1", # fix the version for now as package evolves "tensorflow", + "onnx", ] grid = [ "gridx-prescient>=2.2.1", # idaes.tests.prescient @@ -124,6 +125,7 @@ def __getitem__(self, key): "*.trc", "*.nl", "*.keras", # idaes/core/surrogate/tests/data/keras_models + "*.onnx", ] }, include_package_data=True, From 3e403b78e3aaf5d6c69b1995c05f859548d5a097 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Tue, 5 Nov 2024 12:23:07 -0500 Subject: [PATCH 4/4] Diagnostics tool for ill-posed constraints (#1454) * Working on cancellation detection * Fixing typo * Improving efficiency * More efficiency and testing * Reducing duplicated code in tests * Finishing testing * Adding additioanl tests and clean up * Running pylint * Adding constraint walker to diagnostics toolbox * Addressing simple review comments * Accumulate node rather than str(node) * Clean up some tests * Supporting ranged expressions * Apply suggestions from code review * Regorganising and fixing pylint issues * Fixing false positives for linking constraints * Fixing unused variable warning * Fixing bug related to external functions with string arguments * Addressing comments * removing debugging prints * Addressing more comments * Debugging issue with nested external functions * Adding todo * Relax test tolerance due to Windows failures * Adding catches for string arguments/values * Test for nested external functions * Fixing bug in walker logic for string arguments to external funcitons * Address double counting of constraints * Fixing typo * Bumping Pyomo version * Adding zero tolerance * Adding test for zero tolerance * Fixing typo * Logic to skip scaled equality constraints * Handling mole fraction type constraints * Catching negations * Fixing typos * Limit of cancelling combinations * Fixing typos * Displaying cancelling terms to user * More detail of mismatches * Adding hooks and warnings about new config arguments * Expression aware to string method * Moving expression writer to util/misc * Typos and comments * Using StrEnum --- idaes/core/scaling/custom_scaler_base.py | 4 +- idaes/core/util/misc.py | 66 + idaes/core/util/model_diagnostics.py | 791 +++++++++- idaes/core/util/tests/test_misc.py | 56 +- .../core/util/tests/test_model_diagnostics.py | 1377 ++++++++++++++++- setup.py | 2 +- 6 files changed, 2256 insertions(+), 40 deletions(-) diff --git a/idaes/core/scaling/custom_scaler_base.py b/idaes/core/scaling/custom_scaler_base.py index 8e94cc7c5a..a6b311086d 100644 --- a/idaes/core/scaling/custom_scaler_base.py +++ b/idaes/core/scaling/custom_scaler_base.py @@ -16,7 +16,6 @@ Author: Andrew Lee """ from copy import copy -from enum import Enum from pyomo.environ import ComponentMap, units, value from pyomo.core.base.units_container import UnitsError @@ -26,6 +25,7 @@ from idaes.core.scaling.scaling_base import CONFIG, ScalerBase from idaes.core.scaling.util import get_scaling_factor, NominalValueExtractionVisitor import idaes.logger as idaeslog +from idaes.core.util.misc import StrEnum # Set up logger _log = idaeslog.getLogger(__name__) @@ -41,7 +41,7 @@ } -class ConstraintScalingScheme(str, Enum): +class ConstraintScalingScheme(StrEnum): """ Schemes available for calculating constraint scaling factors. diff --git a/idaes/core/util/misc.py b/idaes/core/util/misc.py index 5daa73676d..78185256cd 100644 --- a/idaes/core/util/misc.py +++ b/idaes/core/util/misc.py @@ -15,9 +15,12 @@ This module contains miscellaneous utility functions for use in IDAES models. """ from enum import Enum +import sys import pyomo.environ as pyo from pyomo.common.config import ConfigBlock +from pyomo.core.expr.visitor import _ToStringVisitor +from pyomo.core.base.expression import ExpressionData import idaes.logger as idaeslog @@ -177,3 +180,66 @@ class StrEnum(str, Enum): def __str__(self): return str(self.value) + + +class _ToExprStringVisitor(_ToStringVisitor): + """ + Derived version of the Pyomo _ToStringVisitor class + which checks for named Expressions in the expression tree + and prints their name instead of expanding the expression tree. + """ + + def visiting_potential_leaf(self, node): + # Check if node is a named Expression + if isinstance(node, ExpressionData): + # If it is, return the name of the Expression + return True, node.name + + # Otherwise, continue descending as normal + return super().visiting_potential_leaf(node) + + +def compact_expression_to_string(expr): + """Return a compact string representation of an expression. + + Unlike the normal Pyomo string representations, this function will + identify any Expressions that appear within the expression tree and + print the name of these rather than expanding the expression tree. + + Args: + expr: ExpressionBase. The root node of an expression tree. + + Returns: + A string representation for the expression. + """ + # Create and execute the visitor pattern + # + visitor = _ToExprStringVisitor(verbose=False, smap=None) + return visitor.dfs_postorder_stack(expr) + + +def print_compact_form(expr, stream=None): + """ + Writes a compact string representation of the given component or + expression to the specified stream. + + Unlike the normal Pyomo string representations, this function will + identify any Expressions that appear within the expression tree and + print the name of these rather than expanding the expression tree. + + Args: + expr: component or expression to print + stream: StringIO object to write to. Default is stdout. + + Returns: + None + """ + if stream is None: + stream = sys.stdout + + if hasattr(expr, "expr"): + # We have a Constraint or Expression + # We want to print the expression, not the object name + expr = expr.expr + + stream.write(compact_expression_to_string(expr)) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index c7f5335c26..f701677d47 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -20,10 +20,11 @@ from operator import itemgetter import sys from inspect import signature -from math import log, isclose, inf, isfinite +from math import log, log10, isclose, inf, isfinite import json from typing import List import logging +from itertools import combinations, chain import numpy as np from scipy.linalg import svd @@ -35,6 +36,7 @@ Integers, Block, check_optimal_termination, + ComponentMap, ConcreteModel, Constraint, Expression, @@ -58,6 +60,7 @@ from pyomo.core.base.block import BlockData from pyomo.core.base.var import VarData from pyomo.core.base.constraint import ConstraintData +from pyomo.core.base.expression import ExpressionData from pyomo.repn.standard_repn import ( # pylint: disable=no-name-in-module generate_standard_repn, ) @@ -67,6 +70,8 @@ ConfigValue, document_kwargs_from_configdict, PositiveInt, + NonNegativeFloat, + NonNegativeInt, ) from pyomo.util.check_units import identify_inconsistent_units from pyomo.contrib.incidence_analysis import IncidenceGraphInterface @@ -78,8 +83,12 @@ from pyomo.common.deprecation import deprecation_warning from pyomo.common.errors import PyomoException from pyomo.common.tempfiles import TempfileManager +from pyomo.core import expr as EXPR +from pyomo.common.numeric_types import native_types +from pyomo.core.base.units_container import _PyomoUnit from idaes.core.solvers.get_solver import get_solver +from idaes.core.util.misc import compact_expression_to_string from idaes.core.util.model_statistics import ( activated_blocks_set, deactivated_blocks_set, @@ -180,7 +189,6 @@ def svd_sparse(jacobian, number_singular_values): """ u, s, vT = svds(jacobian, k=number_singular_values, which="SM") - print(u, s, vT, number_singular_values) return u, s, vT.transpose() @@ -189,7 +197,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_bounds_absolute_tolerance", ConfigValue( default=1e-4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a variable to be close " "to its bounds.", ), @@ -198,7 +206,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_bounds_relative_tolerance", ConfigValue( default=1e-4, - domain=float, + domain=NonNegativeFloat, description="Relative tolerance for considering a variable to be close " "to its bounds.", ), @@ -207,7 +215,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_bounds_violation_tolerance", ConfigValue( default=0, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a variable to violate its bounds.", doc="Absolute tolerance for considering a variable to violate its bounds. " "Some solvers relax bounds on variables thus allowing a small violation to be " @@ -218,15 +226,47 @@ def svd_sparse(jacobian, number_singular_values): "constraint_residual_tolerance", ConfigValue( default=1e-5, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance to use when checking constraint residuals.", ), ) +CONFIG.declare( + "constraint_term_mismatch_tolerance", + ConfigValue( + default=1e6, + domain=NonNegativeFloat, + description="Magnitude difference to use when checking for mismatched additive terms in constraints.", + ), +) +CONFIG.declare( + "constraint_term_cancellation_tolerance", + ConfigValue( + default=1e-4, + domain=NonNegativeFloat, + description="Absolute tolerance to use when checking for canceling additive terms in constraints.", + ), +) +CONFIG.declare( + "max_canceling_terms", + ConfigValue( + default=5, + domain=NonNegativeInt, + description="Maximum number of terms to consider when looking for canceling combinations in expressions.", + ), +) +CONFIG.declare( + "constraint_term_zero_tolerance", + ConfigValue( + default=1e-10, + domain=NonNegativeFloat, + description="Absolute tolerance to use when determining if a constraint term is equal to zero.", + ), +) CONFIG.declare( "variable_large_value_tolerance", ConfigValue( default=1e4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a value to be large.", ), ) @@ -234,7 +274,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_small_value_tolerance", ConfigValue( default=1e-4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a value to be small.", ), ) @@ -242,7 +282,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_zero_value_tolerance", ConfigValue( default=1e-8, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a value to be near to zero.", ), ) @@ -250,7 +290,7 @@ def svd_sparse(jacobian, number_singular_values): "jacobian_large_value_caution", ConfigValue( default=1e4, - domain=float, + domain=NonNegativeFloat, description="Tolerance for raising a caution for large Jacobian values.", ), ) @@ -258,7 +298,7 @@ def svd_sparse(jacobian, number_singular_values): "jacobian_large_value_warning", ConfigValue( default=1e8, - domain=float, + domain=NonNegativeFloat, description="Tolerance for raising a warning for large Jacobian values.", ), ) @@ -266,7 +306,7 @@ def svd_sparse(jacobian, number_singular_values): "jacobian_small_value_caution", ConfigValue( default=1e-4, - domain=float, + domain=NonNegativeFloat, description="Tolerance for raising a caution for small Jacobian values.", ), ) @@ -274,7 +314,7 @@ def svd_sparse(jacobian, number_singular_values): "jacobian_small_value_warning", ConfigValue( default=1e-8, - domain=float, + domain=NonNegativeFloat, description="Tolerance for raising a warning for small Jacobian values.", ), ) @@ -290,7 +330,7 @@ def svd_sparse(jacobian, number_singular_values): "parallel_component_tolerance", ConfigValue( default=1e-8, - domain=float, + domain=NonNegativeFloat, description="Tolerance for identifying near-parallel Jacobian rows/columns", ), ) @@ -298,7 +338,7 @@ def svd_sparse(jacobian, number_singular_values): "absolute_feasibility_tolerance", ConfigValue( default=1e-6, - domain=float, + domain=NonNegativeFloat, description="Feasibility tolerance for identifying infeasible constraints and bounds", ), ) @@ -336,7 +376,7 @@ def svd_sparse(jacobian, number_singular_values): "singular_value_tolerance", ConfigValue( default=1e-6, - domain=float, + domain=NonNegativeFloat, description="Tolerance for defining a small singular value", ), ) @@ -344,7 +384,7 @@ def svd_sparse(jacobian, number_singular_values): "size_cutoff_in_singular_vector", ConfigValue( default=0.1, - domain=float, + domain=NonNegativeFloat, description="Size below which to ignore constraints and variables in " "the singular vector", ), @@ -371,7 +411,7 @@ def svd_sparse(jacobian, number_singular_values): "M", # TODO: Need better name ConfigValue( default=1e5, - domain=float, + domain=NonNegativeFloat, description="Maximum value for nu in MILP models.", ), ) @@ -379,7 +419,7 @@ def svd_sparse(jacobian, number_singular_values): "m_small", # TODO: Need better name ConfigValue( default=1e-5, - domain=float, + domain=NonNegativeFloat, description="Smallest value for nu to be considered non-zero in MILP models.", ), ) @@ -387,7 +427,7 @@ def svd_sparse(jacobian, number_singular_values): "trivial_constraint_tolerance", ConfigValue( default=1e-6, - domain=float, + domain=NonNegativeFloat, description="Tolerance for identifying non-zero rows in Jacobian.", ), ) @@ -1052,6 +1092,239 @@ def display_near_parallel_variables(self, stream=None): # TODO: Block triangularization analysis # Number and size of blocks, polynomial degree of 1x1 blocks, simple pivot check of moderate sized sub-blocks? + def _collect_constraint_mismatches(self, descend_into=True): + """ + Call ConstraintTermAnalysisVisitor on all Constraints in model to walk expression + tree and collect any instances of sum expressions with mismatched terms or potential + cancellations. + + Args: + descend_into: whether to descend_into child_blocks + + Returns: + List of strings summarising constraints with mismatched terms + List of strings summarising constraints with cancellations + List of strings with constraint names where constraint contains no free variables + """ + walker = ConstraintTermAnalysisVisitor( + term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance, + term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance, + term_zero_tolerance=self.config.constraint_term_zero_tolerance, + # for the high level summary, we only need to know if there are any cancellations, + # but don't need to find all of them + max_cancellations_per_node=1, + max_canceling_terms=self.config.max_canceling_terms, + ) + + mismatch = [] + cancellation = [] + constant = [] + + for c in self._model.component_data_objects( + Constraint, descend_into=descend_into + ): + _, expr_mismatch, expr_cancellation, expr_constant, _ = ( + walker.walk_expression(c.expr) + ) + + if len(expr_mismatch) > 0: + mismatch.append(f"{c.name}: {len(expr_mismatch)} mismatched term(s)") + + if len(expr_cancellation) > 0: + cancellation.append( + f"{c.name}: {len(expr_cancellation)} potential canceling term(s)" + ) + + if expr_constant: + constant.append(c.name) + + return mismatch, cancellation, constant + + def display_constraints_with_mismatched_terms(self, stream=None): + """ + Display constraints in model which contain additive terms of significantly different magnitude. + + Args: + stream: I/O object to write report to (default = stdout) + + Returns: + None + + """ + if stream is None: + stream = sys.stdout + + mismatch, _, _ = self._collect_constraint_mismatches() + + # Write the output + _write_report_section( + stream=stream, + lines_list=mismatch, + title="The following constraints have mismatched terms:", + end_line="Call display_problematic_constraint_terms(constraint) for more information.", + header="=", + footer="=", + ) + + def display_constraints_with_canceling_terms(self, stream=None): + """ + Display constraints in model which contain additive terms which potentially cancel each other. + + Note that this method looks a the current state of the constraint, and will flag terms as + cancelling if you have a form A == B + C where C is significantly smaller than A and B. In some + cases this behavior is intended, as C is a correction term which happens to be very + small at the current state. However, you should review these constraints to determine whether + the correction term is important for the situation you are modeling and consider removing the + term if it will never be significant. + + Args: + stream: I/O object to write report to (default = stdout) + + Returns: + None + + """ + if stream is None: + stream = sys.stdout + + _, canceling, _ = self._collect_constraint_mismatches() + + # Write the output + _write_report_section( + stream=stream, + lines_list=canceling, + title="The following constraints have canceling terms:", + end_line="Call display_problematic_constraint_terms(constraint) for more information.", + header="=", + footer="=", + ) + + def display_problematic_constraint_terms( + self, constraint, max_cancellations: int = 5, stream=None + ): + """ + Display a summary of potentially problematic terms in a given constraint. + + Note that this method looks a the current state of the constraint, and will flag terms as + cancelling if you have a form A == B + C where C is significantly smaller than A and B. In some + cases this behavior is intended, as C is a correction term which happens to be very + small at the current state. However, you should review these constraints to determine whether + the correction term is important for the situation you are modeling and consider removing the + term if it will never be significant. + + Args: + constraint: ConstraintData object to be examined + max_cancellations: maximum number of cancellations per node before termination. + None = find all cancellations. + stream: I/O object to write report to (default = stdout) + + Returns: + None + + """ + if stream is None: + stream = sys.stdout + + # Check that constraint is of correct type to give useful error message + if not isinstance(constraint, ConstraintData): + # Wrong type, check if it is an indexed constraint + if isinstance(constraint, Constraint): + raise TypeError( + f"{constraint.name} is an IndexedConstraint. Please provide " + f"an individual element of {constraint.name} (ConstraintData) " + "to be examined for problematic terms." + ) + else: + # Not a constraint + raise TypeError( + f"{constraint.name} is not an instance of a Pyomo Constraint." + ) + + walker = ConstraintTermAnalysisVisitor( + term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance, + term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance, + term_zero_tolerance=self.config.constraint_term_zero_tolerance, + max_cancellations_per_node=max_cancellations, + max_canceling_terms=self.config.max_canceling_terms, + ) + + _, expr_mismatch, expr_cancellation, _, tripped = walker.walk_expression( + constraint.expr + ) + + # Combine mismatches and cancellations into a summary list + issues = [] + for k, v in expr_mismatch.items(): + tag = " " + if isinstance(k, ExpressionData): + # For clarity, if the problem node is a named Expression, note this in output + tag = " Expression " + # Want to show full expression node plus largest and smallest magnitudes + issues.append( + f"Mismatch in{tag}{compact_expression_to_string(k)} (Max {v[0]}, Min {v[1]})" + ) + # Collect summary of cancelling terms for user + # Walker gives us back a list of nodes with cancelling terms + for k, v in expr_cancellation.items(): + # Each node may have multiple cancellations, these are given as separate tuples + tag = " " + if isinstance(k, ExpressionData): + # For clarity, if the problem node is a named Expression, note this in output + tag = " Expression " + for i in v: + # For each cancellation, iterate over contributing terms and write a summary + terms = "" + for j in i: + if len(terms) > 0: + terms += ", " + # +1 to switch from 0-index to 1-index + terms += f"{j[0]+1} ({j[1]})" + issues.append( + f"Cancellation in{tag}{compact_expression_to_string(k)}. Terms {terms}" + ) + + if tripped: + end_line = ( + f"Number of canceling terms per node limited to {max_cancellations}." + ) + else: + end_line = None + + # Write the output + _write_report_section( + stream=stream, + lines_list=issues, + title=f"The following terms in {constraint.name} are potentially problematic:", + end_line=end_line, + header="=", + footer="=", + ) + + def display_constraints_with_no_free_variables(self, stream=None): + """ + Display constraints in model which contain no free variables. + + Args: + stream: I/O object to write report to (default = stdout) + + Returns: + None + + """ + if stream is None: + stream = sys.stdout + + _, _, constant = self._collect_constraint_mismatches() + + # Write the output + _write_report_section( + stream=stream, + lines_list=constant, + title="The following constraints have no free variables:", + header="=", + footer="=", + ) + def _collect_structural_warnings( self, ignore_evaluation_errors=False, ignore_unit_consistency=False ): @@ -1335,6 +1608,28 @@ def _collect_numerical_cautions(self, jac=None, nlp=None): cstring = "Variable" cautions.append(f"Caution: {len(none_value)} {cstring} with None value") + # Constraints with possible ill-posed terms + mismatch, cancellation, constant = self._collect_constraint_mismatches() + if len(mismatch) > 0: + cstring = "Constraints" + if len(mismatch) == 1: + cstring = "Constraint" + cautions.append(f"Caution: {len(mismatch)} {cstring} with mismatched terms") + if len(cancellation) > 0: + cstring = "Constraints" + if len(cancellation) == 1: + cstring = "Constraint" + cautions.append( + f"Caution: {len(cancellation)} {cstring} with potential cancellation of terms" + ) + if len(constant) > 0: + cstring = "Constraints" + if len(constant) == 1: + cstring = "Constraint" + cautions.append( + f"Caution: {len(constant)} {cstring} with no free variables" + ) + # Extreme Jacobian rows and columns jac_col = extreme_jacobian_columns( jac=jac, @@ -2777,8 +3072,6 @@ def eq_degenerate(m_dh, v): # This variable does not appear in any constraint return Constraint.Skip - m_dh.pprint() - m_dh.degenerate = Constraint(m_dh.V, rule=eq_degenerate) # When y_pos = 1, nu >= m_small @@ -4122,3 +4415,457 @@ def _collect_model_statistics(model): ) return stats + + +class ConstraintTermAnalysisVisitor(EXPR.StreamBasedExpressionVisitor): + """ + Expression walker for checking Constraints for problematic terms. + + This walker will walk the expression and look for summation terms + with mismatched magnitudes or potential cancellations. + + Args: + term_mismatch_tolerance: tolerance to use when determining mismatched + terms + term_cancellation_tolerance: tolerance to use when identifying + possible cancellation of terms + term_zero_tolerance: tolerance for considering terms equal to zero + max_canceling_terms: maximum number of terms to consider when looking + for canceling combinations (None = consider all possible combinations) + max_cancellations_per_node: maximum number of cancellations to collect + for a single node. Collection will terminate when this many cancellations + have been identified (None = collect all cancellations) + + Returns: + list of values for top-level summation terms + list of terms with mismatched magnitudes + list of terms with potential cancellations + bool indicating whether expression is a constant + """ + + def __init__( + self, + term_mismatch_tolerance: float = 1e6, + term_cancellation_tolerance: float = 1e-4, + term_zero_tolerance: float = 1e-10, + max_canceling_terms: int = 4, + max_cancellations_per_node: int = 5, + ): + super().__init__() + + # Tolerance attributes + self._log_mm_tol = log10(term_mismatch_tolerance) + self._sum_tol = term_cancellation_tolerance + self._zero_tolerance = term_zero_tolerance + self._max_canceling_terms = max_canceling_terms + self._max_cancellations_per_node = max_cancellations_per_node + + # Placeholders for collecting results + self.canceling_terms = ComponentMap() + self.mismatched_terms = ComponentMap() + + # Flag for if cancellation collection hit limit + self._cancellation_tripped = False + + def _get_value_for_sum_subexpression(self, child_data): + # child_data is a tuple, with the 0-th element being the node values + if isinstance(child_data[0][0], str): + # Values may be a list containing a string in some cases (e.g. external functions) + # Return the string in this case + return child_data[0][0] + return sum(i for i in child_data[0]) + + def _generate_combinations(self, inputs, equality=False): + # We want to test all combinations of terms for cancellation + + # The number of combinations we check depends on whether this is an (in)equality + # expression or a sum node deeper in the expression tree. + # We expect (in)equalities to generally sum to 0 (0 == expr) thus we want to + # check if any subset of the sum terms sum to zero (i.e. are any terms unnecessary). + # For other sum nodes, we need to check for any combination of terms. + + # Maximum number of terms to include in combinations + max_comb = len(inputs) + if equality: + # Subtract 1 if (in)equality node + max_comb += -1 + # We also have a limit on the maximum number of terms to consider + if self._max_canceling_terms is not None: + max_comb = min(max_comb, self._max_canceling_terms) + + # Single terms cannot cancel, thus we want all combinations of length 2 to max terms + # Note the need for +1 due to way range works + combo_range = range(2, max_comb + 1) + + # Collect combinations of terms in an iterator + # In order to identify the terms in each cancellation, we will pair each value + # with its index in the input set as a tuple using enumerate + for i in chain.from_iterable( + combinations(enumerate(inputs), r) for r in combo_range + ): + # Yield each combination in the set + yield i + + def _check_sum_cancellations(self, values_list, equality=False): + # First, strip any terms with value 0 as they do not contribute to cancellation + # We do this to keep the number of possible combinations as small as possible + stripped = [i for i in values_list if abs(i) >= self._zero_tolerance] + + cancellations = [] + + if len(stripped) == 0: + # If the stripped list is empty, there are no non-zero terms + # We can stop here and return False as there are no possible cancellations + return cancellations + + # For scaling of tolerance, we want to compare to the largest absolute value of + # the input values + max_value = abs(max(stripped, key=abs)) + + for i in self._generate_combinations(stripped, equality): + # Generate combinations will return a list of combinations of input terms + # each element of the list will be a 2-tuple representing a term in the + # input list with the first value being the position in the input set and + # the second being the value + + # Check if the sum of values in the combination are below tolerance + if abs(sum(j[1] for j in i)) <= self._sum_tol * max_value: + # If so, record combination as canceling + cancellations.append(i) + + # Terminate loop if we have reached the max cancellations to collect + if len(cancellations) >= self._max_cancellations_per_node: + self._cancellation_tripped = True + break + + return cancellations + + def _perform_checks(self, node, child_data): + # Perform checks for problematic expressions + # First, need to check to see if any child data is a list + # This indicates a sum expression + const = True + + for d in child_data: + # We will check for canceling terms here, rather than the sum itself, to handle special cases + # We want to look for cases where a sum term results in a value much smaller + # than the terms of the sum + # Each element of child_data is a tuple where the 0-th element is the node values + if isinstance(d[0][0], str): + # Values may be a list containing a string in some cases (e.g. external functions) + # Skip if this is the case + pass + else: + for c in self._check_sum_cancellations(d[0]): + if node in self.canceling_terms.keys(): + self.canceling_terms[node].append(c) + else: + self.canceling_terms[node] = [c] + + # Expression is not constant if any child is not constant + # Element 1 is a bool indicating if the child node is constant + if not d[1]: + const = False + + # Return any problematic terms found + return const + + def _check_base_type(self, node): + if isinstance(node, VarData): + const = node.fixed + else: + const = True + return [value(node)], const, False + + def _check_equality_expression(self, node, child_data): + # (In)equality expressions are a special case of sum expressions + # child_data has two elements; 0 is the LHS of the (in)equality and 1 is the RHS + # Each of these then contains three elements; 0 is a list of values for the sum components, + # 1 is a bool indicating if the child node term is constant, and 2 is a bool indicating if + # the child ode is a sum expression. + + # First, to check for cancellation we need to negate one side of the expression + # mdata will contain the new set of child_data with negated values + mdata = [] + # child_data[0][0] has the values of the LHS of the (in)equality, and we will negate these + vals = [] + for j in child_data[0][0]: + vals.append(-j) + # Append the negated values along with whether the node is constant to mdata + mdata.append((vals, child_data[0][1])) + # child_data[1] is the RHS, so we take this as it appears and append to mdata + mdata.append(child_data[1]) + + # Next, call the method to check the sum expression + vals, const, _ = self._check_sum_expression(node, mdata) + + # Next, we need to check for canceling terms. + # child_data[x][2] indicates if a node is a sum expression or not + if not child_data[0][2] and not child_data[1][2]: + # If both sides are not sum expressions, we have the form a == b + # Simple lLinking constraints do not need to be checked + pass + # Next, we can ignore any term that has already been flagged as mismatched + elif node in self.mismatched_terms.keys(): + pass + # We can also ignore any case where one side of the (in)equality is constant + # I.e. if either child_node[x][1] is True + elif any(d[1] for d in mdata): + pass + else: + # Check for cancellation + # First, collect terms from both sides + # Note: outer loop comes first in list comprehension + t = [i for d in mdata for i in d[0]] + + # Then check for cancellations + for c in self._check_sum_cancellations(t, equality=True): + if node in self.canceling_terms.keys(): + self.canceling_terms[node].append(c) + else: + self.canceling_terms[node] = [c] + + return vals, const, False + + def _check_product_expr(self, node, child_data): + # We expect child_data to be a tuple of len 2 (a*b) + # If both or neither a and b are sums (xor), handle like other expressions + if not child_data[0][2] ^ child_data[1][2]: + return self._check_general_expr(node, child_data) + else: + # Here we have the case of a sum and a multiplier + # For this case, we will make the result look like a sum with the + # multiplier applied + # This is important for scaling of the sum terms, as cancellation should be + # Checked on the scaled value + + # First, check if both terms are constants - if so we can just get the value of + # this node and move on. + if child_data[0][1] and child_data[1][1]: + return self._check_general_expr(node, child_data) + + # The length of the values (child_data[i][0]) of the multiplier will be 1 + # We can just iterate over all terms in both value lists and multiply + # each pair + vals = [i * j for i in child_data[0][0] for j in child_data[1][0]] + + # Return the list of values, not constant, is a sum expression (apparent) + return vals, False, True + + def _check_div_expr(self, node, child_data): + # We expect child_data to be a tuple of len 2 (a/b) + # If the numerator is not a sum, handle like general expression + # child_data[0] is numerator, child_data[1] is denominator + if not child_data[0][2]: + return self._check_general_expr(node, child_data) + else: + # If the numerator is a sum, we will treat this as if the division + # were applied to each term in the numerator separately + # This is important for scaling of the sum terms, as cancellation should be + # Checked on the scaled value + + # First, check if the numerator is constant - if so we can just get the value of + # this node and move on. + if child_data[0][1]: + return self._check_general_expr(node, child_data) + + # Next, we need to get the value of the denominator + denom = self._get_value_for_sum_subexpression(child_data[1]) + + try: + vals = [i / denom for i in child_data[0][0]] + except ZeroDivisionError: + raise ZeroDivisionError( + f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 " + f"({str(node)})." + ) + + # Return the list of values, not constant, is a sum expression (apparent) + return vals, False, True + + def _check_negation_expr(self, node, child_data): + # Here we need to defer checking of cancellations too due to different ways + # these can appear in an expression. + # We will simply negate all values on the child node values (child_data[0][0]) + # and pass on the rest. + vals = [-i for i in child_data[0][0]] + + # Return the list of values, not constant, is a sum expression (apparent) + return vals, child_data[0][1], child_data[0][2] + + def _check_general_expr(self, node, child_data): + const = self._perform_checks(node, child_data) + + try: + # pylint: disable-next=protected-access + val = node._apply_operation( + list(map(self._get_value_for_sum_subexpression, child_data)) + ) + except ValueError: + raise ValueError( + f"Error in ConstraintTermAnalysisVisitor: error evaluating {str(node)}." + ) + except ZeroDivisionError: + raise ZeroDivisionError( + f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 " + f"({str(node)})." + ) + except Exception as err: + # Catch and re-raise any other error that occurs + _log.exception(f"Unexpected {err=}, {type(err)=}") + raise + + return [val], const, False + + def _check_other_expression(self, node, child_data): + const = self._perform_checks(node, child_data) + + # First, need to get value of input terms, which may be sub-expressions + input_mag = [] + for i in child_data: + input_mag.append(self._get_value_for_sum_subexpression(i)) + + # Next, create a copy of the function with expected magnitudes as inputs + newfunc = node.create_node_with_local_data(input_mag) + + # Evaluate new function and return the value along with check results + return [value(newfunc)], const, False + + def _check_ranged_expression(self, node, child_data): + # child_data should have 3 elements, LHS, middle term, and RHS + lhs_vals, lhs_const, _ = self._check_equality_expression(node, child_data[:2]) + rhs_vals, rhs_const, _ = self._check_equality_expression(node, child_data[1:]) + + # Constant is a bit vague in terms ranged expressions. + # We will call the term constant only if all parts are constant + const = lhs_const and rhs_const + + # For values, we need to avoid double counting the middle term + # Also for sign convention, we will negate the outer terms + vals = lhs_vals + [-rhs_vals[1]] + + return vals, const, False + + def _check_sum_expression(self, node, child_data): + # Sum expressions need special handling + # For sums, collect all child values into a list + vals = [] + # We will check for cancellation in this node at the next level + # Pyomo is generally good at simplifying compound sums + const = True + # Collect data from child nodes + for d in child_data: + vals.append(self._get_value_for_sum_subexpression(d)) + + # Expression is not constant if any child is not constant + # Element 1 is a bool indicating if node is constant + if not d[1]: + const = False + + # Check for mismatched terms + if len(vals) > 1: + absvals = [abs(v) for v in vals if abs(v) >= self._zero_tolerance] + if len(absvals) > 0: + vl = max(absvals) + vs = min(absvals) + if vl != vs: + if vs == 0: + diff = log10(vl) + else: + diff = log10(vl / vs) + + if diff >= self._log_mm_tol: + self.mismatched_terms[node] = (vl, vs) + + return vals, const, True + + node_type_method_map = { + EXPR.EqualityExpression: _check_equality_expression, + EXPR.InequalityExpression: _check_equality_expression, + EXPR.RangedExpression: _check_ranged_expression, + EXPR.SumExpression: _check_sum_expression, + EXPR.NPV_SumExpression: _check_sum_expression, + EXPR.ProductExpression: _check_product_expr, + EXPR.MonomialTermExpression: _check_product_expr, + EXPR.NPV_ProductExpression: _check_product_expr, + EXPR.DivisionExpression: _check_div_expr, + EXPR.NPV_DivisionExpression: _check_div_expr, + EXPR.PowExpression: _check_general_expr, + EXPR.NPV_PowExpression: _check_general_expr, + EXPR.NegationExpression: _check_negation_expr, + EXPR.NPV_NegationExpression: _check_negation_expr, + EXPR.AbsExpression: _check_general_expr, + EXPR.NPV_AbsExpression: _check_general_expr, + EXPR.UnaryFunctionExpression: _check_general_expr, + EXPR.NPV_UnaryFunctionExpression: _check_general_expr, + EXPR.Expr_ifExpression: _check_other_expression, + EXPR.ExternalFunctionExpression: _check_other_expression, + EXPR.NPV_ExternalFunctionExpression: _check_other_expression, + EXPR.LinearExpression: _check_sum_expression, + } + + def exitNode(self, node, data): + """ + Method to call when exiting node to check for potential issues. + """ + # Return [node values], constant (bool), sum expression (bool) + # first check if the node is a leaf + nodetype = type(node) + + if nodetype in native_types: + return [node], True, False + + node_func = self.node_type_method_map.get(nodetype, None) + if node_func is not None: + return node_func(self, node, data) + + if not node.is_expression_type(): + # this is a leaf, but not a native type + if nodetype is _PyomoUnit: + # Unit have no value, so return 1 as a placeholder + return [1], True, False + + # Var or Param + return self._check_base_type(node) + # might want to add other common types here + + # not a leaf - check if it is a named expression + if ( + hasattr(node, "is_named_expression_type") + and node.is_named_expression_type() + ): + return self._check_other_expression(node, data) + + raise TypeError( + f"An unhandled expression node type: {str(nodetype)} was encountered while " + f"analyzing constraint terms {str(node)}" + ) + + def walk_expression(self, expr): + """ + Main method to call to walk an expression and return analysis. + + Args: + expr - expression to be analyzed + + Returns: + list of values of top-level additive terms + ComponentSet containing any mismatched terms + ComponentSet containing any canceling terms + Bool indicating whether expression is a constant + """ + # Create new holders for collected terms + self.canceling_terms = ComponentMap() + self.mismatched_terms = ComponentMap() + + # Call parent walk_expression method + vals, const, _ = super().walk_expression(expr) + + # Return results + return ( + vals, + self.mismatched_terms, + self.canceling_terms, + const, + self._cancellation_tripped, + ) diff --git a/idaes/core/util/tests/test_misc.py b/idaes/core/util/tests/test_misc.py index 594b790e12..363fcd8df0 100644 --- a/idaes/core/util/tests/test_misc.py +++ b/idaes/core/util/tests/test_misc.py @@ -14,16 +14,19 @@ This module contains miscellaneous utility functions for use in IDAES models. """ +from io import StringIO import pytest import re -from pyomo.environ import ConcreteModel, Set, Block, Var, units +from pyomo.environ import ConcreteModel, Constraint, Expression, Set, Block, Var, units from pyomo.common.config import ConfigBlock from pyomo.core.base.units_container import UnitsError from idaes.core.util.misc import ( add_object_reference, set_param_from_config, + compact_expression_to_string, + print_compact_form, ) import idaes.logger as idaeslog @@ -361,3 +364,54 @@ def test_unitted_default(self, caplog): "b no units provided for parameter test_param - assuming " "default units" in caplog.text ) + + +class TestToExprStringVisitor: + @pytest.mark.unit + def test_compact_expression_to_string(self): + m = ConcreteModel() + m.v1 = Var() + m.v2 = Var() + m.v3 = Var() + + m.e1 = Expression(expr=m.v1 + m.v2) + m.c1 = Constraint(expr=0 == m.v3 * m.e1) + + str = compact_expression_to_string(m.c1.expr) + expected = "v3*e1 == 0" + + assert str == expected + + @pytest.mark.unit + def test_print_compact_form_expr(self): + m = ConcreteModel() + m.v1 = Var() + m.v2 = Var() + m.v3 = Var() + + m.e1 = Expression(expr=m.v1 + m.v2) + m.c1 = Constraint(expr=0 == m.v3 * m.e1) + + expected = "v3*e1 == 0" + + stream = StringIO() + print_compact_form(m.c1.expr, stream=stream) + + assert stream.getvalue() == expected + + @pytest.mark.unit + def test_print_compact_form_component(self): + m = ConcreteModel() + m.v1 = Var() + m.v2 = Var() + m.v3 = Var() + + m.e1 = Expression(expr=m.v1 + m.v2) + m.c1 = Constraint(expr=0 == m.v3 * m.e1) + + expected = "v3*e1 == 0" + + stream = StringIO() + print_compact_form(m.c1, stream=stream) + + assert stream.getvalue() == expected diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index ec5b671e34..e10f8a7c57 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -13,52 +13,67 @@ """ This module contains model diagnostic utility functions for use in IDAES (Pyomo) models. """ +from copy import deepcopy from io import StringIO import math -import numpy as np -import pytest -import re import os -from copy import deepcopy +import re +from unittest import TestCase +import numpy as np from pandas import DataFrame +import pytest from pyomo.environ import ( Block, ConcreteModel, Constraint, Expression, - log, - tan, - asin, - acos, - sqrt, + Expr_if, Objective, PositiveIntegers, Set, SolverFactory, - Suffix, - TransformationFactory, units, value, Var, assert_optimal_termination, Param, Integers, + exp, + log, + log10, + sin, + cos, + tan, + asin, + acos, + atan, + sqrt, + sinh, + cosh, + tanh, + asinh, + acosh, + atanh, ) from pyomo.common.collections import ComponentSet from pyomo.contrib.pynumero.asl import AmplInterface from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.common.fileutils import this_file_dir from pyomo.common.tempfiles import TempfileManager +from pyomo.core import expr as EXPR import idaes.core.util.scaling as iscale import idaes.logger as idaeslog from idaes.core.solvers import get_solver from idaes.core import FlowsheetBlock from idaes.core.util.testing import PhysicalParameterTestBlock -from unittest import TestCase - +from idaes.models.properties.modular_properties.eos.ceos_common import ( + cubic_roots_available, + CubicThermoExpressions, + CubicType as CubicEoS, +) from idaes.core.util.model_diagnostics import ( DiagnosticsToolbox, SVDToolbox, @@ -81,6 +96,7 @@ _collect_model_statistics, check_parallel_jacobian, compute_ill_conditioning_certificate, + ConstraintTermAnalysisVisitor, ) from idaes.core.util.parameter_sweep import ( SequentialSweepRunner, @@ -90,6 +106,7 @@ UniformSampling, ) from idaes.core.util.testing import _enable_scip_solver_for_testing +from idaes.models.properties import iapws95 __author__ = "Alex Dowling, Douglas Allan, Andrew Lee" @@ -999,6 +1016,232 @@ def test_display_near_parallel_variables(self): v1, v4 v2, v4 +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_collect_constraint_mismatches(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + + # Constraint with no free variables + m.c1 = Constraint(expr=m.v1 == m.v2) + m.v1.fix() + m.v2.fix() + + # Constraint with mismatched terms + m.v3 = Var(initialize=10) + m.v4 = Var(initialize=10) + m.v5 = Var(initialize=1e-6) + m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5) + + # Constraint with cancellation + m.v6 = Var(initialize=10) + m.c3 = Constraint(expr=m.v6 == 10 + m.v3 - m.v4) + + dt = DiagnosticsToolbox(model=m) + + mismatch, cancellation, constant = dt._collect_constraint_mismatches() + + assert mismatch == ["c2: 1 mismatched term(s)"] + assert cancellation == [ + "c2: 1 potential canceling term(s)", + "c3: 1 potential canceling term(s)", + ] + assert constant == ["c1"] + + @pytest.mark.component + def test_display_problematic_constraint_terms(self): + m = ConcreteModel() + + # Constraint with mismatched terms + m.v3 = Var(initialize=10) + m.v4 = Var(initialize=10) + m.v5 = Var(initialize=1e-6) + m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5) + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_problematic_constraint_terms(m.c2, stream=stream) + + expected = """==================================================================================== +The following terms in c2 are potentially problematic: + + Mismatch in v4 + v5 (Max 10, Min 1e-06) + Cancellation in v3 == v4 + v5. Terms 1 (-10), 2 (10) + +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_problematic_constraint_terms_limited(self): + m = ConcreteModel() + m.v1 = Var(initialize=10) + m.v2 = Var(initialize=-10) + m.v4 = Var(initialize=10) + + # Constraint with cancellation + m.v6 = Var(initialize=10) + m.c1 = Constraint(expr=m.v6 == 10 - m.v4 + m.v1 + m.v2) + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_problematic_constraint_terms(m.c1, stream=stream) + + expected = """==================================================================================== +The following terms in c1 are potentially problematic: + + Cancellation in v6 == 10 - v4 + v1 + v2. Terms 1 (-10), 2 (10) + Cancellation in v6 == 10 - v4 + v1 + v2. Terms 1 (-10), 4 (10) + Cancellation in v6 == 10 - v4 + v1 + v2. Terms 2 (10), 3 (-10) + Cancellation in v6 == 10 - v4 + v1 + v2. Terms 2 (10), 5 (-10) + Cancellation in v6 == 10 - v4 + v1 + v2. Terms 3 (-10), 4 (10) + +Number of canceling terms per node limited to 5. +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_problematic_constraint_terms_indexed_error(self): + m = ConcreteModel() + m.s = Set(initialize=[1, 2]) + m.v1 = Var(m.s, initialize=2) + m.v2 = Var(m.s, initialize=3) + + def _c_rule(b, i): + return b.v1[i] == b.v2[i] + + m.c1 = Constraint(m.s, rule=_c_rule) + + dt = DiagnosticsToolbox(model=m) + + # Check for indexed constraints + with pytest.raises( + TypeError, + match=re.escape( + "c1 is an IndexedConstraint. Please provide " + "an individual element of c1 (ConstraintData) " + "to be examined for problematic terms." + ), + ): + dt.display_problematic_constraint_terms(m.c1) + + # Check for not a constraints + with pytest.raises( + TypeError, match="v1 is not an instance of a Pyomo Constraint." + ): + dt.display_problematic_constraint_terms(m.v1) + + @pytest.mark.component + def test_display_problematic_constraint_terms_named_expression(self): + m = ConcreteModel() + + # Constraint with mismatched terms + m.v1 = Var(initialize=10) + m.v2 = Var(initialize=10) + m.v3 = Var(initialize=1e-6) + + m.e1 = Expression(expr=(m.v1 - m.v2)) + + m.c2 = Constraint(expr=m.v3 == m.e1**2) + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_problematic_constraint_terms(m.c2, stream=stream) + + expected = """==================================================================================== +The following terms in c2 are potentially problematic: + + Cancellation in Expression e1. Terms 1 (10), 2 (-10) + +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_constraints_with_mismatched_terms(self): + m = ConcreteModel() + # Constraint with mismatched terms + m.v3 = Var(initialize=10) + m.v4 = Var(initialize=10) + m.v5 = Var(initialize=1e-6) + m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5) + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_constraints_with_mismatched_terms(stream=stream) + + expected = """==================================================================================== +The following constraints have mismatched terms: + + c2: 1 mismatched term(s) + +Call display_problematic_constraint_terms(constraint) for more information. +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_constraints_with_canceling_terms(self): + m = ConcreteModel() + # Constraint with mismatched terms + m.v3 = Var(initialize=10) + m.v4 = Var(initialize=10) + + # Constraint with cancellation + m.v6 = Var(initialize=10) + m.c3 = Constraint(expr=m.v6 == 10 + m.v3 - m.v4) + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_constraints_with_canceling_terms(stream=stream) + + expected = """==================================================================================== +The following constraints have canceling terms: + + c3: 1 potential canceling term(s) + +Call display_problematic_constraint_terms(constraint) for more information. +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_display_constraints_with_no_free_variables(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + + # Constraint with no free variables + m.c1 = Constraint(expr=m.v1 == m.v2) + m.v1.fix() + m.v2.fix() + + dt = DiagnosticsToolbox(model=m) + + stream = StringIO() + dt.display_constraints_with_no_free_variables(stream=stream) + + expected = """==================================================================================== +The following constraints have no free variables: + + c1 + ==================================================================================== """ @@ -1151,6 +1394,7 @@ def test_collect_numerical_cautions(self, model): dt = DiagnosticsToolbox(model=model.b) cautions = dt._collect_numerical_cautions() + assert len(cautions) == 5 assert ( "Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" @@ -1440,6 +1684,58 @@ def test_report_numerical_issues(self, model): compute_infeasibility_explanation() display_variables_at_or_outside_bounds() +==================================================================================== +""" + + assert stream.getvalue() == expected + + @pytest.mark.component + def test_report_numerical_issues_cancellation(self): + model = ConcreteModel() + + model.v1 = Var(initialize=1) + model.v2 = Var(initialize=2) + model.v3 = Var(initialize=1e-8) + + # Non-problematic constraints + model.c1 = Constraint(expr=2 * model.v1 == model.v2) + + # Problematic constraints + model.c10 = Constraint(expr=0 == exp(model.v1 - 0.5 * model.v2)) + model.c11 = Constraint(expr=0 == model.v1 - 0.5 * model.v2 + model.v3) + + dt = DiagnosticsToolbox(model=model) + + stream = StringIO() + dt.report_numerical_issues(stream) + + expected = """==================================================================================== +Model Statistics + + Jacobian Condition Number: Undefined (Exactly Singular) + +------------------------------------------------------------------------------------ +3 WARNINGS + + WARNING: 1 Constraint with large residuals (>1.0E-05) + WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-08) + WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08) + +------------------------------------------------------------------------------------ +3 Cautions + + Caution: 1 Variable with value close to zero (tol=1.0E-08) + Caution: 1 Constraint with mismatched terms + Caution: 1 Constraint with potential cancellation of terms + +------------------------------------------------------------------------------------ +Suggested next steps: + + display_constraints_with_large_residuals() + compute_infeasibility_explanation() + display_near_parallel_constraints() + display_near_parallel_variables() + ==================================================================================== """ @@ -2267,7 +2563,6 @@ def test_build_model(self, model): ca = IpoptConvergenceAnalysis(model) clone = ca._build_model() - clone.pprint() assert clone is not model assert isinstance(clone.v1, Var) @@ -4053,3 +4348,1057 @@ def test_output(self, model): Constraints / bounds in guards for stability: """ assert expected in stream.getvalue() + + +class TestConstraintTermAnalysisVisitor: + @pytest.mark.unit + def test_sum_combinations(self): + # Check method to generate sums of all combinations of terms + # excludes single term sums + terms = [1, 2, 3, 4, 5] + visitor = ConstraintTermAnalysisVisitor(max_canceling_terms=None) + sums = [i for i in visitor._generate_combinations(terms)] + + print(sums) + + expected = [ + ((0, 1), (1, 2)), + ((0, 1), (2, 3)), + ((0, 1), (3, 4)), + ((0, 1), (4, 5)), + ((1, 2), (2, 3)), + ((1, 2), (3, 4)), + ((1, 2), (4, 5)), + ((2, 3), (3, 4)), + ((2, 3), (4, 5)), + ((3, 4), (4, 5)), + ((0, 1), (1, 2), (2, 3)), + ((0, 1), (1, 2), (3, 4)), + ((0, 1), (1, 2), (4, 5)), + ((0, 1), (2, 3), (3, 4)), + ((0, 1), (2, 3), (4, 5)), + ((0, 1), (3, 4), (4, 5)), + ((1, 2), (2, 3), (3, 4)), + ((1, 2), (2, 3), (4, 5)), + ((1, 2), (3, 4), (4, 5)), + ((2, 3), (3, 4), (4, 5)), + ((0, 1), (1, 2), (2, 3), (3, 4)), + ((0, 1), (1, 2), (2, 3), (4, 5)), + ((0, 1), (1, 2), (3, 4), (4, 5)), + ((0, 1), (2, 3), (3, 4), (4, 5)), + ((1, 2), (2, 3), (3, 4), (4, 5)), + ((0, 1), (1, 2), (2, 3), (3, 4), (4, 5)), + ] + + assert sums == expected + + @pytest.mark.unit + def test_sum_combinations_limited_default(self): + # Check method to generate sums of all combinations of terms + # excludes single term sums + terms = [1, 2, 3, 4, 5] + visitor = ConstraintTermAnalysisVisitor() + sums = [i for i in visitor._generate_combinations(terms)] + + expected = expected = [ + ((0, 1), (1, 2)), + ((0, 1), (2, 3)), + ((0, 1), (3, 4)), + ((0, 1), (4, 5)), + ((1, 2), (2, 3)), + ((1, 2), (3, 4)), + ((1, 2), (4, 5)), + ((2, 3), (3, 4)), + ((2, 3), (4, 5)), + ((3, 4), (4, 5)), + ((0, 1), (1, 2), (2, 3)), + ((0, 1), (1, 2), (3, 4)), + ((0, 1), (1, 2), (4, 5)), + ((0, 1), (2, 3), (3, 4)), + ((0, 1), (2, 3), (4, 5)), + ((0, 1), (3, 4), (4, 5)), + ((1, 2), (2, 3), (3, 4)), + ((1, 2), (2, 3), (4, 5)), + ((1, 2), (3, 4), (4, 5)), + ((2, 3), (3, 4), (4, 5)), + ((0, 1), (1, 2), (2, 3), (3, 4)), + ((0, 1), (1, 2), (2, 3), (4, 5)), + ((0, 1), (1, 2), (3, 4), (4, 5)), + ((0, 1), (2, 3), (3, 4), (4, 5)), + ((1, 2), (2, 3), (3, 4), (4, 5)), + # 5 term combination should be excluded + ] + + assert sums == expected + + @pytest.mark.unit + def test_check_sum_cancellations(self): + terms = [1, -2, 3, -4, 5] + visitor = ConstraintTermAnalysisVisitor() + cancellations = visitor._check_sum_cancellations(terms) + + # We expect to canceling combinations + # Results should be a list with each entry being a tuple containing a + # canceling combination + # In turn, each element of the tuple should be a 2-tuple with the + # position of the term in the input list and its value + expected = [((0, 1), (2, 3), (3, -4)), ((0, 1), (1, -2), (3, -4), (4, 5))] + + assert cancellations == expected + + @pytest.mark.unit + def test_int(self): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=7) + + assert vv == [7] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_float(self): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=7.7) + + assert vv == [7.7] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_scalar_param(self): + m = ConcreteModel() + m.scalar_param = Param(initialize=1, mutable=True) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_param + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_indexed_param(self): + m = ConcreteModel() + m.indexed_param = Param(["a", "b"], initialize=1, mutable=True) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_param["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_scalar_var(self): + m = ConcreteModel() + m.scalar_var = Var(initialize=1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_var + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_indexed_var(self): + m = ConcreteModel() + m.indexed_var = Var(["a", "b"], initialize=1) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_var["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_scalar_var_fixed(self): + m = ConcreteModel() + m.scalar_var = Var(initialize=1) + m.scalar_var.fix() + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_var + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_indexed_var_fixed(self): + m = ConcreteModel() + m.indexed_var = Var(["a", "b"], initialize=1) + m.indexed_var.fix() + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_var["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_equality_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e-7) + m.v2 = Var(initialize=1e7) + + expr = m.v1 == m.v2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_equality_expr_constant(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e-7) + m.v2 = Var(initialize=1e7) + + # Fix v1, not constant yet as v2 still free + m.v1.fix() + + expr = m.v1 == m.v2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + # Fix v2, now constant + m.v2.fix() + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v1["a"].set_value(1e-7) + + expr = sum(m.v1[i] for i in m.v1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [1e-7, 1e7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_sum_expr_constant(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v1["a"].set_value(1e-7) + m.v1.fix() + + expr = sum(m.v1[i] for i in m.v1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [1e-7, 1e7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.fixture(scope="class") + def model(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5.0000001) + m.v4 = Var(initialize=5) + + return m + + @pytest.mark.unit + def test_product_expr(self, model): + m = ConcreteModel() + expr = model.v1 * model.v2 * model.v3 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(30.0000006, rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_product_sum_expr(self, model): + expr = (model.v1 + model.v2) * (model.v3 + model.v4) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx((2 + 3) * (5.0000001 + 5), rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_product_sum_expr_w_negation(self, model): + expr = (model.v1 + model.v2) * (model.v3 - model.v4) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(0.0000005, rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + def test_division_expr(self, model): + expr = model.v1 / model.v2 / model.v3 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(2 / 3 / 5.0000001, rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_division_sum_expr(self, model): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=(model.v1 + model.v2) / (model.v3 + model.v4) + ) + + assert vv == [ + pytest.approx(2 / (5.0000001 + 5), rel=1e-8), + pytest.approx(3 / (5.0000001 + 5), rel=1e-8), + ] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_division_sum_expr_w_negation(self, model): + expr = (model.v1 - model.v2) / (model.v3 - model.v4) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [ + pytest.approx(2 / 0.0000001, rel=1e-8), + pytest.approx(-3 / 0.0000001, rel=1e-8), + ] + assert len(mm) == 0 + # Check for cancellation should be deferred + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_division_sum_expr_w_zero_denominator(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + m.v4 = Var(initialize=5) + + expr = (m.v1 - m.v2) / (m.v3 - m.v4) + with pytest.raises( + ZeroDivisionError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 " + "((v1 - v2)/(v3 - v4))." + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + @pytest.mark.unit + def test_division_sum_expr_w_negation_deferred(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=2) + m.v3 = Var(initialize=5.0000001) + m.v4 = Var(initialize=5) + + expr = ((m.v1 - m.v2) / (m.v3 - m.v4)) ** 2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(0, abs=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + def test_division_expr_error(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=0) + + with pytest.raises( + ZeroDivisionError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: found division with " + "denominator of 0 (v1/v2)." + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=m.v1 / m.v2) + + @pytest.mark.unit + def test_pow_expr(self, model): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=model.v1**model.v2 + ) + + assert vv == [8] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_pow_sum_expr(self, model): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=(model.v1 + model.v2) ** (model.v3 + model.v4) + ) + + assert vv == [pytest.approx(5**10.0000001, rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_pow_sum_expr_w_negation(self, model): + expr = (model.v1 + model.v2) ** (model.v3 - model.v4) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx((2 + 3) ** (0.0000001), rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.fixture(scope="class") + def func_model(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=0.99999) + m.v3 = Var(initialize=-100) + + m.v2.fix() + + return m + + @pytest.mark.unit + def test_negation_expr(self, func_model): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=-func_model.v1 + ) + + assert vv == [-1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_negation_sum_expr(self, func_model): + expr = -(func_model.v1 - func_model.v2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + # Checking the cancellation should be deferred + # Expect to get two values back + assert vv == [pytest.approx(-1, rel=1e-8), pytest.approx(0.99999, rel=1e-8)] + assert len(mm) == 0 + # Check is deferred, so no cancellations identified + assert len(cc) == 0 + assert not k + assert not tr + + # acosh has bounds that don't fit with other functions - we will assume we caught enough + func_list = [ + exp, + log, + log10, + sqrt, + sin, + cos, + tan, + asin, + acos, + atan, + sinh, + cosh, + tanh, + asinh, + atanh, + ] + func_error_list = [log, log10, sqrt, asin, acos, acosh, atanh] + + @pytest.mark.unit + @pytest.mark.parametrize("func", func_list) + def test_func_expr(self, func_model, func): + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=func(func_model.v2) + ) + + assert vv == [pytest.approx(value(func(0.99999)), rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + @pytest.mark.parametrize("func", func_list) + def test_func_sum_expr(self, func_model, func): + expr = func(func_model.v1 - func_model.v2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + @pytest.mark.parametrize("func", func_list) + def test_func_sum_expr_w_negation(self, func_model, func): + expr = func(func_model.v1 - func_model.v2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + @pytest.mark.parametrize("func", func_error_list) + def test_func_expr_error(self, func_model, func): + with pytest.raises( + ValueError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: error evaluating " + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=func(func_model.v3)) + + @pytest.mark.unit + def test_expr_if(self, model): + model.exprif = Expr_if( + IF=model.v1, + THEN=model.v1 + model.v2, + ELSE=model.v3 - model.v4, + ) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=model.exprif + ) + + assert vv == [pytest.approx(5, rel=1e-8)] + assert len(mm) == 0 + assert model.exprif in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + @pytest.mark.skipif( + not AmplInterface.available(), reason="pynumero_ASL is not available" + ) + @pytest.mark.skipif(not cubic_roots_available, reason="Cubic roots not available") + def test_ext_func(self): + # Use the cubic root external function to test + m = ConcreteModel() + m.a = Var(initialize=1) + m.b = Var(initialize=1) + + m.expr_write = CubicThermoExpressions(m) + Z = m.expr_write.z_liq(eos=CubicEoS.PR, A=m.a, B=m.b) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=Z) + + assert vv == [pytest.approx(-2.1149075414767577, rel=1e-8)] + assert len(mm) == 0 + assert Z in cc + assert len(cc) == 1 + assert not k + assert not tr + + # Check that model state did not change + assert value(m.a) == 1 + assert value(m.b) == 1 + assert value(Z) == pytest.approx(-2.1149075414767577, rel=1e-8) + + @pytest.mark.unit + def test_equality_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v2 = Var(initialize=1e-7) + + expr = m.v2 == sum(m.v1[i] for i in m.v1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 3e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_inequality_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v2 = Var(initialize=1e-7) + + expr = m.v2 <= sum(m.v1[i] for i in m.v1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 3e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_compound_equality_expr_1(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v2 = Var(initialize=1e-7) + + expr = 6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-6e-7, 2.4e8] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_ranged_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e7) + m.v2 = Var(initialize=1e-7) + m.v3 = Var(initialize=1e7) + + m.expr = EXPR.RangedExpression(args=(m.v1, m.v2, m.v3), strict=True) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr) + + assert vv == [-1e7, 1e-7, -1e7] + assert m.expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + # Fix v1 and v2 to make first two terms constant + m.v1.fix() + m.v2.fix() + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr) + + # Should not be flagged as constant due to v3 + assert vv == [-1e7, 1e-7, -1e7] + assert m.expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert not k + assert not tr + + # Fix v3 to make all terms constant + m.v3.fix() + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr) + + # Should now be constant + assert vv == [-1e7, 1e-7, -1e7] + assert m.expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert k + assert not tr + + @pytest.mark.unit + def test_compound_equality_expr_2(self): + m = ConcreteModel() + m.v1 = Var(["a", "b", "c"], initialize=1e7) + m.v2 = Var(initialize=1e-7) + m.v3 = Var(initialize=1e3) + + # Set this small so we get two mismatched warnings + m.v1["a"].set_value(1e-7) + + expr1 = sum(m.v1[i] for i in m.v1) + expr = 6 * m.v2 == 8 * expr1 + m.v3 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-6e-7, pytest.approx(8 * (2e7 + 1e-7) + 1000, rel=1e-8)] + assert expr in mm + assert expr1 in mm + assert len(mm) == 2 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_canceling_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 - m.v2 + ) + + assert vv == [2, -2] + assert len(mm) == 0 + # We do not check cancellation at the sum, so this should be empty + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_expr_w_canceling_sum(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=2) + m.v3 = Var(initialize=3) + + expr = m.v3 * (m.v1 - m.v2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [6, -6] + assert len(mm) == 0 + # Check for cancellation should be deferred + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_expr_w_deferred_canceling_sum(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=2) + m.v3 = Var(initialize=3) + + expr = (m.v3 * (m.v1 - m.v2)) ** 2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0] + assert len(mm) == 0 + # We should get a warning about canceling sums here + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + # Check for tolerance of sum cancellation + m.v2.set_value(2.00022) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor( + term_cancellation_tolerance=1e-4 + ).walk_expression(expr=expr) + + assert vv == [pytest.approx((3 * -0.00022) ** 2, rel=1e-8)] + assert len(mm) == 0 + # This should pass as the difference is greater than tol + assert len(cc) == 0 + assert not k + assert not tr + + # Change tolerance so it should identify cancellation + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor( + term_cancellation_tolerance=1e-3 + ).walk_expression(expr=expr) + + assert vv == [pytest.approx((3 * -0.00022) ** 2, rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + def test_canceling_equality_expr_safe(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e7) + m.v2 = Var(initialize=1e7) + + # This is a standard constraint form, so should have no issues despite cancellation + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=0 == m.v1 - m.v2 + ) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_canceling_equality_expr_zero_term(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e7) + m.v2 = Var(initialize=1e7) + m.v3 = Var(initialize=0) + + expr = m.v3 == m.v1 - m.v2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + # No canceling terms as v3=0 is ignored thus we have a=b form + assert len(cc) == 0 + assert not k + assert not tr + + # Set v3 above zero tolerance + m.v3.set_value(1e-4) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + assert vv == [-1e-4, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + @pytest.mark.unit + def test_canceling_equality_expr_compound(self): + m = ConcreteModel() + m.v1 = Var(["a", "b"], initialize=5e6) + m.v2 = Var(initialize=1e7) + m.v3 = Var(initialize=0) + + expr = m.v3 == sum(v for v in m.v1.values()) - m.v2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + # No canceling terms as v3=0 is ignored thus we have a=b form + assert len(cc) == 0 + assert not k + assert not tr + + # Set v3 above zero tolerance + m.v3.set_value(1e-4) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + assert vv == [-1e-4, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + # Double check for a+eps=c form gets flagged in some way + @pytest.mark.unit + def test_canceling_equality_expr_canceling_sides(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=1e-8) + + expr1 = m.v1 + m.v3 + expr = m.v2 == expr1 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, pytest.approx(1 + 1e-8, abs=1e-8)] + assert expr1 in mm + assert len(mm) == 1 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + # Check to make sure simple linking constraints are not flagged as canceling + @pytest.mark.unit + def test_linking_equality_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + + expr = m.v1 == m.v2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, 1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_linking_equality_expr_compound(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=1) + + expr = m.v1 == m.v2 * m.v3 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, 1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.component + def test_external_function_w_string_argument(self): + m = ConcreteModel() + m.properties = iapws95.Iapws95ParameterBlock() + m.state = m.properties.build_state_block([0]) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.state[0].enth_mol + ) + + assert vv == [pytest.approx(1.1021387e-2, rel=1e-6)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + # Test nested external functions + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.state[0].temperature + ) + + assert vv == [pytest.approx(270.4877, rel=1e-6)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_zero_tolerance(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e-12) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=1e-12) + + expr1 = m.v1 - m.v3 + expr2 = expr1 + 1 + expr = m.v2 == expr2 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, pytest.approx(1, abs=1e-8)] + # We expect no mismatches, as smallest terms are below zero tolerance + assert len(mm) == 0 + # We expect no canceling terms, as v1 and v3 are ignored (zero tolerance), leaving + # 1 == 1 + assert len(cc) == 0 + assert not k + assert not tr + + # Set v1 and v3 above zero tolerance + m.v1.set_value(1e-8) + m.v3.set_value(1e-8) + + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, pytest.approx(1, abs=1e-8)] + assert expr2 in mm + assert len(mm) == 1 + assert expr in cc + assert len(cc) == 1 + assert not k + assert not tr + + # Test to make sure scaled constraints are not flagged as issues + @pytest.mark.unit + def test_scaled_equality_expr_product(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=2) + + expr = 0 == m.v3 * (m.v1 - m.v2) + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, 0] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_scaled_equality_expr_division(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=2) + + expr = 0 == (m.v1 - m.v2) / m.v3 + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, 0] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_mole_fraction_constraint(self): + m = ConcreteModel() + m.mole_frac_a = Var(initialize=0.5) + m.flow_a = Var(initialize=100) + m.flow_b = Var(initialize=100) + + expr = m.mole_frac_a * (m.flow_a + m.flow_b) == m.flow_a + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(-100, rel=1e-5), pytest.approx(100, rel=1e-5)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr + + @pytest.mark.unit + def test_mole_fraction_constraint_trace(self): + m = ConcreteModel() + m.mole_frac_a = Var(initialize=0.999810) + m.flow_a = Var(initialize=122.746) + m.flow_b = Var(initialize=0.0233239) + + expr = m.mole_frac_a * (m.flow_a + m.flow_b) == m.flow_a + vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [ + pytest.approx(-122.746, rel=1e-5), + pytest.approx(122.746, rel=1e-5), + ] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + assert not tr diff --git a/setup.py b/setup.py index d5303b756a..19868f5d1d 100644 --- a/setup.py +++ b/setup.py @@ -82,7 +82,7 @@ def __getitem__(self, key): # Put abstract (non-versioned) deps here. # Concrete dependencies go in requirements[-dev].txt install_requires=[ - "pyomo >= 6.8.0", + "pyomo @ https://github.com/IDAES/pyomo/archive/6.8.0.idaes.2024.10.25.zip", "pint >= 0.24.1", # required to use Pyomo units. Pint 0.24.1 needed for Python 3.9 support "networkx", # required to use Pyomo network "numpy>=1,<3",