From 80fabccf7e30e9e4e80b73d1d410fc2dd169110e Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 12 Sep 2023 09:58:17 -0600 Subject: [PATCH 01/13] diagnostic tool to look for potential evaluation errors --- idaes/core/util/model_diagnostics.py | 214 ++++++++++++++++++++++++++- 1 file changed, 213 insertions(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 92127f7c2e..411e8c8af7 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -19,7 +19,9 @@ from operator import itemgetter from sys import stdout +import math from math import log +from typing import List, Sequence import numpy as np from scipy.linalg import svd @@ -28,6 +30,7 @@ from pyomo.environ import ( Binary, + Integers, Block, check_optimal_termination, ConcreteModel, @@ -38,15 +41,28 @@ SolverFactory, value, Var, + is_fixed, +) +from pyomo.core.expr.numeric_expr import ( + DivisionExpression, + NPV_DivisionExpression, + PowExpression, + NPV_PowExpression, + UnaryFunctionExpression, + NPV_UnaryFunctionExpression, + NumericExpression, ) from pyomo.core.base.block import _BlockData +from pyomo.core.base.var import _GeneralVarData +from pyomo.repn.standard_repn import generate_standard_repn from pyomo.common.collections import ComponentSet from pyomo.common.config import ConfigDict, ConfigValue, document_kwargs_from_configdict from pyomo.util.check_units import identify_inconsistent_units from pyomo.contrib.incidence_analysis import IncidenceGraphInterface -from pyomo.core.expr.visitor import identify_variables +from pyomo.core.expr.visitor import identify_variables, StreamBasedExpressionVisitor from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.asl import AmplInterface +from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr from pyomo.common.deprecation import deprecation_warning from idaes.core.util.model_statistics import ( @@ -1045,6 +1061,202 @@ def report_numerical_issues(self, stream=stdout): footer="=", ) + def _collect_potential_eval_errors(self): + res = list() + warnings = list() + cautions = list() + for con in self.model.component_data_objects(Constraint, active=True, descend_into=True): + walker = _EvalErrorWalker() + con_warnings, con_cautions = walker.walk_expression(con.body) + for msg in con_warnings: + msg = f'{con.name}: ' + msg + warnings.append(msg) + for msg in con_cautions: + msg = f'{con.name}: ' + msg + cautions.append(msg) + for obj in self.model.component_data_objects(Objective, active=True, descend_into=True): + walker = _EvalErrorWalker() + obj_warnings, obj_cautions = walker.walk_expression(obj.expr) + for msg in obj_warnings: + msg = f'{obj.name}: ' + msg + warnings.append(msg) + for msg in obj_cautions: + msg = f'{obj.name}: ' + msg + cautions.append(msg) + + return warnings, cautions + + def report_potential_evaluation_errors(self, stream=stdout): + warnings, cautions = self._collect_potential_eval_errors() + _write_report_section( + stream=stream, + lines_list=warnings, + title=f"{len(warnings)} WARNINGS", + line_if_empty="No warnings found!", + header="=", + ) + _write_report_section( + stream=stream, + lines_list=cautions, + title=f"{len(cautions)} Cautions", + line_if_empty="No cautions found!", + footer="=", + ) + + +def _get_bounds_with_inf(node: NumericExpression): + lb, ub = compute_bounds_on_expr(node) + if lb is None: + lb = -math.inf + if ub is None: + ub = math.inf + return lb, ub + + +def _caution_expression_argument( + node: NumericExpression, + args_to_check: Sequence[NumericExpression], + caution_list: List[str] +): + should_caution = False + for arg in args_to_check: + if is_fixed(arg): + continue + if isinstance(arg, _GeneralVarData): + continue + should_caution = True + break + if should_caution: + msg = f'Potential evaluation error in {node}; ' + msg += 'arguments are expressions with bounds that are not strictly ' + msg += 'enforced; try making the argument a variable' + caution_list.append(msg) + + +def _check_eval_error_division(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, [node.args[1]], caution_list) + lb, ub = _get_bounds_with_inf(node.args[1]) + if lb <= 0 <= ub: + msg = f'Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})' + warn_list.append(msg) + + +def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + arg1, arg2 = node.args + + integer_domains = ComponentSet([Binary, Integers]) + + # if the exponent is an integer, there should not be any evaluation errors + if isinstance(arg2, _GeneralVarData) and arg2.domain in integer_domains: + # life is good. The exponent is an integer variable + return None + lb2, ub2 = _get_bounds_with_inf(arg2) + if lb2 == ub2 and lb2 == round(lb2): + # life is good. The exponent is fixed to an integer + return None + repn = generate_standard_repn(arg2, quadratic=True) + if ( + repn.nonlinear_expr is None + and repn.constant == round(repn.constant) + and all(i.domain in integer_domains for i in repn.linear_vars) + and all(i[0].domain in integer_domains for i in repn.quadratic_vars) + and all(i[1].domain in integer_domains for i in repn.quadratic_vars) + and all(i == round(i) for i in repn.linear_coefs) + and all(i == round(i) for i in repn.quadratic_coefs) + ): + # Life is good. The exponent is a linear or quadratic expression containing + # only integer variables with integer coefficients + return None + + _caution_expression_argument(node, node.args, caution_list) + + # if the base is positive, there should not be any evaluation errors + lb1, ub1 = _get_bounds_with_inf(arg1) + if lb1 > 0: + return None + if lb1 >= 0 and lb2 >= 0: + return None + + msg = f'Potential evaluation error in {node}; ' + msg += f'base bounds are ({lb1}, {ub1}); ' + msg += f'exponent bounds are ({lb2}, {ub2})' + warn_list.append(msg) + + +def _check_eval_error_log(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, node.args, caution_list) + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb <= 0: + msg = f'Potential log of a negative number in {node}; Argument bounds are ({lb}, {ub})' + warn_list.append(msg) + + +def _check_eval_error_tan(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, node.args, caution_list) + lb, ub = _get_bounds_with_inf(node) + if not (math.isfinite(lb) and math.isfinite(ub)): + msg = f'{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}' + warn_list.append(msg) + + +def _check_eval_error_asin(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, node.args, caution_list) + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < -1 or ub > 1: + msg = f'Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})' + warn_list.append(msg) + + +def _check_eval_error_acos(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, node.args, caution_list) + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < -1 or ub > 1: + msg = f'Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})' + warn_list.append(msg) + + +def _check_eval_error_sqrt(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + _caution_expression_argument(node, node.args, caution_list) + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < 0: + msg = f'Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})' + warn_list.append(msg) + + +_unary_eval_err_handler = dict() +_unary_eval_err_handler['log'] = _check_eval_error_log +_unary_eval_err_handler['log10'] = _check_eval_error_log +_unary_eval_err_handler['tan'] = _check_eval_error_tan +_unary_eval_err_handler['asin'] = _check_eval_error_asin +_unary_eval_err_handler['acos'] = _check_eval_error_acos +_unary_eval_err_handler['sqrt'] = _check_eval_error_sqrt + + +def _check_eval_error_unary(node: NumericExpression, warn_list: List[str], caution_list: List[str]): + if node.getname() in _unary_eval_err_handler: + _unary_eval_err_handler[node.getname()](node, warn_list, caution_list) + + +_eval_err_handler = dict() +_eval_err_handler[DivisionExpression] = _check_eval_error_division +_eval_err_handler[NPV_DivisionExpression] = _check_eval_error_division +_eval_err_handler[PowExpression] = _check_eval_error_pow +_eval_err_handler[NPV_PowExpression] = _check_eval_error_pow +_eval_err_handler[UnaryFunctionExpression] = _check_eval_error_unary +_eval_err_handler[NPV_UnaryFunctionExpression] = _check_eval_error_unary + + +class _EvalErrorWalker(StreamBasedExpressionVisitor): + def __init__(self): + super().__init__() + self._warn_list = list() + self._caution_list = list() + + def exitNode(self, node, data): + if type(node) in _eval_err_handler: + _eval_err_handler[type(node)](node, self._warn_list, self._caution_list) + return self._warn_list, self._caution_list + class DegeneracyHunter: """ From 8fb3b339bd285a5e1b694d7a65c54d7e953b14c2 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 11:50:27 -0600 Subject: [PATCH 02/13] working on diagnostic tool to detect possible evaluation errors --- idaes/core/util/model_diagnostics.py | 116 +++++++++++++++------------ 1 file changed, 64 insertions(+), 52 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index fc8c308d48..e3a62bab9c 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1318,6 +1318,51 @@ def report_numerical_issues(self, stream=None): footer="=", ) + def _collect_potential_eval_errors(self): + res = list() + warnings = list() + cautions = list() + for con in self.model.component_data_objects(Constraint, active=True, descend_into=True): + walker = _EvalErrorWalker() + con_warnings, con_cautions = walker.walk_expression(con.body) + for msg in con_warnings: + msg = f'{con.name}: ' + msg + warnings.append(msg) + for msg in con_cautions: + msg = f'{con.name}: ' + msg + cautions.append(msg) + for obj in self.model.component_data_objects(Objective, active=True, descend_into=True): + walker = _EvalErrorWalker() + obj_warnings, obj_cautions = walker.walk_expression(obj.expr) + for msg in obj_warnings: + msg = f'{obj.name}: ' + msg + warnings.append(msg) + for msg in obj_cautions: + msg = f'{obj.name}: ' + msg + cautions.append(msg) + + return warnings, cautions + + def report_potential_evaluation_errors(self, stream=None): + if stream is None: + stream = sys.stdout + + warnings, cautions = self._collect_potential_eval_errors() + _write_report_section( + stream=stream, + lines_list=warnings, + title=f"{len(warnings)} WARNINGS", + line_if_empty="No warnings found!", + header="=", + ) + _write_report_section( + stream=stream, + lines_list=cautions, + title=f"{len(cautions)} Cautions", + line_if_empty="No cautions found!", + footer="=", + ) + @document_kwargs_from_configdict(SVDCONFIG) def prepare_svd_toolbox(self, **kwargs): """ @@ -1637,48 +1682,6 @@ def display_variables_in_constraint(self, constraint, stream=None): footer="=", ) - def _collect_potential_eval_errors(self): - res = list() - warnings = list() - cautions = list() - for con in self.model.component_data_objects(Constraint, active=True, descend_into=True): - walker = _EvalErrorWalker() - con_warnings, con_cautions = walker.walk_expression(con.body) - for msg in con_warnings: - msg = f'{con.name}: ' + msg - warnings.append(msg) - for msg in con_cautions: - msg = f'{con.name}: ' + msg - cautions.append(msg) - for obj in self.model.component_data_objects(Objective, active=True, descend_into=True): - walker = _EvalErrorWalker() - obj_warnings, obj_cautions = walker.walk_expression(obj.expr) - for msg in obj_warnings: - msg = f'{obj.name}: ' + msg - warnings.append(msg) - for msg in obj_cautions: - msg = f'{obj.name}: ' + msg - cautions.append(msg) - - return warnings, cautions - - def report_potential_evaluation_errors(self, stream=stdout): - warnings, cautions = self._collect_potential_eval_errors() - _write_report_section( - stream=stream, - lines_list=warnings, - title=f"{len(warnings)} WARNINGS", - line_if_empty="No warnings found!", - header="=", - ) - _write_report_section( - stream=stream, - lines_list=cautions, - title=f"{len(cautions)} Cautions", - line_if_empty="No cautions found!", - footer="=", - ) - def _get_bounds_with_inf(node: NumericExpression): lb, ub = compute_bounds_on_expr(node) @@ -1705,7 +1708,7 @@ def _caution_expression_argument( if should_caution: msg = f'Potential evaluation error in {node}; ' msg += 'arguments are expressions with bounds that are not strictly ' - msg += 'enforced; try making the argument a variable' + msg += 'enforced;' caution_list.append(msg) @@ -1719,17 +1722,20 @@ def _check_eval_error_division(node: NumericExpression, warn_list: List[str], ca def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution_list: List[str]): arg1, arg2 = node.args + lb1, ub1 = _get_bounds_with_inf(arg1) + lb2, ub2 = _get_bounds_with_inf(arg2) integer_domains = ComponentSet([Binary, Integers]) + integer_exponent = False # if the exponent is an integer, there should not be any evaluation errors if isinstance(arg2, _GeneralVarData) and arg2.domain in integer_domains: - # life is good. The exponent is an integer variable - return None - lb2, ub2 = _get_bounds_with_inf(arg2) + # The exponent is an integer variable + # check if the base can be zero + integer_exponent = True if lb2 == ub2 and lb2 == round(lb2): - # life is good. The exponent is fixed to an integer - return None + # The exponent is fixed to an integer + integer_exponent = True repn = generate_standard_repn(arg2, quadratic=True) if ( repn.nonlinear_expr is None @@ -1740,17 +1746,23 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution and all(i == round(i) for i in repn.linear_coefs) and all(i == round(i) for i in repn.quadratic_coefs) ): - # Life is good. The exponent is a linear or quadratic expression containing + # The exponent is a linear or quadratic expression containing # only integer variables with integer coefficients - return None + integer_exponent = True - _caution_expression_argument(node, node.args, caution_list) + if integer_exponent and (lb1 > 0 or ub1 < 0): + # life is good; the exponent is an integer and the base is nonzero + return None + elif integer_exponent and lb2 >= 0: + # life is good; the exponent is a nonnegative integer + return None # if the base is positive, there should not be any evaluation errors - lb1, ub1 = _get_bounds_with_inf(arg1) if lb1 > 0: + _caution_expression_argument(node, node.args, caution_list) return None if lb1 >= 0 and lb2 >= 0: + _caution_expression_argument(node, node.args, caution_list) return None msg = f'Potential evaluation error in {node}; ' From 8b52c20e969c376f44ab31300b00ac5113729040 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 12:06:38 -0600 Subject: [PATCH 03/13] working on diagnostic tool to detect possible evaluation errors --- idaes/core/util/model_diagnostics.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index e3a62bab9c..90211bec9d 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1244,6 +1244,10 @@ def report_structural_issues(self, stream=None): warnings, next_steps = self._collect_structural_warnings() cautions = self._collect_structural_cautions() + eval_error_warnings, eval_error_cautions = self._collect_potential_eval_errors() + warnings.extend(eval_error_warnings) + cautions.extend(eval_error_cautions) + _write_report_section( stream=stream, lines_list=stats, title="Model Statistics", header="=" ) @@ -1343,7 +1347,7 @@ def _collect_potential_eval_errors(self): return warnings, cautions - def report_potential_evaluation_errors(self, stream=None): + def display_potential_evaluation_errors(self, stream=None): if stream is None: stream = sys.stdout @@ -1775,7 +1779,7 @@ def _check_eval_error_log(node: NumericExpression, warn_list: List[str], caution _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb <= 0: - msg = f'Potential log of a negative number in {node}; Argument bounds are ({lb}, {ub})' + msg = f'Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})' warn_list.append(msg) From 048dfd31f9d9db3f9247c9df472ab438fa0e418c Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 12:14:10 -0600 Subject: [PATCH 04/13] run black --- idaes/core/util/model_diagnostics.py | 104 ++++++++++++++++----------- 1 file changed, 62 insertions(+), 42 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 90211bec9d..36416ab472 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1326,27 +1326,31 @@ def _collect_potential_eval_errors(self): res = list() warnings = list() cautions = list() - for con in self.model.component_data_objects(Constraint, active=True, descend_into=True): + for con in self.model.component_data_objects( + Constraint, active=True, descend_into=True + ): walker = _EvalErrorWalker() con_warnings, con_cautions = walker.walk_expression(con.body) for msg in con_warnings: - msg = f'{con.name}: ' + msg + msg = f"{con.name}: " + msg warnings.append(msg) for msg in con_cautions: - msg = f'{con.name}: ' + msg + msg = f"{con.name}: " + msg cautions.append(msg) - for obj in self.model.component_data_objects(Objective, active=True, descend_into=True): + for obj in self.model.component_data_objects( + Objective, active=True, descend_into=True + ): walker = _EvalErrorWalker() obj_warnings, obj_cautions = walker.walk_expression(obj.expr) for msg in obj_warnings: - msg = f'{obj.name}: ' + msg + msg = f"{obj.name}: " + msg warnings.append(msg) for msg in obj_cautions: - msg = f'{obj.name}: ' + msg + msg = f"{obj.name}: " + msg cautions.append(msg) return warnings, cautions - + def display_potential_evaluation_errors(self, stream=None): if stream is None: stream = sys.stdout @@ -1699,7 +1703,7 @@ def _get_bounds_with_inf(node: NumericExpression): def _caution_expression_argument( node: NumericExpression, args_to_check: Sequence[NumericExpression], - caution_list: List[str] + caution_list: List[str], ): should_caution = False for arg in args_to_check: @@ -1710,21 +1714,25 @@ def _caution_expression_argument( should_caution = True break if should_caution: - msg = f'Potential evaluation error in {node}; ' - msg += 'arguments are expressions with bounds that are not strictly ' - msg += 'enforced;' + msg = f"Potential evaluation error in {node}; " + msg += "arguments are expressions with bounds that are not strictly " + msg += "enforced;" caution_list.append(msg) -def _check_eval_error_division(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_division( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, [node.args[1]], caution_list) lb, ub = _get_bounds_with_inf(node.args[1]) if lb <= 0 <= ub: - msg = f'Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})' + msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_pow( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): arg1, arg2 = node.args lb1, ub1 = _get_bounds_with_inf(arg1) lb2, ub2 = _get_bounds_with_inf(arg2) @@ -1742,13 +1750,13 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution integer_exponent = True repn = generate_standard_repn(arg2, quadratic=True) if ( - repn.nonlinear_expr is None - and repn.constant == round(repn.constant) - and all(i.domain in integer_domains for i in repn.linear_vars) - and all(i[0].domain in integer_domains for i in repn.quadratic_vars) - and all(i[1].domain in integer_domains for i in repn.quadratic_vars) - and all(i == round(i) for i in repn.linear_coefs) - and all(i == round(i) for i in repn.quadratic_coefs) + repn.nonlinear_expr is None + and repn.constant == round(repn.constant) + and all(i.domain in integer_domains for i in repn.linear_vars) + and all(i[0].domain in integer_domains for i in repn.quadratic_vars) + and all(i[1].domain in integer_domains for i in repn.quadratic_vars) + and all(i == round(i) for i in repn.linear_coefs) + and all(i == round(i) for i in repn.quadratic_coefs) ): # The exponent is a linear or quadratic expression containing # only integer variables with integer coefficients @@ -1769,62 +1777,74 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str], caution _caution_expression_argument(node, node.args, caution_list) return None - msg = f'Potential evaluation error in {node}; ' - msg += f'base bounds are ({lb1}, {ub1}); ' - msg += f'exponent bounds are ({lb2}, {ub2})' + msg = f"Potential evaluation error in {node}; " + msg += f"base bounds are ({lb1}, {ub1}); " + msg += f"exponent bounds are ({lb2}, {ub2})" warn_list.append(msg) -def _check_eval_error_log(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_log( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb <= 0: - msg = f'Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})' + msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_tan(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_tan( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node) if not (math.isfinite(lb) and math.isfinite(ub)): - msg = f'{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}' + msg = f"{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}" warn_list.append(msg) -def _check_eval_error_asin(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_asin( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: - msg = f'Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})' + msg = f"Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_acos(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_acos( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: - msg = f'Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})' + msg = f"Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_sqrt(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_sqrt( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < 0: - msg = f'Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})' + msg = f"Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) _unary_eval_err_handler = dict() -_unary_eval_err_handler['log'] = _check_eval_error_log -_unary_eval_err_handler['log10'] = _check_eval_error_log -_unary_eval_err_handler['tan'] = _check_eval_error_tan -_unary_eval_err_handler['asin'] = _check_eval_error_asin -_unary_eval_err_handler['acos'] = _check_eval_error_acos -_unary_eval_err_handler['sqrt'] = _check_eval_error_sqrt +_unary_eval_err_handler["log"] = _check_eval_error_log +_unary_eval_err_handler["log10"] = _check_eval_error_log +_unary_eval_err_handler["tan"] = _check_eval_error_tan +_unary_eval_err_handler["asin"] = _check_eval_error_asin +_unary_eval_err_handler["acos"] = _check_eval_error_acos +_unary_eval_err_handler["sqrt"] = _check_eval_error_sqrt -def _check_eval_error_unary(node: NumericExpression, warn_list: List[str], caution_list: List[str]): +def _check_eval_error_unary( + node: NumericExpression, warn_list: List[str], caution_list: List[str] +): if node.getname() in _unary_eval_err_handler: _unary_eval_err_handler[node.getname()](node, warn_list, caution_list) @@ -1843,7 +1863,7 @@ def __init__(self): super().__init__() self._warn_list = list() self._caution_list = list() - + def exitNode(self, node, data): if type(node) in _eval_err_handler: _eval_err_handler[type(node)](node, self._warn_list, self._caution_list) From 09bc21b6ad0ae07d342f4a9c141b8943e4d17511 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 12:19:26 -0600 Subject: [PATCH 05/13] working on tests for evaluation error detection --- idaes/core/util/tests/test_model_diagnostics.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index e91a69160b..edff4a6f45 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -2579,3 +2579,15 @@ def test_ipopt_solve_halt_on_error(capsys): captured = capsys.readouterr() assert "c: can't evaluate log(-5)." in captured.out + + +@pytest.mark.unit +def test_eval_error_detection_div(): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == 1/m.x) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + assert len(warnings) == 0 + assert len(cautions) == 0 From 6f3cae0680737167ab8722270f4d26884db1d51d Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 13:27:35 -0600 Subject: [PATCH 06/13] working on tests for evaluation error detection --- idaes/core/util/model_diagnostics.py | 18 +- .../core/util/tests/test_model_diagnostics.py | 241 +++++++++++++++++- 2 files changed, 243 insertions(+), 16 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 36416ab472..074415d363 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1723,11 +1723,12 @@ def _caution_expression_argument( def _check_eval_error_division( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, [node.args[1]], caution_list) lb, ub = _get_bounds_with_inf(node.args[1]) if lb <= 0 <= ub: msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) + else: + _caution_expression_argument(node, [node.args[1]], caution_list) def _check_eval_error_pow( @@ -1786,51 +1787,56 @@ def _check_eval_error_pow( def _check_eval_error_log( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb <= 0: msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) + else: + _caution_expression_argument(node, node.args, caution_list) def _check_eval_error_tan( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node) if not (math.isfinite(lb) and math.isfinite(ub)): msg = f"{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}" warn_list.append(msg) + else: + _caution_expression_argument(node, node.args, caution_list) def _check_eval_error_asin( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) + else: + _caution_expression_argument(node, node.args, caution_list) def _check_eval_error_acos( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) + else: + _caution_expression_argument(node, node.args, caution_list) def _check_eval_error_sqrt( node: NumericExpression, warn_list: List[str], caution_list: List[str] ): - _caution_expression_argument(node, node.args, caution_list) lb, ub = _get_bounds_with_inf(node.args[0]) if lb < 0: msg = f"Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) + else: + _caution_expression_argument(node, node.args, caution_list) _unary_eval_err_handler = dict() diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index edff4a6f45..df4459e956 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -14,6 +14,7 @@ This module contains model diagnostic utility functions for use in IDAES (Pyomo) models. """ from io import StringIO +import math import numpy as np import pytest @@ -23,6 +24,10 @@ Constraint, Expression, log, + tan, + asin, + acos, + sqrt, Objective, Set, SolverFactory, @@ -31,6 +36,8 @@ units, value, Var, + Param, + Integers, ) from pyomo.common.collections import ComponentSet from pyomo.contrib.pynumero.asl import AmplInterface @@ -41,6 +48,7 @@ from idaes.core.solvers import get_solver from idaes.core import FlowsheetBlock from idaes.core.util.testing import PhysicalParameterTestBlock +from unittest import TestCase # TODO: Add pyomo.dae test case """ @@ -2581,13 +2589,226 @@ def test_ipopt_solve_halt_on_error(capsys): assert "c: can't evaluate log(-5)." in captured.out -@pytest.mark.unit -def test_eval_error_detection_div(): - m = ConcreteModel() - m.x = Var(bounds=(1, None)) - m.y = Var() - m.c = Constraint(expr=m.y == 1/m.x) - dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() - assert len(warnings) == 0 - assert len(cautions) == 0 +class TestEvalErrorDetection(TestCase): + @pytest.mark.unit + def test_div(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == 1/m.x) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(-1) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential division by 0 in 1/x; Denominator bounds are (-1, inf)') + + @pytest.mark.unit + def test_pow1(self): + m = ConcreteModel() + m.x = Var(bounds=(None, None)) + m.y = Var() + m.p = Param(initialize=2, mutable=True) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.p.value = 2.5 + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)') + + m.x.setlb(1) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + @pytest.mark.unit + def test_pow2(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.p = Var(domain=Integers) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(None) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (-inf, inf)') + + @pytest.mark.unit + def test_pow3(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.p = Var(bounds=(0, None)) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + @pytest.mark.unit + def test_log(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == log(m.x)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(-1) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential log of a non-positive number in log(x); Argument bounds are (-1, inf)') + + @pytest.mark.unit + def test_tan(self): + m = ConcreteModel() + m.x = Var(bounds=(-math.pi/4, math.pi/4)) + m.y = Var() + m.c = Constraint(expr=m.y == tan(m.x)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(-math.pi) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: tan(x) may evaluate to -inf or inf; Argument bounds are (-3.141592653589793, 0.7853981633974483)') + + @pytest.mark.unit + def test_asin(self): + m = ConcreteModel() + m.x = Var(bounds=(-0.5, 0.5)) + m.y = Var() + m.c = Constraint(expr=m.y == asin(m.x)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(None) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-inf, 0.5)') + + m.x.setlb(-0.5) + m.x.setub(None) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-0.5, inf)') + + @pytest.mark.unit + def test_acos(self): + m = ConcreteModel() + m.x = Var(bounds=(-0.5, 0.5)) + m.y = Var() + m.c = Constraint(expr=m.y == acos(m.x)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(None) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-inf, 0.5)') + + m.x.setlb(-0.5) + m.x.setub(None) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + self.assertEqual(w, 'c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-0.5, inf)') + + @pytest.mark.unit + def test_sqrt(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == sqrt(m.x)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 0) + + m.x.setlb(-1) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + print(w) + self.assertEqual(w, 'c: Potential square root of a negative number in sqrt(x); Argument bounds are (-1, inf)') + + @pytest.mark.unit + def test_expression(self): + m = ConcreteModel() + m.x = Var(bounds=(3, None)) + m.y = Var() + m.c = Constraint(expr=m.y == log(m.x - 1)) + dtb = DiagnosticsToolbox(m) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + self.assertEqual(len(cautions), 1) + w = cautions[0] + self.assertEqual(w, 'c: Potential evaluation error in log(x - 1); arguments are expressions with bounds that are not strictly enforced;') + + m.x.setlb(1) + warnings, cautions = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + self.assertEqual(len(cautions), 0) + w = warnings[0] + print(w) + self.assertEqual(w, 'c: Potential log of a non-positive number in log(x - 1); Argument bounds are (0, inf)') + + @pytest.mark.unit + def test_display(self): + stream = StringIO() + m = ConcreteModel() + m.x = Var() + m.y = Var() + m.obj = Objective(expr=m.x**2 + m.y**2.5) + m.c1 = Constraint(expr=m.y >= log(m.x)) + m.c2 = Constraint(expr=m.y >= (m.x - 1)**2.5) + m.c3 = Constraint(expr=m.x - 1 >= 0) + dtb = DiagnosticsToolbox(m) + dtb.display_potential_evaluation_errors(stream=stream) + expected = "====================================================================================\n3 WARNINGS\n\n c1: Potential log of a non-positive number in log(x); Argument bounds are (-inf, inf)\n c2: Potential evaluation error in (x - 1)**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n obj: Potential evaluation error in y**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n\n------------------------------------------------------------------------------------\n0 Cautions\n\n No cautions found!\n\n====================================================================================\n" + got = stream.getvalue() + exp_list = expected.split('\n') + got_list = got.split('\n') + print(exp_list) + print(got_list) + self.assertEqual(len(exp_list), len(got_list)) + for _exp, _got in zip(exp_list, got_list): + self.assertEqual(_exp, _got) From 2a135a1f350e4bbab40190a3501f3db9b93b02f1 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Mon, 30 Oct 2023 13:28:46 -0600 Subject: [PATCH 07/13] run black --- .../core/util/tests/test_model_diagnostics.py | 69 ++++++++++++++----- 1 file changed, 52 insertions(+), 17 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index df4459e956..fb9a1fa7b4 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -2595,7 +2595,7 @@ def test_div(self): m = ConcreteModel() m.x = Var(bounds=(1, None)) m.y = Var() - m.c = Constraint(expr=m.y == 1/m.x) + m.c = Constraint(expr=m.y == 1 / m.x) dtb = DiagnosticsToolbox(m) warnings, cautions = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) @@ -2606,7 +2606,9 @@ def test_div(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential division by 0 in 1/x; Denominator bounds are (-1, inf)') + self.assertEqual( + w, "c: Potential division by 0 in 1/x; Denominator bounds are (-1, inf)" + ) @pytest.mark.unit def test_pow1(self): @@ -2625,7 +2627,10 @@ def test_pow1(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)') + self.assertEqual( + w, + "c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)", + ) m.x.setlb(1) warnings, cautions = dtb._collect_potential_eval_errors() @@ -2649,7 +2654,10 @@ def test_pow2(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (-inf, inf)') + self.assertEqual( + w, + "c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (-inf, inf)", + ) @pytest.mark.unit def test_pow3(self): @@ -2679,12 +2687,15 @@ def test_log(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential log of a non-positive number in log(x); Argument bounds are (-1, inf)') + self.assertEqual( + w, + "c: Potential log of a non-positive number in log(x); Argument bounds are (-1, inf)", + ) @pytest.mark.unit def test_tan(self): m = ConcreteModel() - m.x = Var(bounds=(-math.pi/4, math.pi/4)) + m.x = Var(bounds=(-math.pi / 4, math.pi / 4)) m.y = Var() m.c = Constraint(expr=m.y == tan(m.x)) dtb = DiagnosticsToolbox(m) @@ -2697,7 +2708,10 @@ def test_tan(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: tan(x) may evaluate to -inf or inf; Argument bounds are (-3.141592653589793, 0.7853981633974483)') + self.assertEqual( + w, + "c: tan(x) may evaluate to -inf or inf; Argument bounds are (-3.141592653589793, 0.7853981633974483)", + ) @pytest.mark.unit def test_asin(self): @@ -2715,7 +2729,10 @@ def test_asin(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-inf, 0.5)') + self.assertEqual( + w, + "c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-inf, 0.5)", + ) m.x.setlb(-0.5) m.x.setub(None) @@ -2723,7 +2740,10 @@ def test_asin(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-0.5, inf)') + self.assertEqual( + w, + "c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-0.5, inf)", + ) @pytest.mark.unit def test_acos(self): @@ -2741,7 +2761,10 @@ def test_acos(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-inf, 0.5)') + self.assertEqual( + w, + "c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-inf, 0.5)", + ) m.x.setlb(-0.5) m.x.setub(None) @@ -2749,7 +2772,10 @@ def test_acos(self): self.assertEqual(len(warnings), 1) self.assertEqual(len(cautions), 0) w = warnings[0] - self.assertEqual(w, 'c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-0.5, inf)') + self.assertEqual( + w, + "c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-0.5, inf)", + ) @pytest.mark.unit def test_sqrt(self): @@ -2768,7 +2794,10 @@ def test_sqrt(self): self.assertEqual(len(cautions), 0) w = warnings[0] print(w) - self.assertEqual(w, 'c: Potential square root of a negative number in sqrt(x); Argument bounds are (-1, inf)') + self.assertEqual( + w, + "c: Potential square root of a negative number in sqrt(x); Argument bounds are (-1, inf)", + ) @pytest.mark.unit def test_expression(self): @@ -2781,7 +2810,10 @@ def test_expression(self): self.assertEqual(len(warnings), 0) self.assertEqual(len(cautions), 1) w = cautions[0] - self.assertEqual(w, 'c: Potential evaluation error in log(x - 1); arguments are expressions with bounds that are not strictly enforced;') + self.assertEqual( + w, + "c: Potential evaluation error in log(x - 1); arguments are expressions with bounds that are not strictly enforced;", + ) m.x.setlb(1) warnings, cautions = dtb._collect_potential_eval_errors() @@ -2789,7 +2821,10 @@ def test_expression(self): self.assertEqual(len(cautions), 0) w = warnings[0] print(w) - self.assertEqual(w, 'c: Potential log of a non-positive number in log(x - 1); Argument bounds are (0, inf)') + self.assertEqual( + w, + "c: Potential log of a non-positive number in log(x - 1); Argument bounds are (0, inf)", + ) @pytest.mark.unit def test_display(self): @@ -2799,14 +2834,14 @@ def test_display(self): m.y = Var() m.obj = Objective(expr=m.x**2 + m.y**2.5) m.c1 = Constraint(expr=m.y >= log(m.x)) - m.c2 = Constraint(expr=m.y >= (m.x - 1)**2.5) + m.c2 = Constraint(expr=m.y >= (m.x - 1) ** 2.5) m.c3 = Constraint(expr=m.x - 1 >= 0) dtb = DiagnosticsToolbox(m) dtb.display_potential_evaluation_errors(stream=stream) expected = "====================================================================================\n3 WARNINGS\n\n c1: Potential log of a non-positive number in log(x); Argument bounds are (-inf, inf)\n c2: Potential evaluation error in (x - 1)**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n obj: Potential evaluation error in y**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n\n------------------------------------------------------------------------------------\n0 Cautions\n\n No cautions found!\n\n====================================================================================\n" got = stream.getvalue() - exp_list = expected.split('\n') - got_list = got.split('\n') + exp_list = expected.split("\n") + got_list = got.split("\n") print(exp_list) print(got_list) self.assertEqual(len(exp_list), len(got_list)) From e9e8fc8fc5ffd97e407168e0ab53b40c1cc41bf0 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 31 Oct 2023 06:48:47 -0600 Subject: [PATCH 08/13] evaluation error detection --- idaes/core/util/model_diagnostics.py | 107 ++++-------------- .../core/util/tests/test_model_diagnostics.py | 92 ++++----------- 2 files changed, 44 insertions(+), 155 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 074415d363..5d20656119 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -951,6 +951,13 @@ def _collect_structural_warnings(self): if any(len(x) > 0 for x in [oc_var, oc_con]): next_steps.append(self.display_overconstrained_set.__name__ + "()") + eval_warnings = self._collect_potential_eval_errors() + if len(eval_warnings) > 0: + warnings.append( + f"WARNING: Found {len(eval_warnings)} potential evaluation errors." + ) + next_steps.append(self.display_potential_evaluation_errors.__name__ + "()") + return warnings, next_steps def _collect_structural_cautions(self): @@ -1244,10 +1251,6 @@ def report_structural_issues(self, stream=None): warnings, next_steps = self._collect_structural_warnings() cautions = self._collect_structural_cautions() - eval_error_warnings, eval_error_cautions = self._collect_potential_eval_errors() - warnings.extend(eval_error_warnings) - cautions.extend(eval_error_cautions) - _write_report_section( stream=stream, lines_list=stats, title="Model Statistics", header="=" ) @@ -1322,52 +1325,39 @@ def report_numerical_issues(self, stream=None): footer="=", ) - def _collect_potential_eval_errors(self): + def _collect_potential_eval_errors(self) -> List[str]: res = list() warnings = list() - cautions = list() for con in self.model.component_data_objects( Constraint, active=True, descend_into=True ): walker = _EvalErrorWalker() - con_warnings, con_cautions = walker.walk_expression(con.body) + con_warnings = walker.walk_expression(con.body) for msg in con_warnings: msg = f"{con.name}: " + msg warnings.append(msg) - for msg in con_cautions: - msg = f"{con.name}: " + msg - cautions.append(msg) for obj in self.model.component_data_objects( Objective, active=True, descend_into=True ): walker = _EvalErrorWalker() - obj_warnings, obj_cautions = walker.walk_expression(obj.expr) + obj_warnings = walker.walk_expression(obj.expr) for msg in obj_warnings: msg = f"{obj.name}: " + msg warnings.append(msg) - for msg in obj_cautions: - msg = f"{obj.name}: " + msg - cautions.append(msg) - return warnings, cautions + return warnings def display_potential_evaluation_errors(self, stream=None): if stream is None: stream = sys.stdout - warnings, cautions = self._collect_potential_eval_errors() + warnings = self._collect_potential_eval_errors() _write_report_section( stream=stream, lines_list=warnings, title=f"{len(warnings)} WARNINGS", line_if_empty="No warnings found!", header="=", - ) - _write_report_section( - stream=stream, - lines_list=cautions, - title=f"{len(cautions)} Cautions", - line_if_empty="No cautions found!", footer="=", ) @@ -1700,40 +1690,14 @@ def _get_bounds_with_inf(node: NumericExpression): return lb, ub -def _caution_expression_argument( - node: NumericExpression, - args_to_check: Sequence[NumericExpression], - caution_list: List[str], -): - should_caution = False - for arg in args_to_check: - if is_fixed(arg): - continue - if isinstance(arg, _GeneralVarData): - continue - should_caution = True - break - if should_caution: - msg = f"Potential evaluation error in {node}; " - msg += "arguments are expressions with bounds that are not strictly " - msg += "enforced;" - caution_list.append(msg) - - -def _check_eval_error_division( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_division(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node.args[1]) if lb <= 0 <= ub: msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) - else: - _caution_expression_argument(node, [node.args[1]], caution_list) -def _check_eval_error_pow( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_pow(node: NumericExpression, warn_list: List[str]): arg1, arg2 = node.args lb1, ub1 = _get_bounds_with_inf(arg1) lb2, ub2 = _get_bounds_with_inf(arg2) @@ -1772,10 +1736,8 @@ def _check_eval_error_pow( # if the base is positive, there should not be any evaluation errors if lb1 > 0: - _caution_expression_argument(node, node.args, caution_list) return None if lb1 >= 0 and lb2 >= 0: - _caution_expression_argument(node, node.args, caution_list) return None msg = f"Potential evaluation error in {node}; " @@ -1784,59 +1746,39 @@ def _check_eval_error_pow( warn_list.append(msg) -def _check_eval_error_log( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_log(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node.args[0]) if lb <= 0: msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) - else: - _caution_expression_argument(node, node.args, caution_list) -def _check_eval_error_tan( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_tan(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node) if not (math.isfinite(lb) and math.isfinite(ub)): msg = f"{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}" warn_list.append(msg) - else: - _caution_expression_argument(node, node.args, caution_list) -def _check_eval_error_asin( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_asin(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) - else: - _caution_expression_argument(node, node.args, caution_list) -def _check_eval_error_acos( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_acos(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) - else: - _caution_expression_argument(node, node.args, caution_list) -def _check_eval_error_sqrt( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_sqrt(node: NumericExpression, warn_list: List[str]): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < 0: msg = f"Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) - else: - _caution_expression_argument(node, node.args, caution_list) _unary_eval_err_handler = dict() @@ -1848,11 +1790,9 @@ def _check_eval_error_sqrt( _unary_eval_err_handler["sqrt"] = _check_eval_error_sqrt -def _check_eval_error_unary( - node: NumericExpression, warn_list: List[str], caution_list: List[str] -): +def _check_eval_error_unary(node: NumericExpression, warn_list: List[str]): if node.getname() in _unary_eval_err_handler: - _unary_eval_err_handler[node.getname()](node, warn_list, caution_list) + _unary_eval_err_handler[node.getname()](node, warn_list) _eval_err_handler = dict() @@ -1868,12 +1808,11 @@ class _EvalErrorWalker(StreamBasedExpressionVisitor): def __init__(self): super().__init__() self._warn_list = list() - self._caution_list = list() def exitNode(self, node, data): if type(node) in _eval_err_handler: - _eval_err_handler[type(node)](node, self._warn_list, self._caution_list) - return self._warn_list, self._caution_list + _eval_err_handler[type(node)](node, self._warn_list) + return self._warn_list # TODO: Rename and redirect once old DegeneracyHunter is removed diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index fb9a1fa7b4..663136c194 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -2597,14 +2597,12 @@ def test_div(self): m.y = Var() m.c = Constraint(expr=m.y == 1 / m.x) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(-1) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, "c: Potential division by 0 in 1/x; Denominator bounds are (-1, inf)" @@ -2618,14 +2616,12 @@ def test_pow1(self): m.p = Param(initialize=2, mutable=True) m.c = Constraint(expr=m.y == m.x**m.p) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.p.value = 2.5 - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2633,9 +2629,8 @@ def test_pow1(self): ) m.x.setlb(1) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) @pytest.mark.unit def test_pow2(self): @@ -2645,14 +2640,12 @@ def test_pow2(self): m.p = Var(domain=Integers) m.c = Constraint(expr=m.y == m.x**m.p) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(None) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2667,9 +2660,8 @@ def test_pow3(self): m.p = Var(bounds=(0, None)) m.c = Constraint(expr=m.y == m.x**m.p) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) @pytest.mark.unit def test_log(self): @@ -2678,14 +2670,12 @@ def test_log(self): m.y = Var() m.c = Constraint(expr=m.y == log(m.x)) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(-1) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2699,14 +2689,12 @@ def test_tan(self): m.y = Var() m.c = Constraint(expr=m.y == tan(m.x)) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(-math.pi) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2720,14 +2708,12 @@ def test_asin(self): m.y = Var() m.c = Constraint(expr=m.y == asin(m.x)) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(None) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2736,9 +2722,8 @@ def test_asin(self): m.x.setlb(-0.5) m.x.setub(None) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2752,14 +2737,12 @@ def test_acos(self): m.y = Var() m.c = Constraint(expr=m.y == acos(m.x)) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(None) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2768,9 +2751,8 @@ def test_acos(self): m.x.setlb(-0.5) m.x.setub(None) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] self.assertEqual( w, @@ -2784,48 +2766,18 @@ def test_sqrt(self): m.y = Var() m.c = Constraint(expr=m.y == sqrt(m.x)) dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 0) m.x.setlb(-1) - warnings, cautions = dtb._collect_potential_eval_errors() + warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) w = warnings[0] - print(w) self.assertEqual( w, "c: Potential square root of a negative number in sqrt(x); Argument bounds are (-1, inf)", ) - @pytest.mark.unit - def test_expression(self): - m = ConcreteModel() - m.x = Var(bounds=(3, None)) - m.y = Var() - m.c = Constraint(expr=m.y == log(m.x - 1)) - dtb = DiagnosticsToolbox(m) - warnings, cautions = dtb._collect_potential_eval_errors() - self.assertEqual(len(warnings), 0) - self.assertEqual(len(cautions), 1) - w = cautions[0] - self.assertEqual( - w, - "c: Potential evaluation error in log(x - 1); arguments are expressions with bounds that are not strictly enforced;", - ) - - m.x.setlb(1) - warnings, cautions = dtb._collect_potential_eval_errors() - self.assertEqual(len(warnings), 1) - self.assertEqual(len(cautions), 0) - w = warnings[0] - print(w) - self.assertEqual( - w, - "c: Potential log of a non-positive number in log(x - 1); Argument bounds are (0, inf)", - ) - @pytest.mark.unit def test_display(self): stream = StringIO() @@ -2838,12 +2790,10 @@ def test_display(self): m.c3 = Constraint(expr=m.x - 1 >= 0) dtb = DiagnosticsToolbox(m) dtb.display_potential_evaluation_errors(stream=stream) - expected = "====================================================================================\n3 WARNINGS\n\n c1: Potential log of a non-positive number in log(x); Argument bounds are (-inf, inf)\n c2: Potential evaluation error in (x - 1)**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n obj: Potential evaluation error in y**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n\n------------------------------------------------------------------------------------\n0 Cautions\n\n No cautions found!\n\n====================================================================================\n" + expected = "====================================================================================\n3 WARNINGS\n\n c1: Potential log of a non-positive number in log(x); Argument bounds are (-inf, inf)\n c2: Potential evaluation error in (x - 1)**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n obj: Potential evaluation error in y**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n\n====================================================================================\n" got = stream.getvalue() exp_list = expected.split("\n") got_list = got.split("\n") - print(exp_list) - print(got_list) self.assertEqual(len(exp_list), len(got_list)) for _exp, _got in zip(exp_list, got_list): self.assertEqual(_exp, _got) From 246b92a72db38a4fbcb3a749409e6be37a2669b5 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 31 Oct 2023 07:39:59 -0600 Subject: [PATCH 09/13] evaluation error detection --- idaes/core/util/model_diagnostics.py | 62 ++++++++++---- .../core/util/tests/test_model_diagnostics.py | 81 +++++++++++++++++++ 2 files changed, 126 insertions(+), 17 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 5d20656119..27bfa68a60 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -266,6 +266,14 @@ def svd_sparse(jacobian, number_singular_values): description="Tolerance for raising a warning for small Jacobian values.", ), ) +CONFIG.declare( + "strict_evaluation_error_detection", + ConfigValue( + default=True, + domain=bool, + description="If False, warnings will not be generated for things like log(x) with x >= 0", + ), +) SVDCONFIG = ConfigDict() @@ -1331,7 +1339,7 @@ def _collect_potential_eval_errors(self) -> List[str]: for con in self.model.component_data_objects( Constraint, active=True, descend_into=True ): - walker = _EvalErrorWalker() + walker = _EvalErrorWalker(self.config) con_warnings = walker.walk_expression(con.body) for msg in con_warnings: msg = f"{con.name}: " + msg @@ -1339,7 +1347,7 @@ def _collect_potential_eval_errors(self) -> List[str]: for obj in self.model.component_data_objects( Objective, active=True, descend_into=True ): - walker = _EvalErrorWalker() + walker = _EvalErrorWalker(self.config) obj_warnings = walker.walk_expression(obj.expr) for msg in obj_warnings: msg = f"{obj.name}: " + msg @@ -1690,14 +1698,18 @@ def _get_bounds_with_inf(node: NumericExpression): return lb, ub -def _check_eval_error_division(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_division( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node.args[1]) - if lb <= 0 <= ub: + if (config.strict_evaluation_error_detection and (lb <= 0 <= ub)) or (lb < 0 < ub): msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_pow(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_pow( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): arg1, arg2 = node.args lb1, ub1 = _get_bounds_with_inf(arg1) lb2, ub2 = _get_bounds_with_inf(arg2) @@ -1727,7 +1739,10 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str]): # only integer variables with integer coefficients integer_exponent = True - if integer_exponent and (lb1 > 0 or ub1 < 0): + if integer_exponent and ( + (lb1 > 0 or ub1 < 0) + or (not config.strict_evaluation_error_detection and (lb1 >= 0 or ub1 <= 0)) + ): # life is good; the exponent is an integer and the base is nonzero return None elif integer_exponent and lb2 >= 0: @@ -1735,7 +1750,7 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str]): return None # if the base is positive, there should not be any evaluation errors - if lb1 > 0: + if lb1 > 0 or (not config.strict_evaluation_error_detection and lb1 >= 0): return None if lb1 >= 0 and lb2 >= 0: return None @@ -1746,35 +1761,45 @@ def _check_eval_error_pow(node: NumericExpression, warn_list: List[str]): warn_list.append(msg) -def _check_eval_error_log(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_log( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node.args[0]) - if lb <= 0: + if (config.strict_evaluation_error_detection and lb <= 0) or lb < 0: msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_tan(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_tan( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node) if not (math.isfinite(lb) and math.isfinite(ub)): msg = f"{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}" warn_list.append(msg) -def _check_eval_error_asin(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_asin( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_acos(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_acos( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < -1 or ub > 1: msg = f"Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) -def _check_eval_error_sqrt(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_sqrt( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): lb, ub = _get_bounds_with_inf(node.args[0]) if lb < 0: msg = f"Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})" @@ -1790,9 +1815,11 @@ def _check_eval_error_sqrt(node: NumericExpression, warn_list: List[str]): _unary_eval_err_handler["sqrt"] = _check_eval_error_sqrt -def _check_eval_error_unary(node: NumericExpression, warn_list: List[str]): +def _check_eval_error_unary( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): if node.getname() in _unary_eval_err_handler: - _unary_eval_err_handler[node.getname()](node, warn_list) + _unary_eval_err_handler[node.getname()](node, warn_list, config) _eval_err_handler = dict() @@ -1805,13 +1832,14 @@ def _check_eval_error_unary(node: NumericExpression, warn_list: List[str]): class _EvalErrorWalker(StreamBasedExpressionVisitor): - def __init__(self): + def __init__(self, config: ConfigDict): super().__init__() self._warn_list = list() + self._config = config def exitNode(self, node, data): if type(node) in _eval_err_handler: - _eval_err_handler[type(node)](node, self._warn_list) + _eval_err_handler[type(node)](node, self._warn_list, self._config) return self._warn_list diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 663136c194..b86f35dd1b 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -2600,6 +2600,18 @@ def test_div(self): warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) + m.x.setlb(0) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, "c: Potential division by 0 in 1/x; Denominator bounds are (0, inf)" + ) + + dtb.config.strict_evaluation_error_detection = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + m.x.setlb(-1) warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) @@ -2663,6 +2675,62 @@ def test_pow3(self): warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) + @pytest.mark.unit + def test_pow4(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.c = Constraint(expr=m.y == m.x ** (-2)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2; base bounds are (0, inf); exponent bounds are (-2, -2)", + ) + + dtb.config.strict_evaluation_error_detection = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2; base bounds are (-1, inf); exponent bounds are (-2, -2)", + ) + + @pytest.mark.unit + def test_pow5(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.c = Constraint(expr=m.y == m.x ** (-2.5)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2.5; base bounds are (0, inf); exponent bounds are (-2.5, -2.5)", + ) + + dtb.config.strict_evaluation_error_detection = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2.5; base bounds are (-1, inf); exponent bounds are (-2.5, -2.5)", + ) + @pytest.mark.unit def test_log(self): m = ConcreteModel() @@ -2673,6 +2741,19 @@ def test_log(self): warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) + m.x.setlb(0) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential log of a non-positive number in log(x); Argument bounds are (0, inf)", + ) + + dtb.config.strict_evaluation_error_detection = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + m.x.setlb(-1) warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 1) From df094c5cc31b6ef27f308050bdc5ce8110304747 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 14 Nov 2023 11:38:44 -0700 Subject: [PATCH 10/13] fix pylint issues --- idaes/core/util/model_diagnostics.py | 34 ++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 27bfa68a60..ca1d2e2f31 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -22,7 +22,7 @@ from inspect import signature import math from math import log -from typing import List, Sequence +from typing import List import numpy as np from scipy.linalg import svd @@ -42,7 +42,6 @@ SolverFactory, value, Var, - is_fixed, ) from pyomo.core.expr.numeric_expr import ( DivisionExpression, @@ -56,7 +55,7 @@ from pyomo.core.base.block import _BlockData from pyomo.core.base.var import _GeneralVarData, _VarData from pyomo.core.base.constraint import _ConstraintData -from pyomo.repn.standard_repn import generate_standard_repn +from pyomo.repn.standard_repn import generate_standard_repn # pylint: disable=no-name-in-module from pyomo.common.collections import ComponentSet from pyomo.common.config import ( ConfigDict, @@ -267,7 +266,7 @@ def svd_sparse(jacobian, number_singular_values): ), ) CONFIG.declare( - "strict_evaluation_error_detection", + "warn_for_evaluation_error_at_bounds", ConfigValue( default=True, domain=bool, @@ -1334,7 +1333,6 @@ def report_numerical_issues(self, stream=None): ) def _collect_potential_eval_errors(self) -> List[str]: - res = list() warnings = list() for con in self.model.component_data_objects( Constraint, active=True, descend_into=True @@ -1356,6 +1354,16 @@ def _collect_potential_eval_errors(self) -> List[str]: return warnings def display_potential_evaluation_errors(self, stream=None): + """ + Prints constraints that may be prone to evaluation errors + (e.g., log of a negative number) based on variable bounds. + + Args: + stream: an I/O object to write the output to (default = stdout) + + Returns: + None + """ if stream is None: stream = sys.stdout @@ -1702,7 +1710,7 @@ def _check_eval_error_division( node: NumericExpression, warn_list: List[str], config: ConfigDict ): lb, ub = _get_bounds_with_inf(node.args[1]) - if (config.strict_evaluation_error_detection and (lb <= 0 <= ub)) or (lb < 0 < ub): + if (config.warn_for_evaluation_error_at_bounds and (lb <= 0 <= ub)) or (lb < 0 < ub): msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) @@ -1741,7 +1749,7 @@ def _check_eval_error_pow( if integer_exponent and ( (lb1 > 0 or ub1 < 0) - or (not config.strict_evaluation_error_detection and (lb1 >= 0 or ub1 <= 0)) + or (not config.warn_for_evaluation_error_at_bounds and (lb1 >= 0 or ub1 <= 0)) ): # life is good; the exponent is an integer and the base is nonzero return None @@ -1750,7 +1758,7 @@ def _check_eval_error_pow( return None # if the base is positive, there should not be any evaluation errors - if lb1 > 0 or (not config.strict_evaluation_error_detection and lb1 >= 0): + if lb1 > 0 or (not config.warn_for_evaluation_error_at_bounds and lb1 >= 0): return None if lb1 >= 0 and lb2 >= 0: return None @@ -1765,7 +1773,7 @@ def _check_eval_error_log( node: NumericExpression, warn_list: List[str], config: ConfigDict ): lb, ub = _get_bounds_with_inf(node.args[0]) - if (config.strict_evaluation_error_detection and lb <= 0) or lb < 0: + if (config.warn_for_evaluation_error_at_bounds and lb <= 0) or lb < 0: msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" warn_list.append(msg) @@ -1838,6 +1846,14 @@ def __init__(self, config: ConfigDict): self._config = config def exitNode(self, node, data): + """ + callback to be called as the visitor moves from the leaf + nodes back to the root node. + + Args: + node: a pyomo expression node + data: not used in this walker + """ if type(node) in _eval_err_handler: _eval_err_handler[type(node)](node, self._warn_list, self._config) return self._warn_list From 261b29ebfee8a714d62703b19dbfe2554d283c72 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 14 Nov 2023 11:39:18 -0700 Subject: [PATCH 11/13] run black --- idaes/core/util/model_diagnostics.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index ca1d2e2f31..9279bf510b 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -55,7 +55,9 @@ from pyomo.core.base.block import _BlockData from pyomo.core.base.var import _GeneralVarData, _VarData from pyomo.core.base.constraint import _ConstraintData -from pyomo.repn.standard_repn import generate_standard_repn # pylint: disable=no-name-in-module +from pyomo.repn.standard_repn import ( + generate_standard_repn, +) # pylint: disable=no-name-in-module from pyomo.common.collections import ComponentSet from pyomo.common.config import ( ConfigDict, @@ -1355,7 +1357,7 @@ def _collect_potential_eval_errors(self) -> List[str]: def display_potential_evaluation_errors(self, stream=None): """ - Prints constraints that may be prone to evaluation errors + Prints constraints that may be prone to evaluation errors (e.g., log of a negative number) based on variable bounds. Args: @@ -1710,7 +1712,9 @@ def _check_eval_error_division( node: NumericExpression, warn_list: List[str], config: ConfigDict ): lb, ub = _get_bounds_with_inf(node.args[1]) - if (config.warn_for_evaluation_error_at_bounds and (lb <= 0 <= ub)) or (lb < 0 < ub): + if (config.warn_for_evaluation_error_at_bounds and (lb <= 0 <= ub)) or ( + lb < 0 < ub + ): msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" warn_list.append(msg) @@ -1847,7 +1851,7 @@ def __init__(self, config: ConfigDict): def exitNode(self, node, data): """ - callback to be called as the visitor moves from the leaf + callback to be called as the visitor moves from the leaf nodes back to the root node. Args: From 6a1b57e118030390b44758c887e5ab6af07786b3 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 14 Nov 2023 11:53:12 -0700 Subject: [PATCH 12/13] fix pylint issues --- idaes/core/util/model_diagnostics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 9279bf510b..5194ee0064 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -55,9 +55,9 @@ from pyomo.core.base.block import _BlockData from pyomo.core.base.var import _GeneralVarData, _VarData from pyomo.core.base.constraint import _ConstraintData -from pyomo.repn.standard_repn import ( +from pyomo.repn.standard_repn import ( # pylint: disable=no-name-in-module generate_standard_repn, -) # pylint: disable=no-name-in-module +) from pyomo.common.collections import ComponentSet from pyomo.common.config import ( ConfigDict, From a8c0e0e92526c7404b1d07367b319b06c5100907 Mon Sep 17 00:00:00 2001 From: Michael Bynum Date: Tue, 14 Nov 2023 13:14:00 -0700 Subject: [PATCH 13/13] fix tests --- idaes/core/util/tests/test_model_diagnostics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index b86f35dd1b..fe59206ad9 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -2608,7 +2608,7 @@ def test_div(self): w, "c: Potential division by 0 in 1/x; Denominator bounds are (0, inf)" ) - dtb.config.strict_evaluation_error_detection = False + dtb.config.warn_for_evaluation_error_at_bounds = False warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) @@ -2690,7 +2690,7 @@ def test_pow4(self): "c: Potential evaluation error in x**-2; base bounds are (0, inf); exponent bounds are (-2, -2)", ) - dtb.config.strict_evaluation_error_detection = False + dtb.config.warn_for_evaluation_error_at_bounds = False warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) @@ -2718,7 +2718,7 @@ def test_pow5(self): "c: Potential evaluation error in x**-2.5; base bounds are (0, inf); exponent bounds are (-2.5, -2.5)", ) - dtb.config.strict_evaluation_error_detection = False + dtb.config.warn_for_evaluation_error_at_bounds = False warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0) @@ -2750,7 +2750,7 @@ def test_log(self): "c: Potential log of a non-positive number in log(x); Argument bounds are (0, inf)", ) - dtb.config.strict_evaluation_error_detection = False + dtb.config.warn_for_evaluation_error_at_bounds = False warnings = dtb._collect_potential_eval_errors() self.assertEqual(len(warnings), 0)