From f4a17e1debd5e69ae0cc0c1372dc70192910cf8e Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sat, 8 Jun 2024 17:11:31 -0600
Subject: [PATCH 01/69] Update ConstraintData to only store the expression (and
 not lower, body, upper)

---
 pyomo/core/base/constraint.py     | 296 ++++++++++--------------------
 pyomo/core/kernel/constraint.py   |   3 +
 pyomo/core/tests/unit/test_con.py | 248 +++++++++++++++++--------
 3 files changed, 273 insertions(+), 274 deletions(-)

diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index e12860991c2..b3846b44a0a 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -64,6 +64,7 @@
     InequalityExpression,
     RangedExpression,
 }
+_strict_relational_exprs = {True, (False, True), (True, False), (True, True)}
 _rule_returned_none_error = """Constraint '%s': rule returned None.
 
 Constraint rules must return either a valid expression, a 2- or 3-member
@@ -151,7 +152,7 @@ class ConstraintData(ActiveComponentData):
         _active         A boolean that indicates whether this data is active
     """
 
-    __slots__ = ('_body', '_lower', '_upper', '_expr')
+    __slots__ = ('_expr',)
 
     # Set to true when a constraint class stores its expression
     # in linear canonical form
@@ -167,126 +168,139 @@ def __init__(self, expr=None, component=None):
         self._component = weakref_ref(component) if (component is not None) else None
         self._active = True
 
-        self._body = None
-        self._lower = None
-        self._upper = None
         self._expr = None
         if expr is not None:
             self.set_value(expr)
 
     def __call__(self, exception=True):
         """Compute the value of the body of this constraint."""
-        return value(self.body, exception=exception)
+        body = self.normalize_constraint()[1]
+        if body.__class__ not in native_numeric_types:
+            body = value(self.body, exception=exception)
+        return body
+
+    def normalize_constraint(self):
+        expr = self._expr
+        if expr.__class__ is RangedExpression:
+            lb, body, ub = ans = expr.args
+            if (
+                lb.__class__ not in native_types
+                and lb.is_potentially_variable()
+                and not lb.is_fixed()
+            ):
+                raise ValueError(
+                    f"Constraint '{self.name}' is a Ranged Inequality with a "
+                    "variable lower bound.  Cannot normalize the "
+                    "constraint or send it to a solver."
+                )
+            if (
+                ub.__class__ not in native_types
+                and ub.is_potentially_variable()
+                and not ub.is_fixed()
+            ):
+                raise ValueError(
+                    f"Constraint '{self.name}' is a Ranged Inequality with a "
+                    "variable upper bound.  Cannot normalize the "
+                    "constraint or send it to a solver."
+                )
+            return ans
+        elif expr is not None:
+            lhs, rhs = expr.args
+            if rhs.__class__ in native_types or not rhs.is_potentially_variable():
+                return rhs if expr.__class__ is EqualityExpression else None, lhs, rhs
+            if lhs.__class__ in native_types or not lhs.is_potentially_variable():
+                return lhs, rhs, lhs if expr.__class__ is EqualityExpression else None
+            return 0 if expr.__class__ is EqualityExpression else None, lhs - rhs, 0
+        return None, None, None
 
     @property
     def body(self):
         """Access the body of a constraint expression."""
-        if self._body is not None:
-            return self._body
-        # The incoming RangedInequality had a potentially variable
-        # bound.  The "body" is fine, but the bounds may not be
-        # (although the responsibility for those checks lies with the
-        # lower/upper properties)
-        body = self._expr.arg(1)
-        if body.__class__ in native_types and body is not None:
-            return as_numeric(body)
-        return body
-
-    def _get_range_bound(self, range_arg):
-        # Equalities and simple inequalities can always be (directly)
-        # reformulated at construction time to force constant bounds.
-        # The only time we need to defer the determination of bounds is
-        # for ranged inequalities that contain non-constant bounds (so
-        # we *know* that the expr will have 3 args)
-        #
-        # It is possible that there is no expression at all (so catch that)
-        if self._expr is None:
-            return None
-        bound = self._expr.arg(range_arg)
-        if not is_fixed(bound):
-            raise ValueError(
-                "Constraint '%s' is a Ranged Inequality with a "
-                "variable %s bound.  Cannot normalize the "
-                "constraint or send it to a solver."
-                % (self.name, {0: 'lower', 2: 'upper'}[range_arg])
-            )
-        return bound
+        try:
+            ans = self.normalize_constraint()[1]
+        except ValueError:
+            if self._expr.__class__ is RangedExpression:
+                _, ans, _ = self._expr.args
+            else:
+                raise
+        if ans.__class__ in native_numeric_types:
+            return as_numeric(ans)
+        return ans
 
     @property
     def lower(self):
         """Access the lower bound of a constraint expression."""
-        bound = self._lower if self._body is not None else self._get_range_bound(0)
-        # Historically, constraint.lower was guaranteed to return a type
-        # derived from Pyomo NumericValue (or None).  Replicate that
-        # functionality, although clients should in almost all cases
-        # move to using ConstraintData.lb instead of accessing
-        # lower/body/upper to avoid the unnecessary creation (and
-        # inevitable destruction) of the NumericConstant wrappers.
-        if bound is None:
-            return None
-        return as_numeric(bound)
+        ans = self.normalize_constraint()[0]
+        if ans.__class__ in native_numeric_types:
+            # Historically, constraint.lower was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that
+            # functionality, although clients should in almost all cases
+            # move to using ConstraintData.lb instead of accessing
+            # lower/body/upper to avoid the unnecessary creation (and
+            # inevitable destruction) of the NumericConstant wrappers.
+            return as_numeric(ans)
+        return ans
 
     @property
     def upper(self):
         """Access the upper bound of a constraint expression."""
-        bound = self._upper if self._body is not None else self._get_range_bound(2)
-        # Historically, constraint.upper was guaranteed to return a type
-        # derived from Pyomo NumericValue (or None).  Replicate that
-        # functionality, although clients should in almost all cases
-        # move to using ConstraintData.ub instead of accessing
-        # lower/body/upper to avoid the unnecessary creation (and
-        # inevitable destruction) of the NumericConstant wrappers.
-        if bound is None:
-            return None
-        return as_numeric(bound)
+        ans = self.normalize_constraint()[2]
+        if ans.__class__ in native_numeric_types:
+            # Historically, constraint.lower was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that
+            # functionality, although clients should in almost all cases
+            # move to using ConstraintData.lb instead of accessing
+            # lower/body/upper to avoid the unnecessary creation (and
+            # inevitable destruction) of the NumericConstant wrappers.
+            return as_numeric(ans)
+        return ans
 
     @property
     def lb(self):
         """Access the value of the lower bound of a constraint expression."""
-        bound = self._lower if self._body is not None else self._get_range_bound(0)
+        bound = self.normalize_constraint()[0]
+        if bound is None:
+            return None
         if bound.__class__ not in native_numeric_types:
-            if bound is None:
-                return None
             bound = float(value(bound))
+        # Note that "bound != bound" catches float('nan')
         if bound in _nonfinite_values or bound != bound:
-            # Note that "bound != bound" catches float('nan')
             if bound == -_inf:
                 return None
-            else:
-                raise ValueError(
-                    "Constraint '%s' created with an invalid non-finite "
-                    "lower bound (%s)." % (self.name, bound)
-                )
+            raise ValueError(
+                f"Constraint '{self.name}' created with an invalid non-finite "
+                f"lower bound ({bound})."
+            )
         return bound
 
     @property
     def ub(self):
         """Access the value of the upper bound of a constraint expression."""
-        bound = self._upper if self._body is not None else self._get_range_bound(2)
+        bound = self.normalize_constraint()[2]
+        if bound is None:
+            return None
         if bound.__class__ not in native_numeric_types:
-            if bound is None:
-                return None
             bound = float(value(bound))
+        # Note that "bound != bound" catches float('nan')
         if bound in _nonfinite_values or bound != bound:
-            # Note that "bound != bound" catches float('nan')
             if bound == _inf:
                 return None
-            else:
-                raise ValueError(
-                    "Constraint '%s' created with an invalid non-finite "
-                    "upper bound (%s)." % (self.name, bound)
-                )
+            raise ValueError(
+                f"Constraint '{self.name}' created with an invalid non-finite "
+                f"upper bound ({bound})."
+            )
         return bound
 
     @property
     def equality(self):
         """A boolean indicating whether this is an equality constraint."""
-        if self._expr.__class__ is EqualityExpression:
+        expr = self.expr
+        if expr.__class__ is EqualityExpression:
             return True
-        elif self._expr.__class__ is RangedExpression:
+        elif expr.__class__ is RangedExpression:
             # TODO: this is a very restrictive form of structural equality.
-            lb = self._expr.arg(0)
-            if lb is not None and lb is self._expr.arg(2):
+            lb = expr.arg(0)
+            if lb is not None and lb is expr.arg(2):
                 return True
         return False
 
@@ -317,15 +331,22 @@ def expr(self):
 
     def get_value(self):
         """Get the expression on this constraint."""
-        return self._expr
+        return self.expr
 
     def set_value(self, expr):
         """Set the expression on this constraint."""
         # Clear any previously-cached normalized constraint
-        self._lower = self._upper = self._body = self._expr = None
-
+        self._expr = None
         if expr.__class__ in _known_relational_expressions:
+            if getattr(expr, 'strict', False) in _strict_relational_exprs:
+                raise ValueError(
+                    "Constraint '%s' encountered a strict "
+                    "inequality expression ('>' or '< '). All"
+                    " constraints must be formulated using "
+                    "using '<=', '>=', or '=='." % (self.name,)
+                )
             self._expr = expr
+
         elif expr.__class__ is tuple:  # or expr_type is list:
             for arg in expr:
                 if (
@@ -422,120 +443,6 @@ def set_value(self, expr):
                     "\n   (0, model.price[item], 50)" % (self.name, str(expr))
                 )
                 raise ValueError(msg)
-        #
-        # Normalize the incoming expressions, if we can
-        #
-        args = self._expr.args
-        if self._expr.__class__ is InequalityExpression:
-            if self._expr.strict:
-                raise ValueError(
-                    "Constraint '%s' encountered a strict "
-                    "inequality expression ('>' or '< '). All"
-                    " constraints must be formulated using "
-                    "using '<=', '>=', or '=='." % (self.name,)
-                )
-            if (
-                args[1] is None
-                or args[1].__class__ in native_numeric_types
-                or not args[1].is_potentially_variable()
-            ):
-                self._body = args[0]
-                self._upper = args[1]
-            elif (
-                args[0] is None
-                or args[0].__class__ in native_numeric_types
-                or not args[0].is_potentially_variable()
-            ):
-                self._lower = args[0]
-                self._body = args[1]
-            else:
-                self._body = args[0] - args[1]
-                self._upper = 0
-        elif self._expr.__class__ is EqualityExpression:
-            if args[0] is None or args[1] is None:
-                # Error check: ensure equality does not have infinite RHS
-                raise ValueError(
-                    "Equality constraint '%s' defined with "
-                    "non-finite term (%sHS == None)."
-                    % (self.name, 'L' if args[0] is None else 'R')
-                )
-            if (
-                args[0].__class__ in native_numeric_types
-                or not args[0].is_potentially_variable()
-            ):
-                self._lower = self._upper = args[0]
-                self._body = args[1]
-            elif (
-                args[1].__class__ in native_numeric_types
-                or not args[1].is_potentially_variable()
-            ):
-                self._lower = self._upper = args[1]
-                self._body = args[0]
-            else:
-                self._lower = self._upper = 0
-                self._body = args[0] - args[1]
-            # The following logic is caught below when checking for
-            # invalid non-finite bounds:
-            #
-            # if self._lower.__class__ in native_numeric_types and \
-            #    not math.isfinite(self._lower):
-            #     raise ValueError(
-            #         "Equality constraint '%s' defined with "
-            #         "non-finite term." % (self.name))
-        elif self._expr.__class__ is RangedExpression:
-            if any(self._expr.strict):
-                raise ValueError(
-                    "Constraint '%s' encountered a strict "
-                    "inequality expression ('>' or '< '). All"
-                    " constraints must be formulated using "
-                    "using '<=', '>=', or '=='." % (self.name,)
-                )
-            if all(
-                (
-                    arg is None
-                    or arg.__class__ in native_numeric_types
-                    or not arg.is_potentially_variable()
-                )
-                for arg in (args[0], args[2])
-            ):
-                self._lower, self._body, self._upper = args
-        else:
-            # Defensive programming: we currently only support three
-            # relational expression types.  This will only be hit if
-            # someone defines a fourth...
-            raise DeveloperError(
-                "Unrecognized relational expression type: %s"
-                % (self._expr.__class__.__name__,)
-            )
-
-        # We have historically forced the body to be a numeric expression.
-        # TODO: remove this requirement
-        if self._body.__class__ in native_types and self._body is not None:
-            self._body = as_numeric(self._body)
-
-        # We have historically mapped incoming inf to None
-        if self._lower.__class__ in native_numeric_types:
-            bound = self._lower
-            if bound in _nonfinite_values or bound != bound:
-                # Note that "bound != bound" catches float('nan')
-                if bound == -_inf:
-                    self._lower = None
-                else:
-                    raise ValueError(
-                        "Constraint '%s' created with an invalid non-finite "
-                        "lower bound (%s)." % (self.name, self._lower)
-                    )
-        if self._upper.__class__ in native_numeric_types:
-            bound = self._upper
-            if bound in _nonfinite_values or bound != bound:
-                # Note that "bound != bound" catches float('nan')
-                if bound == _inf:
-                    self._upper = None
-                else:
-                    raise ValueError(
-                        "Constraint '%s' created with an invalid non-finite "
-                        "upper bound (%s)." % (self.name, self._upper)
-                    )
 
     def lslack(self):
         """
@@ -911,6 +818,7 @@ class SimpleConstraint(metaclass=RenamedClass):
     {
         'add',
         'set_value',
+        'normalize_constraint',
         'body',
         'lower',
         'upper',
diff --git a/pyomo/core/kernel/constraint.py b/pyomo/core/kernel/constraint.py
index 6aa4abc4bfe..6b8c4c619f5 100644
--- a/pyomo/core/kernel/constraint.py
+++ b/pyomo/core/kernel/constraint.py
@@ -177,6 +177,9 @@ class _MutableBoundsConstraintMixin(object):
     # Define some of the IConstraint abstract methods
     #
 
+    def normalize_constraint(self):
+        return self.lower, self.body, self.upper
+
     @property
     def lower(self):
         """The expression for the lower bound of the constraint"""
diff --git a/pyomo/core/tests/unit/test_con.py b/pyomo/core/tests/unit/test_con.py
index 15f190e281e..7274adae397 100644
--- a/pyomo/core/tests/unit/test_con.py
+++ b/pyomo/core/tests/unit/test_con.py
@@ -84,21 +84,55 @@ def rule(model):
         self.assertEqual(model.c.upper, 0)
 
     def test_tuple_construct_inf_equality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return (model.x, float('inf'))
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return (float('inf'), model.x)
+        model = self.create_model(abstract=True).create_instance()
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.c = Constraint(expr=(model.x, float('inf')))
+        self.assertEqual(model.c.equality, True)
+        self.assertEqual(model.c.lower, float('inf'))
+        self.assertIs(model.c.body, model.x)
+        self.assertEqual(model.c.upper, float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.c.lb
+        self.assertEqual(model.c.ub, None)
+
+        model.d = Constraint(expr=(float('inf'), model.x))
+        self.assertEqual(model.d.equality, True)
+        self.assertEqual(model.d.lower, float('inf'))
+        self.assertIs(model.d.body, model.x)
+        self.assertEqual(model.d.upper, float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.d.lb
+        self.assertEqual(model.d.ub, None)
+
+        model.e = Constraint(expr=(model.x, float('-inf')))
+        self.assertEqual(model.e.equality, True)
+        self.assertEqual(model.e.lower, float('-inf'))
+        self.assertIs(model.e.body, model.x)
+        self.assertEqual(model.e.upper, float('-inf'))
+        self.assertEqual(model.e.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'e' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.e.ub
+
+        model.f = Constraint(expr=(float('-inf'), model.x))
+        self.assertEqual(model.f.equality, True)
+        self.assertEqual(model.f.lower, float('-inf'))
+        self.assertIs(model.f.body, model.x)
+        self.assertEqual(model.f.upper, float('-inf'))
+        self.assertEqual(model.f.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'f' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.f.ub
 
     def test_tuple_construct_1sided_inequality(self):
         model = self.create_model()
@@ -134,9 +168,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, 1)
+        self.assertEqual(model.c.lb, None)
+        self.assertEqual(model.c.ub, 1)
 
         model = self.create_model()
 
@@ -148,7 +184,9 @@ def rule(model):
         self.assertEqual(model.c.equality, False)
         self.assertEqual(model.c.lower, 0)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertEqual(model.c.lb, 0)
+        self.assertEqual(model.c.ub, None)
 
     def test_tuple_construct_unbounded_inequality(self):
         model = self.create_model()
@@ -171,9 +209,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertEqual(model.c.lb, None)
+        self.assertEqual(model.c.ub, None)
 
     def test_tuple_construct_invalid_1sided_inequality(self):
         model = self.create_model(abstract=True)
@@ -229,7 +269,11 @@ def rule(model):
         ):
             instance.c.lower
         self.assertIs(instance.c.body, instance.y)
-        self.assertEqual(instance.c.upper, 1)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            instance.c.upper
         instance.x.fix(3)
         self.assertEqual(value(instance.c.lower), 3)
 
@@ -240,7 +284,11 @@ def rule(model):
 
         model.c = Constraint(rule=rule)
         instance = model.create_instance()
-        self.assertEqual(instance.c.lower, 0)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.lower
         self.assertIs(instance.c.body, instance.y)
         with self.assertRaisesRegex(
             ValueError,
@@ -276,21 +324,23 @@ def rule(model):
         self.assertEqual(model.c.upper, 0)
 
     def test_expr_construct_inf_equality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.x == float('inf')
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model = self.create_model(abstract=True).create_instance()
 
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return float('inf') == model.x
+        model.c = Constraint(expr=model.x == float('inf'))
+        self.assertEqual(model.c.ub, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.c.lb
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.d = Constraint(expr=model.x == float('-inf'))
+        self.assertEqual(model.d.lb, None)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.d.ub
 
     def test_expr_construct_1sided_inequality(self):
         model = self.create_model()
@@ -350,9 +400,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertIs(model.c.lower, None)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -362,9 +414,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, None)
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -374,9 +428,11 @@ def rule(model):
         model.c = Constraint(rule=rule)
 
         self.assertEqual(model.c.equality, False)
-        self.assertEqual(model.c.lower, None)
+        self.assertEqual(model.c.lower, float('-inf'))
         self.assertIs(model.c.body, model.y)
         self.assertEqual(model.c.upper, None)
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
         model = self.create_model()
 
@@ -388,40 +444,40 @@ def rule(model):
         self.assertEqual(model.c.equality, False)
         self.assertEqual(model.c.lower, None)
         self.assertIs(model.c.body, model.y)
-        self.assertEqual(model.c.upper, None)
+        self.assertEqual(model.c.upper, float('inf'))
+        self.assertIs(model.c.ub, None)
+        self.assertIs(model.c.lb, None)
 
     def test_expr_construct_invalid_unbounded_inequality(self):
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.y <= float('-inf')
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return float('inf') <= model.y
-
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
-
-        model = self.create_model(abstract=True)
-
-        def rule(model):
-            return model.y >= float('inf')
+        model = self.create_model(abstract=True).create_instance()
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.c = Constraint(expr=model.y <= float('-inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.c.ub
 
-        model = self.create_model(abstract=True)
+        model.d = Constraint(expr=float('inf') <= model.y)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'd' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.d.lb
 
-        def rule(model):
-            return float('-inf') >= model.y
+        model.e = Constraint(expr=model.y >= float('inf'))
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'e' created with an invalid non-finite lower bound \(inf\).",
+        ):
+            model.e.lb
 
-        model.c = Constraint(rule=rule)
-        self.assertRaises(ValueError, model.create_instance)
+        model.f = Constraint(expr=float('-inf') >= model.y)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'f' created with an invalid non-finite upper bound \(-inf\).",
+        ):
+            model.f.ub
 
     def test_expr_construct_invalid(self):
         m = ConcreteModel()
@@ -484,9 +540,6 @@ def test_nondata_bounds(self):
         model.e2 = Expression()
         model.e3 = Expression()
         model.c.set_value((model.e1, model.e2, model.e3))
-        self.assertIsNone(model.c._lower)
-        self.assertIsNone(model.c._body)
-        self.assertIsNone(model.c._upper)
         self.assertIs(model.c.lower, model.e1)
         self.assertIs(model.c.body, model.e2)
         self.assertIs(model.c.upper, model.e3)
@@ -507,7 +560,7 @@ def test_nondata_bounds(self):
         self.assertIs(model.c.body.expr, model.v[2])
         with self.assertRaisesRegex(
             ValueError,
-            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             model.c.upper
 
@@ -1574,10 +1627,30 @@ def rule1(model):
         self.assertIs(instance.c.body, instance.x)
         with self.assertRaisesRegex(
             ValueError,
-            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+            "Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             instance.c.upper
         #
+        def rule1(model):
+            return (0, model.x, model.z)
+
+        model = AbstractModel()
+        model.x = Var()
+        model.z = Var()
+        model.c = Constraint(rule=rule1)
+        instance = model.create_instance()
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.lower
+        self.assertIs(instance.c.body, instance.x)
+        with self.assertRaisesRegex(
+            ValueError,
+            "Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            instance.c.upper
+
 
     def test_expression_constructor_coverage(self):
         def rule1(model):
@@ -1807,23 +1880,39 @@ def test_potentially_variable_bounds(self):
             r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             m.c.lower
-        self.assertIs(m.c.upper, m.u)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            self.assertIs(m.c.upper, m.u)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             m.c.lb
-        self.assertEqual(m.c.ub, 10)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable lower bound",
+        ):
+            self.assertEqual(m.c.ub, 10)
 
         m.l = 15
         m.u.expr = m.x
-        self.assertIs(m.c.lower, m.l)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            self.assertIs(m.c.lower, m.l)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
         ):
             m.c.upper
-        self.assertEqual(m.c.lb, 15)
+        with self.assertRaisesRegex(
+            ValueError,
+            r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
+        ):
+            self.assertEqual(m.c.lb, 15)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' is a Ranged Inequality with a variable upper bound",
@@ -1890,17 +1979,16 @@ def test_tuple_expression(self):
         ):
             m.c = (m.x, None)
 
+        # You can create it with an infinite value, but then one of the
+        # bounds will fail:
+        m.c = (m.x, float('inf'))
+        self.assertIsNone(m.c.ub)
         with self.assertRaisesRegex(
             ValueError,
             r"Constraint 'c' created with an invalid "
             r"non-finite lower bound \(inf\)",
         ):
-            m.c = (m.x, float('inf'))
-
-        with self.assertRaisesRegex(
-            ValueError, r"Equality constraint 'c' defined with non-finite term"
-        ):
-            m.c = EqualityExpression((m.x, None))
+            m.c.lb
 
 
 if __name__ == "__main__":

From 25027f9fe76e6513a3491976b5c4a2e2f5a14f4a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sat, 8 Jun 2024 17:30:43 -0600
Subject: [PATCH 02/69] Update FBBT to work with raw relational expressions

---
 .../contrib/fbbt/expression_bounds_walker.py  |   6 +-
 pyomo/contrib/fbbt/fbbt.py                    | 105 ++++++++++++++----
 pyomo/contrib/fbbt/interval.py                |  18 +--
 3 files changed, 94 insertions(+), 35 deletions(-)

diff --git a/pyomo/contrib/fbbt/expression_bounds_walker.py b/pyomo/contrib/fbbt/expression_bounds_walker.py
index cb287d54df5..3cb32fcbf29 100644
--- a/pyomo/contrib/fbbt/expression_bounds_walker.py
+++ b/pyomo/contrib/fbbt/expression_bounds_walker.py
@@ -232,15 +232,15 @@ def _handle_unknowable_bounds(visitor, node, arg):
 
 
 def _handle_equality(visitor, node, arg1, arg2):
-    return eq(*arg1, *arg2)
+    return eq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_inequality(visitor, node, arg1, arg2):
-    return ineq(*arg1, *arg2)
+    return ineq(*arg1, *arg2, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_ranged(visitor, node, arg1, arg2, arg3):
-    return ranged(*arg1, *arg2, *arg3)
+    return ranged(*arg1, *arg2, *arg3, feasibility_tol=visitor.feasibility_tol)
 
 
 def _handle_expr_if(visitor, node, arg1, arg2, arg3):
diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 1507c4a3cc5..62cd90d9c87 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -12,6 +12,7 @@
 from collections import defaultdict
 from pyomo.common.collections import ComponentMap, ComponentSet
 from pyomo.contrib.fbbt.expression_bounds_walker import ExpressionBoundsVisitor
+import pyomo.core.expr.relational_expr as relational_expr
 import pyomo.core.expr.numeric_expr as numeric_expr
 from pyomo.core.expr.visitor import (
     ExpressionValueVisitor,
@@ -80,6 +81,24 @@ class FBBTException(PyomoException):
     pass
 
 
+def _prop_bnds_leaf_to_root_equality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.eq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+def _prop_bnds_leaf_to_root_inequality(visitor, node, arg1, arg2):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ineq(
+        *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
+    )
+
+def _prop_bnds_leaf_to_root_ranged(visitor, node, arg1, arg2, arg3):
+    bnds_dict = visitor.bnds_dict
+    bnds_dict[node] = interval.ranged(
+        *bnds_dict[arg1], *bnds_dict[arg2], *bnds_dict[arg3], visitor.feasibility_tol
+    )
+
 def _prop_bnds_leaf_to_root_ProductExpression(visitor, node, arg1, arg2):
     """
 
@@ -367,6 +386,9 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
         numeric_expr.UnaryFunctionExpression: _prop_bnds_leaf_to_root_UnaryFunctionExpression,
         numeric_expr.LinearExpression: _prop_bnds_leaf_to_root_SumExpression,
         numeric_expr.AbsExpression: _prop_bnds_leaf_to_root_abs,
+        relational_expr.EqualityExpression: _prop_bnds_leaf_to_root_equality,
+        relational_expr.InequalityExpression: _prop_bnds_leaf_to_root_inequality,
+        relational_expr.RangedExpression: _prop_bnds_leaf_to_root_ranged,
         ExpressionData: _prop_bnds_leaf_to_root_NamedExpression,
         ScalarExpression: _prop_bnds_leaf_to_root_NamedExpression,
         ObjectiveData: _prop_bnds_leaf_to_root_NamedExpression,
@@ -375,6 +397,51 @@ def _prop_bnds_leaf_to_root_NamedExpression(visitor, node, expr):
 )
 
 
+def _prop_bnds_root_to_leaf_equality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    bnds_dict[arg1] = bnds_dict[arg2] = max(lb1, lb2), min(ub1, ub2)
+
+
+def _prop_bnds_root_to_leaf_inequality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
+def _prop_bnds_root_to_leaf_ranged(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2, arg3 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    lb3, ub3 = bnds_dict[arg3]
+    if lb1 > lb2:
+        bnds_dict[arg2] = lb1, ub2
+        lb2 = lb1
+    if lb2 > lb3:
+        bnds_dict[arg3] = lb2, ub3
+    if ub2 > ub3:
+        bnds_dict[arg2] = lb2, ub3
+        ub2 = ub3
+    if ub1 > ub2:
+        bnds_dict[arg1] = lb1, ub2
+
+
+def _prop_bnds_root_to_leaf_equality(node, bnds_dict, feasibility_tol):
+    assert bnds_dict[node][1]  # This expression is feasible
+    arg1, arg2 = node.args
+    lb1, ub1 = bnds_dict[arg1]
+    lb2, ub2 = bnds_dict[arg2]
+    bnds_dict[arg1] = bnds_dict[arg2] = max(lb1, lb2), min(ub1, ub2)
+
+
 def _prop_bnds_root_to_leaf_ProductExpression(node, bnds_dict, feasibility_tol):
     """
 
@@ -953,6 +1020,15 @@ def _prop_bnds_root_to_leaf_NamedExpression(node, bnds_dict, feasibility_tol):
 _prop_bnds_root_to_leaf_map[ObjectiveData] = _prop_bnds_root_to_leaf_NamedExpression
 _prop_bnds_root_to_leaf_map[ScalarObjective] = _prop_bnds_root_to_leaf_NamedExpression
 
+_prop_bnds_root_to_leaf_map[relational_expr.EqualityExpression] = (
+    _prop_bnds_root_to_leaf_equality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.InequalityExpression] = (
+    _prop_bnds_root_to_leaf_inequality
+)
+_prop_bnds_root_to_leaf_map[relational_expr.RangedExpression] = (
+    _prop_bnds_root_to_leaf_ranged
+)
 
 def _check_and_reset_bounds(var, lb, ub):
     """
@@ -1250,36 +1326,19 @@ def _fbbt_con(con, config):
 
     # a walker to propagate bounds from the variables to the root
     visitorA = _FBBTVisitorLeafToRoot(bnds_dict, feasibility_tol=config.feasibility_tol)
-    visitorA.walk_expression(con.body)
-
-    # Now we need to replace the bounds in bnds_dict for the root
-    # node with the bounds on the constraint (if those bounds are
-    # better).
-    _lb = value(con.lower)
-    _ub = value(con.upper)
-    if _lb is None:
-        _lb = -interval.inf
-    if _ub is None:
-        _ub = interval.inf
+    visitorA.walk_expression(con.expr)
 
-    lb, ub = bnds_dict[con.body]
+    always_feasible, feasible = bnds_dict[con.expr]
 
     # check if the constraint is infeasible
-    if lb > _ub + config.feasibility_tol or ub < _lb - config.feasibility_tol:
+    if not feasible:
         raise InfeasibleConstraintException(
             'Detected an infeasible constraint during FBBT: {0}'.format(str(con))
         )
 
     # check if the constraint is always satisfied
-    if config.deactivate_satisfied_constraints:
-        if lb >= _lb - config.feasibility_tol and ub <= _ub + config.feasibility_tol:
-            con.deactivate()
-
-    if _lb > lb:
-        lb = _lb
-    if _ub < ub:
-        ub = _ub
-    bnds_dict[con.body] = (lb, ub)
+    if config.deactivate_satisfied_constraints and always_feasible:
+        con.deactivate()
 
     # Now, propagate bounds back from the root to the variables
     visitorB = _FBBTVisitorRootToLeaf(
@@ -1287,7 +1346,7 @@ def _fbbt_con(con, config):
         integer_tol=config.integer_tol,
         feasibility_tol=config.feasibility_tol,
     )
-    visitorB.dfs_postorder_stack(con.body)
+    visitorB.dfs_postorder_stack(con.expr)
 
     new_var_bounds = ComponentMap()
     for _node, _bnds in bnds_dict.items():
diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py
index a12d1a4529f..b9117961990 100644
--- a/pyomo/contrib/fbbt/interval.py
+++ b/pyomo/contrib/fbbt/interval.py
@@ -57,7 +57,7 @@ def BoolFlag(val):
     return _true if val else _false
 
 
-def ineq(xl, xu, yl, yu):
+def ineq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an InequalityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -67,9 +67,9 @@ def ineq(xl, xu, yl, yu):
 
     """
     ans = []
-    if yl < xu:
+    if yl < xu - feasibility_tol:
         ans.append(_false)
-    if xl <= yu:
+    if xl <= yu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
@@ -77,7 +77,7 @@ def ineq(xl, xu, yl, yu):
     return tuple(ans)
 
 
-def eq(xl, xu, yl, yu):
+def eq(xl, xu, yl, yu, feasibility_tol):
     """Compute the "bounds" on an EqualityExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -87,9 +87,9 @@ def eq(xl, xu, yl, yu):
 
     """
     ans = []
-    if xl != xu or yl != yu or xl != yl:
+    if abs(xl - xu) > feasibility_tol or abs(yl - yu) > feasibility_tol or abs(xl - yl) > feasibility_tol:
         ans.append(_false)
-    if xl <= yu and yl <= xu:
+    if xl <= yu + feasibility_tol and yl <= xu + feasibility_tol:
         ans.append(_true)
     assert ans
     if len(ans) == 1:
@@ -97,7 +97,7 @@ def eq(xl, xu, yl, yu):
     return tuple(ans)
 
 
-def ranged(xl, xu, yl, yu, zl, zu):
+def ranged(xl, xu, yl, yu, zl, zu, feasibility_tol):
     """Compute the "bounds" on a RangedExpression
 
     Note this is *not* performing interval arithmetic: we are
@@ -106,8 +106,8 @@ def ranged(xl, xu, yl, yu, zl, zu):
     `z` and `z`, `y` can be outside the range `x` and `z`, or both.
 
     """
-    lb = ineq(xl, xu, yl, yu)
-    ub = ineq(yl, yu, zl, zu)
+    lb = ineq(xl, xu, yl, yu, feasibility_tol)
+    ub = ineq(yl, yu, zl, zu, feasibility_tol)
     ans = []
     if not lb[0] or not ub[0]:
         ans.append(_false)

From 1b20539f12f4eefc7f2679de1f83f59bff0ebd4a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sat, 8 Jun 2024 17:34:54 -0600
Subject: [PATCH 03/69] Update transformations to not rely on Constraint
 lower/body/upper

---
 pyomo/core/plugins/transform/add_slack_vars.py | 13 ++++++++-----
 pyomo/gdp/plugins/bilinear.py                  |  7 ++++---
 pyomo/gdp/plugins/cuttingplane.py              |  6 ++++--
 3 files changed, 16 insertions(+), 10 deletions(-)

diff --git a/pyomo/core/plugins/transform/add_slack_vars.py b/pyomo/core/plugins/transform/add_slack_vars.py
index 39903384729..31c1107d692 100644
--- a/pyomo/core/plugins/transform/add_slack_vars.py
+++ b/pyomo/core/plugins/transform/add_slack_vars.py
@@ -150,26 +150,29 @@ def _apply_to_impl(self, instance, **kwds):
             if not cons.active:
                 continue
             cons_name = cons.getname(fully_qualified=True)
-            if cons.lower is not None:
+            lower = cons.lower
+            body = cons.body
+            upper = cons.upper
+            if lower is not None:
                 # we add positive slack variable to body:
                 # declare positive slack
                 varName = "_slack_plus_" + cons_name
                 posSlack = Var(within=NonNegativeReals)
                 xblock.add_component(varName, posSlack)
                 # add positive slack to body expression
-                cons._body += posSlack
+                body += posSlack
                 # penalize slack in objective
                 obj_expr += posSlack
-            if cons.upper is not None:
+            if upper is not None:
                 # we subtract a positive slack variable from the body:
                 # declare slack
                 varName = "_slack_minus_" + cons_name
                 negSlack = Var(within=NonNegativeReals)
                 xblock.add_component(varName, negSlack)
                 # add negative slack to body expression
-                cons._body -= negSlack
+                body -= negSlack
                 # add slack to objective
                 obj_expr += negSlack
-
+            cons.set_value((lower, body, upper))
         # make a new objective that minimizes sum of slack variables
         xblock._slack_objective = Objective(expr=obj_expr)
diff --git a/pyomo/gdp/plugins/bilinear.py b/pyomo/gdp/plugins/bilinear.py
index 67390801348..70b6e83b52f 100644
--- a/pyomo/gdp/plugins/bilinear.py
+++ b/pyomo/gdp/plugins/bilinear.py
@@ -77,9 +77,10 @@ def _transformBlock(self, block, instance):
         for component in block.component_data_objects(
             Constraint, active=True, descend_into=False
         ):
-            expr = self._transformExpression(component.body, instance)
-            instance.bilinear_data_.c_body[id(component)] = component.body
-            component._body = expr
+            lb, body, ub = component.normalize_constraint()
+            expr = self._transformExpression(body, instance)
+            instance.bilinear_data_.c_body[id(component)] = body
+            component.set_value((lb, expr, ub))
 
     def _transformExpression(self, expr, instance):
         if expr.polynomial_degree() > 2:
diff --git a/pyomo/gdp/plugins/cuttingplane.py b/pyomo/gdp/plugins/cuttingplane.py
index 6c77a582987..a757f23c826 100644
--- a/pyomo/gdp/plugins/cuttingplane.py
+++ b/pyomo/gdp/plugins/cuttingplane.py
@@ -400,7 +400,8 @@ def back_off_constraint_with_calculated_cut_violation(
     val = value(transBlock_rHull.infeasibility_objective) - TOL
     if val <= 0:
         logger.info("\tBacking off cut by %s" % val)
-        cut._body += abs(val)
+        lb, body, ub = cut.normalize_constraint()
+        cut.set_value((lb, body + abs(val), ub))
     # else there is nothing to do: restore the objective
     transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective)
     transBlock_rHull.separation_objective.activate()
@@ -424,7 +425,8 @@ def back_off_constraint_by_fixed_tolerance(
                    this callback
     TOL: An absolute tolerance to be added to make cut more conservative.
     """
-    cut._body += TOL
+    lb, body, ub = cut.normalize_constraint()
+    cut.set_value((lb, body + TOL, ub))
 
 
 @TransformationFactory.register(

From 8714e4ca7653da39488a198f7c3271d516b20427 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sat, 8 Jun 2024 17:35:37 -0600
Subject: [PATCH 04/69] Update solver interfaces to not rely on Constraint
 lower/body/upper

---
 pyomo/contrib/solver/persistent.py            | 37 +++----------------
 .../plugins/solvers/persistent_solver.py      |  4 +-
 2 files changed, 8 insertions(+), 33 deletions(-)

diff --git a/pyomo/contrib/solver/persistent.py b/pyomo/contrib/solver/persistent.py
index 71322b7043e..65da81a0c08 100644
--- a/pyomo/contrib/solver/persistent.py
+++ b/pyomo/contrib/solver/persistent.py
@@ -111,8 +111,8 @@ def add_constraints(self, cons: List[ConstraintData]):
                 raise ValueError(
                     'constraint {name} has already been added'.format(name=con.name)
                 )
-            self._active_constraints[con] = (con.lower, con.body, con.upper)
-            tmp = collect_vars_and_named_exprs(con.body)
+            self._active_constraints[con] = con.expr
+            tmp = collect_vars_and_named_exprs(con.expr)
             named_exprs, variables, fixed_vars, external_functions = tmp
             self._check_for_new_vars(variables)
             self._named_expressions[con] = [(e, e.expr) for e in named_exprs]
@@ -417,40 +417,13 @@ def update(self, timer: HierarchicalTimer = None):
         cons_to_remove_and_add = {}
         need_to_set_objective = False
         if config.update_constraints:
-            cons_to_update = []
-            sos_to_update = []
             for c in current_cons_dict.keys():
-                if c not in new_cons_set:
-                    cons_to_update.append(c)
+                if c not in new_cons_set and c.expr is not self._active_constraints[c]:
+                    cons_to_remove_and_add[c] = None
+            sos_to_update = []
             for c in current_sos_dict.keys():
                 if c not in new_sos_set:
                     sos_to_update.append(c)
-            for c in cons_to_update:
-                lower, body, upper = self._active_constraints[c]
-                new_lower, new_body, new_upper = c.lower, c.body, c.upper
-                if new_body is not body:
-                    cons_to_remove_and_add[c] = None
-                    continue
-                if new_lower is not lower:
-                    if (
-                        type(new_lower) is NumericConstant
-                        and type(lower) is NumericConstant
-                        and new_lower.value == lower.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
-                if new_upper is not upper:
-                    if (
-                        type(new_upper) is NumericConstant
-                        and type(upper) is NumericConstant
-                        and new_upper.value == upper.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
             self.remove_sos_constraints(sos_to_update)
             self.add_sos_constraints(sos_to_update)
         timer.stop('cons')
diff --git a/pyomo/solvers/plugins/solvers/persistent_solver.py b/pyomo/solvers/plugins/solvers/persistent_solver.py
index 3c2a9e52eab..ef96bfa339f 100644
--- a/pyomo/solvers/plugins/solvers/persistent_solver.py
+++ b/pyomo/solvers/plugins/solvers/persistent_solver.py
@@ -262,7 +262,9 @@ def _add_and_collect_column_data(self, var, obj_coef, constraints, coefficients)
         coeff_list = list()
         constr_list = list()
         for val, c in zip(coefficients, constraints):
-            c._body += val * var
+            lb, body, ub = c.normalize_constraint()
+            body += val * var
+            c.set_value((lb, body, ub))
             self._vars_referenced_by_con[c].add(var)
 
             cval = _convert_to_const(val)

From 7aee9e04c907cae3c936954110be1e664a736160 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sun, 16 Jun 2024 14:36:28 -0600
Subject: [PATCH 05/69] NFC: apply black

---
 pyomo/contrib/fbbt/fbbt.py        | 4 ++++
 pyomo/contrib/fbbt/interval.py    | 6 +++++-
 pyomo/core/tests/unit/test_con.py | 2 +-
 3 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 62cd90d9c87..593d875ca6f 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -87,18 +87,21 @@ def _prop_bnds_leaf_to_root_equality(visitor, node, arg1, arg2):
         *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
     )
 
+
 def _prop_bnds_leaf_to_root_inequality(visitor, node, arg1, arg2):
     bnds_dict = visitor.bnds_dict
     bnds_dict[node] = interval.ineq(
         *bnds_dict[arg1], *bnds_dict[arg2], visitor.feasibility_tol
     )
 
+
 def _prop_bnds_leaf_to_root_ranged(visitor, node, arg1, arg2, arg3):
     bnds_dict = visitor.bnds_dict
     bnds_dict[node] = interval.ranged(
         *bnds_dict[arg1], *bnds_dict[arg2], *bnds_dict[arg3], visitor.feasibility_tol
     )
 
+
 def _prop_bnds_leaf_to_root_ProductExpression(visitor, node, arg1, arg2):
     """
 
@@ -1030,6 +1033,7 @@ def _prop_bnds_root_to_leaf_NamedExpression(node, bnds_dict, feasibility_tol):
     _prop_bnds_root_to_leaf_ranged
 )
 
+
 def _check_and_reset_bounds(var, lb, ub):
     """
     This function ensures that lb is not less than var.lb and that ub is not greater than var.ub.
diff --git a/pyomo/contrib/fbbt/interval.py b/pyomo/contrib/fbbt/interval.py
index b9117961990..4b93d6e3f31 100644
--- a/pyomo/contrib/fbbt/interval.py
+++ b/pyomo/contrib/fbbt/interval.py
@@ -87,7 +87,11 @@ def eq(xl, xu, yl, yu, feasibility_tol):
 
     """
     ans = []
-    if abs(xl - xu) > feasibility_tol or abs(yl - yu) > feasibility_tol or abs(xl - yl) > feasibility_tol:
+    if (
+        abs(xl - xu) > feasibility_tol
+        or abs(yl - yu) > feasibility_tol
+        or abs(xl - yl) > feasibility_tol
+    ):
         ans.append(_false)
     if xl <= yu + feasibility_tol and yl <= xu + feasibility_tol:
         ans.append(_true)
diff --git a/pyomo/core/tests/unit/test_con.py b/pyomo/core/tests/unit/test_con.py
index 7274adae397..07c7eb3af8e 100644
--- a/pyomo/core/tests/unit/test_con.py
+++ b/pyomo/core/tests/unit/test_con.py
@@ -1630,6 +1630,7 @@ def rule1(model):
             "Constraint 'c' is a Ranged Inequality with a variable lower bound",
         ):
             instance.c.upper
+
         #
         def rule1(model):
             return (0, model.x, model.z)
@@ -1651,7 +1652,6 @@ def rule1(model):
         ):
             instance.c.upper
 
-
     def test_expression_constructor_coverage(self):
         def rule1(model):
             expr = model.x

From 3d77bc9b400600669d8445f2c1d9ee36779e6fb1 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Sun, 16 Jun 2024 23:10:21 -0600
Subject: [PATCH 06/69] Switch from caching .body to caching .expr

---
 pyomo/contrib/appsi/base.py | 35 ++++-------------------------------
 1 file changed, 4 insertions(+), 31 deletions(-)

diff --git a/pyomo/contrib/appsi/base.py b/pyomo/contrib/appsi/base.py
index 6d2b5ccfcd4..9c7da1eb60b 100644
--- a/pyomo/contrib/appsi/base.py
+++ b/pyomo/contrib/appsi/base.py
@@ -1007,7 +1007,7 @@ def add_constraints(self, cons: List[ConstraintData]):
                 raise ValueError(
                     'constraint {name} has already been added'.format(name=con.name)
                 )
-            self._active_constraints[con] = (con.lower, con.body, con.upper)
+            self._active_constraints[con] = con.expr
             if self.use_extensions and cmodel_available:
                 tmp = cmodel.prep_for_repn(con.body, self._expr_types)
             else:
@@ -1363,40 +1363,13 @@ def update(self, timer: HierarchicalTimer = None):
         cons_to_remove_and_add = dict()
         need_to_set_objective = False
         if config.update_constraints:
-            cons_to_update = list()
-            sos_to_update = list()
             for c in current_cons_dict.keys():
-                if c not in new_cons_set:
-                    cons_to_update.append(c)
+                if c not in new_cons_set and c.expr is not self._active_constraints[c]:
+                    cons_to_remove_and_add[c] = None
+            sos_to_update = []
             for c in current_sos_dict.keys():
                 if c not in new_sos_set:
                     sos_to_update.append(c)
-            for c in cons_to_update:
-                lower, body, upper = self._active_constraints[c]
-                new_lower, new_body, new_upper = c.lower, c.body, c.upper
-                if new_body is not body:
-                    cons_to_remove_and_add[c] = None
-                    continue
-                if new_lower is not lower:
-                    if (
-                        type(new_lower) is NumericConstant
-                        and type(lower) is NumericConstant
-                        and new_lower.value == lower.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
-                if new_upper is not upper:
-                    if (
-                        type(new_upper) is NumericConstant
-                        and type(upper) is NumericConstant
-                        and new_upper.value == upper.value
-                    ):
-                        pass
-                    else:
-                        cons_to_remove_and_add[c] = None
-                        continue
             self.remove_sos_constraints(sos_to_update)
             self.add_sos_constraints(sos_to_update)
         timer.stop('cons')

From d8de8100bb08124173e3ab125555a3c8785e7ad2 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 17 Jun 2024 05:15:44 -0600
Subject: [PATCH 07/69] Update APPSI cmodel to call normalize_constraint()

---
 pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp | 2 +-
 pyomo/contrib/appsi/cmodel/src/lp_writer.cpp  | 2 +-
 pyomo/contrib/appsi/cmodel/src/nl_writer.cpp  | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
index bd8d7dbf854..708cfd9e073 100644
--- a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
@@ -205,7 +205,7 @@ void process_fbbt_constraints(FBBTModel *model, PyomoExprTypes &expr_types,
   py::handle con_body;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("normalize_constraint")();
     con_lb = lower_body_upper[0];
     con_body = lower_body_upper[1];
     con_ub = lower_body_upper[2];
diff --git a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
index 68baf2b8ae8..996bb34f564 100644
--- a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
@@ -289,7 +289,7 @@ void process_lp_constraints(py::list cons, py::object writer) {
   py::object nonlinear_expr;
   PyomoExprTypes expr_types = PyomoExprTypes();
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("normalize_constraint")();
     cname = getSymbol(c, labeler);
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = true);
diff --git a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
index 8de6cc74ab4..477bdd87aee 100644
--- a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
@@ -527,7 +527,7 @@ void process_nl_constraints(NLWriter *nl_writer, PyomoExprTypes &expr_types,
   py::handle repn_nonlinear_expr;
 
   for (py::handle c : cons) {
-    lower_body_upper = active_constraints[c];
+    lower_body_upper = c.attr("normalize_constraint")();
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = false);
     _const = appsi_expr_from_pyomo_expr(repn.attr("constant"), var_map,

From 496a6638f93041055dff951b6a7ec1d94dd44d6a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 17 Jun 2024 08:06:07 -0600
Subject: [PATCH 08/69] Remove repeated function

---
 pyomo/contrib/fbbt/fbbt.py | 8 --------
 1 file changed, 8 deletions(-)

diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 593d875ca6f..39df91675a8 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -437,14 +437,6 @@ def _prop_bnds_root_to_leaf_ranged(node, bnds_dict, feasibility_tol):
         bnds_dict[arg1] = lb1, ub2
 
 
-def _prop_bnds_root_to_leaf_equality(node, bnds_dict, feasibility_tol):
-    assert bnds_dict[node][1]  # This expression is feasible
-    arg1, arg2 = node.args
-    lb1, ub1 = bnds_dict[arg1]
-    lb2, ub2 = bnds_dict[arg2]
-    bnds_dict[arg1] = bnds_dict[arg2] = max(lb1, lb2), min(ub1, ub2)
-
-
 def _prop_bnds_root_to_leaf_ProductExpression(node, bnds_dict, feasibility_tol):
     """
 

From 0b0833f1149467ecf2c0c9c2539b9ca08e5fa602 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 17 Jun 2024 08:06:22 -0600
Subject: [PATCH 09/69] Switch logic for mapping lower/body/upper to match
 previous behavior

---
 pyomo/core/base/constraint.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index b3846b44a0a..3ee5a82ef58 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -223,7 +223,13 @@ def body(self):
                 _, ans, _ = self._expr.args
             else:
                 raise
-        if ans.__class__ in native_numeric_types:
+        if ans.__class__ in native_types and ans is not None:
+            # Historically, constraint.lower was guaranteed to return a type
+            # derived from Pyomo NumericValue (or None).  Replicate that.
+            #
+            # [JDS 6/2024: it would be nice to remove this behavior,
+            # although possibly unnecessary, as people should use
+            # normalize_constraint() instead]
             return as_numeric(ans)
         return ans
 
@@ -231,7 +237,7 @@ def body(self):
     def lower(self):
         """Access the lower bound of a constraint expression."""
         ans = self.normalize_constraint()[0]
-        if ans.__class__ in native_numeric_types:
+        if ans.__class__ in native_types and ans is not None:
             # Historically, constraint.lower was guaranteed to return a type
             # derived from Pyomo NumericValue (or None).  Replicate that
             # functionality, although clients should in almost all cases
@@ -245,8 +251,8 @@ def lower(self):
     def upper(self):
         """Access the upper bound of a constraint expression."""
         ans = self.normalize_constraint()[2]
-        if ans.__class__ in native_numeric_types:
-            # Historically, constraint.lower was guaranteed to return a type
+        if ans.__class__ in native_types and ans is not None:
+            # Historically, constraint.upper was guaranteed to return a type
             # derived from Pyomo NumericValue (or None).  Replicate that
             # functionality, although clients should in almost all cases
             # move to using ConstraintData.lb instead of accessing

From cdfb2b17d5d78147a0a335b996934405cbc5c8c0 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 17 Jun 2024 08:09:08 -0600
Subject: [PATCH 10/69] Update baseline due to changes in Constraint

(something caused gdpopt to set variables to int instead of floats)
---
 doc/OnlineDocs/contributed_packages/gdpopt.rst | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/OnlineDocs/contributed_packages/gdpopt.rst b/doc/OnlineDocs/contributed_packages/gdpopt.rst
index d550b0ced76..670d7633f6d 100644
--- a/doc/OnlineDocs/contributed_packages/gdpopt.rst
+++ b/doc/OnlineDocs/contributed_packages/gdpopt.rst
@@ -93,10 +93,10 @@ An example that includes the modeling approach may be found below.
     Variables:
       x : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :  -1.2 :   0.0 :     2 : False : False :  Reals
+          None :  -1.2 :     0 :     2 : False : False :  Reals
       y : Size=1, Index=None
           Key  : Lower : Value : Upper : Fixed : Stale : Domain
-          None :   -10 :   1.0 :    10 : False : False :  Reals
+          None :   -10 :     1 :    10 : False : False :  Reals
   <BLANKLINE>
     Objectives:
       objective : Size=1, Index=None, Active=True
@@ -106,7 +106,7 @@ An example that includes the modeling approach may be found below.
     Constraints:
       c : Size=1
           Key  : Lower : Body : Upper
-          None :   1.0 :  1.0 :   1.0
+          None :   1.0 :    1 :   1.0
 
 .. note:: 
 

From 64e01bbfcf84d3664925768943ceae166ae61636 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 17 Jun 2024 12:17:11 -0600
Subject: [PATCH 11/69] NFC: improve variable naming

---
 pyomo/contrib/fbbt/fbbt.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 39df91675a8..60ac0603388 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -1324,10 +1324,10 @@ def _fbbt_con(con, config):
     visitorA = _FBBTVisitorLeafToRoot(bnds_dict, feasibility_tol=config.feasibility_tol)
     visitorA.walk_expression(con.expr)
 
-    always_feasible, feasible = bnds_dict[con.expr]
+    always_feasible, possibly_feasible = bnds_dict[con.expr]
 
     # check if the constraint is infeasible
-    if not feasible:
+    if not possibly_feasible:
         raise InfeasibleConstraintException(
             'Detected an infeasible constraint during FBBT: {0}'.format(str(con))
         )

From 55e022461e14d659a0700c74625119cc019b4a9a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:05:13 -0600
Subject: [PATCH 12/69] Add FreeBSD to the set o known OSes

---
 pyomo/common/fileutils.py | 9 ++++++++-
 1 file changed, 8 insertions(+), 1 deletion(-)

diff --git a/pyomo/common/fileutils.py b/pyomo/common/fileutils.py
index 7b6520327a0..8c3c6dfecaa 100644
--- a/pyomo/common/fileutils.py
+++ b/pyomo/common/fileutils.py
@@ -286,10 +286,17 @@ def find_dir(
     )
 
 
-_exeExt = {'linux': None, 'windows': '.exe', 'cygwin': '.exe', 'darwin': None}
+_exeExt = {
+    'linux': None,
+    'freebsd': None,
+    'windows': '.exe',
+    'cygwin': '.exe',
+    'darwin': None,
+}
 
 _libExt = {
     'linux': ('.so', '.so.*'),
+    'freebsd': ('.so', '.so.*'),
     'windows': ('.dll', '.pyd'),
     'cygwin': ('.dll', '.so', '.so.*'),
     'darwin': ('.dylib', '.so', '.so.*'),

From 615c646b4510654a7d9d83ca762c15b5867a9fbd Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:05:39 -0600
Subject: [PATCH 13/69] Fix typo in test (exercised by unknown playform)

---
 pyomo/common/tests/test_download.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/common/tests/test_download.py b/pyomo/common/tests/test_download.py
index 8fee0ba7e31..4fde029f1b1 100644
--- a/pyomo/common/tests/test_download.py
+++ b/pyomo/common/tests/test_download.py
@@ -206,7 +206,7 @@ def test_get_os_version(self):
             self.assertEqual(_os, 'win')
             self.assertEqual(_norm, _os + ''.join(_ver.split('.')[:2]))
         else:
-            self.assertEqual(ans, '')
+            self.assertEqual(_os, '')
 
         self.assertEqual((_os, _ver), FileDownloader._os_version)
         # Exercise the fetch from CACHE

From e2941a433fff6043c8a2a6fa990939ee0d6df00b Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:06:34 -0600
Subject: [PATCH 14/69] Add `has_linear_solver()` to the IPOPT solver
 interfaces

---
 pyomo/contrib/appsi/solvers/ipopt.py   | 20 ++++++++++++++++++++
 pyomo/contrib/solver/ipopt.py          | 15 +++++++++++++++
 pyomo/solvers/plugins/solvers/IPOPT.py | 11 +++++++++++
 3 files changed, 46 insertions(+)

diff --git a/pyomo/contrib/appsi/solvers/ipopt.py b/pyomo/contrib/appsi/solvers/ipopt.py
index 76cd204e36d..97e8122fe78 100644
--- a/pyomo/contrib/appsi/solvers/ipopt.py
+++ b/pyomo/contrib/appsi/solvers/ipopt.py
@@ -567,3 +567,23 @@ def get_reduced_costs(
             return ComponentMap((k, v) for k, v in self._reduced_costs.items())
         else:
             return ComponentMap((v, self._reduced_costs[v]) for v in vars_to_load)
+
+    def has_linear_solver(self, linear_solver):
+        import pyomo.core as AML
+        from pyomo.common.tee import capture_output
+        m = AML.ConcreteModel()
+        m.x = AML.Var()
+        m.o = AML.Objective(expr=(m.x-2)**2)
+        with capture_output() as OUT:
+            solver = self.__class__()
+            solver.config.stream_solver = True
+            solver.config.load_solution = False
+            solver.ipopt_options['linear_solver'] = linear_solver
+            try:
+                solver.solve(m)
+            except FileNotFoundError:
+                # The APPSI interface always tries to open the SOL file,
+                # and will generate a FileNotFoundError if ipopt didn't
+                # generate one
+                return False
+        return 'running with linear solver' in OUT.getvalue()
diff --git a/pyomo/contrib/solver/ipopt.py b/pyomo/contrib/solver/ipopt.py
index c88696f531b..c467d283d9b 100644
--- a/pyomo/contrib/solver/ipopt.py
+++ b/pyomo/contrib/solver/ipopt.py
@@ -238,6 +238,21 @@ def version(self, config=None):
                 self._version_cache = (pth, version)
         return self._version_cache[1]
 
+    def has_linear_solver(self, linear_solver):
+        import pyomo.core as AML
+
+        m = AML.ConcreteModel()
+        m.x = AML.Var()
+        m.o = AML.Objective(expr=(m.x - 2) ** 2)
+        results = self.solve(
+            m,
+            tee=False,
+            raise_exception_on_nonoptimal_result=False,
+            load_solutions=False,
+            solver_options={'linear_solver': linear_solver},
+        )
+        return 'running with linear solver' in results.solver_log
+
     def _write_options_file(self, filename: str, options: Mapping):
         # First we need to determine if we even need to create a file.
         # If options is empty, then we return False
diff --git a/pyomo/solvers/plugins/solvers/IPOPT.py b/pyomo/solvers/plugins/solvers/IPOPT.py
index 21045cb7b4f..be0f143ea46 100644
--- a/pyomo/solvers/plugins/solvers/IPOPT.py
+++ b/pyomo/solvers/plugins/solvers/IPOPT.py
@@ -14,6 +14,7 @@
 
 from pyomo.common import Executable
 from pyomo.common.collections import Bunch
+from pyomo.common.tee import capture_output
 from pyomo.common.tempfiles import TempfileManager
 
 from pyomo.opt.base import ProblemFormat, ResultsFormat
@@ -207,3 +208,13 @@ def process_output(self, rc):
                             res.solver.message = line.split(':')[2].strip()
                             assert "degrees of freedom" in res.solver.message
             return res
+
+    def has_linear_solver(self, linear_solver):
+        import pyomo.core as AML
+
+        m = AML.ConcreteModel()
+        m.x = AML.Var()
+        m.o = AML.Objective(expr=(m.x - 2) ** 2)
+        with capture_output() as OUT:
+            self.solve(m, tee=True, options={'linear_solver': linear_solver})
+        return 'running with linear solver' in OUT.getvalue()

From f6fbfd47cad42dd6597d74621d1bcfa659208162 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:19:51 -0600
Subject: [PATCH 15/69] Fix log message string formatting

---
 pyomo/contrib/doe/doe.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index a120add4200..e92aefee651 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -995,7 +995,7 @@ def run_grid_search(
                 )
                 count += 1
                 failed_count += 1
-                self.logger.warning("failed count:", failed_count)
+                self.logger.warning("failed count: %s", failed_count)
                 result_combine[tuple(design_set_iter)] = None
 
         # For user's access

From 800fe28475d3abae6a8983525c88b29bf6225747 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:23:27 -0600
Subject: [PATCH 16/69] DoE: only specify linear solver if it is available.

---
 pyomo/contrib/doe/doe.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/doe/doe.py b/pyomo/contrib/doe/doe.py
index e92aefee651..90818ddf622 100644
--- a/pyomo/contrib/doe/doe.py
+++ b/pyomo/contrib/doe/doe.py
@@ -97,7 +97,7 @@ def __init__(
             A Python ``function`` that returns a Concrete Pyomo model, similar to the interface for ``parmest``
         solver:
             A ``solver`` object that User specified, default=None.
-            If not specified, default solver is IPOPT MA57.
+            If not specified, default solver is IPOPT (with MA57, if available).
         prior_FIM:
             A 2D numpy array containing Fisher information matrix (FIM) for prior experiments.
             The default None means there is no prior information.
@@ -1387,7 +1387,10 @@ def _fix_design(self, m, design_val, fix_opt=True, optimize_option=None):
     def _get_default_ipopt_solver(self):
         """Default solver"""
         solver = SolverFactory("ipopt")
-        solver.options["linear_solver"] = "ma57"
+        for linear_solver in ('ma57', 'ma27', 'ma97'):
+            if solver.has_linear_solver(linear_solver):
+                solver.options["linear_solver"] = linear_solver
+                break
         solver.options["halt_on_ampl_error"] = "yes"
         solver.options["max_iter"] = 3000
         return solver

From 308d2a95bb3322ad7d52ff263c0c1b9625ec0b8f Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 07:26:35 -0600
Subject: [PATCH 17/69] ensure testing global sets are unregistered

---
 pyomo/core/tests/unit/test_set.py | 52 ++++++++++++++++++-------------
 1 file changed, 30 insertions(+), 22 deletions(-)

diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index f62589a6873..8b4c567b9ca 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -3518,21 +3518,25 @@ def test_iteration(self):
     def test_declare(self):
         NS = {}
         DeclareGlobalSet(RangeSet(name='TrinarySet', ranges=(NR(0, 2, 1),)), NS)
-        self.assertEqual(list(NS['TrinarySet']), [0, 1, 2])
-        a = pickle.loads(pickle.dumps(NS['TrinarySet']))
-        self.assertIs(a, NS['TrinarySet'])
-        with self.assertRaisesRegex(NameError, "name 'TrinarySet' is not defined"):
-            TrinarySet
-        del SetModule.GlobalSets['TrinarySet']
-        del NS['TrinarySet']
+        try:
+            self.assertEqual(list(NS['TrinarySet']), [0, 1, 2])
+            a = pickle.loads(pickle.dumps(NS['TrinarySet']))
+            self.assertIs(a, NS['TrinarySet'])
+            with self.assertRaisesRegex(NameError, "name 'TrinarySet' is not defined"):
+                TrinarySet
+        finally:
+            del SetModule.GlobalSets['TrinarySet']
+            del NS['TrinarySet']
 
         # Now test the automatic identification of the globals() scope
         DeclareGlobalSet(RangeSet(name='TrinarySet', ranges=(NR(0, 2, 1),)))
-        self.assertEqual(list(TrinarySet), [0, 1, 2])
-        a = pickle.loads(pickle.dumps(TrinarySet))
-        self.assertIs(a, TrinarySet)
-        del SetModule.GlobalSets['TrinarySet']
-        del globals()['TrinarySet']
+        try:
+            self.assertEqual(list(TrinarySet), [0, 1, 2])
+            a = pickle.loads(pickle.dumps(TrinarySet))
+            self.assertIs(a, TrinarySet)
+        finally:
+            del SetModule.GlobalSets['TrinarySet']
+            del globals()['TrinarySet']
         with self.assertRaisesRegex(NameError, "name 'TrinarySet' is not defined"):
             TrinarySet
 
@@ -3551,18 +3555,22 @@ def test_exceptions(self):
 
         NS = {}
         ts = DeclareGlobalSet(RangeSet(name='TrinarySet', ranges=(NR(0, 2, 1),)), NS)
-        self.assertIs(NS['TrinarySet'], ts)
+        try:
+            self.assertIs(NS['TrinarySet'], ts)
 
-        # Repeat declaration is OK
-        DeclareGlobalSet(ts, NS)
-        self.assertIs(NS['TrinarySet'], ts)
+            # Repeat declaration is OK
+            DeclareGlobalSet(ts, NS)
+            self.assertIs(NS['TrinarySet'], ts)
 
-        # but conflicting one raises exception
-        NS['foo'] = None
-        with self.assertRaisesRegex(
-            RuntimeError, "Refusing to overwrite global object, foo"
-        ):
-            DeclareGlobalSet(RangeSet(name='foo', ranges=(NR(0, 2, 1),)), NS)
+            # but conflicting one raises exception
+            NS['foo'] = None
+            with self.assertRaisesRegex(
+                RuntimeError, "Refusing to overwrite global object, foo"
+            ):
+                DeclareGlobalSet(RangeSet(name='foo', ranges=(NR(0, 2, 1),)), NS)
+        finally:
+            del SetModule.GlobalSets['TrinarySet']
+            del NS['TrinarySet']
 
     def test_RealSet_IntegerSet(self):
         output = StringIO()

From 3f843402f9a671d443871cc8b3f497e77ffff5d8 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 08:18:47 -0600
Subject: [PATCH 18/69] Fix state bleedover between tests

---
 pyomo/solvers/tests/mip/test_factory.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/pyomo/solvers/tests/mip/test_factory.py b/pyomo/solvers/tests/mip/test_factory.py
index 31d47486aa4..f69fd198009 100644
--- a/pyomo/solvers/tests/mip/test_factory.py
+++ b/pyomo/solvers/tests/mip/test_factory.py
@@ -53,8 +53,8 @@ def setUpClass(cls):
 
     def tearDown(self):
         ReaderFactory.unregister('rtest3')
-        ReaderFactory.unregister('stest3')
-        ReaderFactory.unregister('wtest3')
+        SolverFactory.unregister('stest3')
+        WriterFactory.unregister('wtest3')
 
     def test_solver_factory(self):
         """
@@ -119,6 +119,9 @@ def test_writer_instance(self):
         ans = WriterFactory("none")
         self.assertEqual(ans, None)
         ans = WriterFactory("wtest3")
+        self.assertEqual(ans, None)
+        WriterFactory.register('wtest3')(MockWriter)
+        ans = WriterFactory("wtest3")
         self.assertNotEqual(ans, None)
 
     def test_writer_registration(self):

From c61d2d252300035c51a0eb05908337058379d787 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Thu, 20 Jun 2024 08:40:06 -0600
Subject: [PATCH 19/69] Guard tests that require k_aug

---
 pyomo/contrib/doe/tests/test_example.py         | 2 ++
 pyomo/contrib/doe/tests/test_reactor_example.py | 2 ++
 2 files changed, 4 insertions(+)

diff --git a/pyomo/contrib/doe/tests/test_example.py b/pyomo/contrib/doe/tests/test_example.py
index e4ffbe89142..47ce39d596a 100644
--- a/pyomo/contrib/doe/tests/test_example.py
+++ b/pyomo/contrib/doe/tests/test_example.py
@@ -39,6 +39,7 @@
 from pyomo.opt import SolverFactory
 
 ipopt_available = SolverFactory("ipopt").available()
+k_aug_available = SolverFactory("k_aug").available(exception_flag=False)
 
 
 class TestReactorExamples(unittest.TestCase):
@@ -57,6 +58,7 @@ def test_reactor_optimize_doe(self):
 
         reactor_optimize_doe.main()
 
+    @unittest.skipIf(not k_aug_available, "The 'k_aug' command is not available")
     @unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
     @unittest.skipIf(not pandas_available, "pandas is not available")
     @unittest.skipIf(not numpy_available, "Numpy is not available")
diff --git a/pyomo/contrib/doe/tests/test_reactor_example.py b/pyomo/contrib/doe/tests/test_reactor_example.py
index 19fb4e61820..f88ae48db1a 100644
--- a/pyomo/contrib/doe/tests/test_reactor_example.py
+++ b/pyomo/contrib/doe/tests/test_reactor_example.py
@@ -35,6 +35,7 @@
 from pyomo.opt import SolverFactory
 
 ipopt_available = SolverFactory("ipopt").available()
+k_aug_available = SolverFactory("k_aug").available(exception_flag=False)
 
 
 class Test_Reaction_Kinetics_Example(unittest.TestCase):
@@ -133,6 +134,7 @@ def test_kinetics_example_sequential_finite_then_optimize(self):
         # self.assertAlmostEqual(value(optimize_result.model.T[0.5]), 300, places=2)
         self.assertAlmostEqual(np.log10(optimize_result.trace), 3.340, places=2)
 
+    @unittest.skipIf(not k_aug_available, "The 'k_aug' solver is not available")
     @unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available")
     @unittest.skipIf(not numpy_available, "Numpy is not available")
     @unittest.skipIf(not pandas_available, "Pandas is not available")

From 2455454a98ede2340c6f59545bc047c54c490170 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 21 Jun 2024 10:45:55 -0600
Subject: [PATCH 20/69] Update test guards to correctly skip when pynumero_ASL
 missing

---
 pyomo/contrib/parmest/tests/test_examples.py |  3 +-
 pyomo/contrib/parmest/tests/test_parmest.py  | 74 ++++++++------------
 2 files changed, 31 insertions(+), 46 deletions(-)

diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index dca05026e80..450863a08a4 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -12,10 +12,11 @@
 import pyomo.common.unittest as unittest
 import pyomo.contrib.parmest.parmest as parmest
 from pyomo.contrib.parmest.graphics import matplotlib_available, seaborn_available
+from pyomo.contrib.pynumero.asl import AmplInterface
 from pyomo.opt import SolverFactory
 
 ipopt_available = SolverFactory("ipopt").available()
-
+pynumero_ASL_available = AmplInterface.available()
 
 @unittest.skipIf(
     not parmest.parmest_available,
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index 65e2e4a3b06..1ff42a38d9e 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -9,44 +9,35 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-from pyomo.common.dependencies import (
-    numpy as np,
-    numpy_available,
-    pandas as pd,
-    pandas_available,
-    scipy,
-    scipy_available,
-    matplotlib,
-    matplotlib_available,
-)
-
 import platform
-
-is_osx = platform.mac_ver()[0] != ""
-
-import pyomo.common.unittest as unittest
 import sys
 import os
 import subprocess
 from itertools import product
 
+import pyomo.common.unittest as unittest
 import pyomo.contrib.parmest.parmest as parmest
 import pyomo.contrib.parmest.graphics as graphics
 import pyomo.contrib.parmest as parmestbase
-from pyomo.contrib.parmest.experiment import Experiment
 import pyomo.environ as pyo
 import pyomo.dae as dae
 
+from pyomo.common.dependencies import (
+    numpy as np,
+    pandas as pd,
+    scipy,
+    matplotlib,
+)
+from pyomo.common.fileutils import this_file_dir
+from pyomo.contrib.parmest.experiment import Experiment
+from pyomo.contrib.pynumero.asl import AmplInterface
 from pyomo.opt import SolverFactory
 
+is_osx = platform.mac_ver()[0] != ""
 ipopt_available = SolverFactory("ipopt").available()
-
-from pyomo.common.fileutils import find_library
-
-pynumero_ASL_available = False if find_library("pynumero_ASL") is None else True
-
-testdir = os.path.dirname(os.path.abspath(__file__))
-
+k_aug_available = SolverFactory("k_aug").available(exception_flag=False)
+pynumero_ASL_available = AmplInterface.available()
+testdir = this_file_dir()
 
 @unittest.skipIf(
     not parmest.parmest_available,
@@ -208,17 +199,13 @@ def test_parallel_parmest(self):
             retcode = subprocess.call(rlist)
         self.assertEqual(retcode, 0)
 
-    @unittest.skip("Most folks don't have k_aug installed")
+    @unittest.skipUnless(k_aug_available, "k_aug solver not found")
     def test_theta_k_aug_for_Hessian(self):
         # this will fail if k_aug is not installed
         objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug")
         self.assertAlmostEqual(objval, 4.4675, places=2)
 
-    @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
-    @unittest.skipIf(
-        not parmest.inverse_reduced_hessian_available,
-        "Cannot test covariance matrix: required ASL dependency is missing",
-    )
+    @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6)
 
@@ -568,11 +555,7 @@ def SSE(model):
             },
         }
 
-    @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
-    @unittest.skipIf(
-        not parmest.inverse_reduced_hessian_available,
-        "Cannot test covariance matrix: required ASL dependency is missing",
-    )
+    @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def check_rooney_biegler_results(self, objval, cov):
 
         # get indices in covariance matrix
@@ -596,6 +579,7 @@ def check_rooney_biegler_results(self, objval, cov):
             cov.iloc[rate_constant_index, rate_constant_index], 0.04193591, places=2
         )  # 0.04124 from paper
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics(self):
 
         for model_type, parmest_input in self.input.items():
@@ -609,6 +593,7 @@ def test_parmest_basics(self):
             obj_at_theta = pest.objective_at_theta(parmest_input["theta_vals"])
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_initialize_parmest_model_option(self):
 
         for model_type, parmest_input in self.input.items():
@@ -625,6 +610,7 @@ def test_parmest_basics_with_initialize_parmest_model_option(self):
 
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_square_problem_solve(self):
 
         for model_type, parmest_input in self.input.items():
@@ -641,6 +627,7 @@ def test_parmest_basics_with_square_problem_solve(self):
 
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_square_problem_solve_no_theta_vals(self):
 
         for model_type, parmest_input in self.input.items():
@@ -923,6 +910,7 @@ def test_return_continuous_set_multiple_datasets(self):
         self.assertAlmostEqual(return_vals1["time"].loc[1][18], 2.368, places=3)
         self.assertAlmostEqual(return_vals2["time"].loc[1][18], 2.368, places=3)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_covariance(self):
         from pyomo.contrib.interior_point.inverse_reduced_hessian import (
             inv_reduced_hessian_barrier,
@@ -1217,17 +1205,13 @@ def test_parallel_parmest(self):
             retcode = subprocess.call(rlist)
         assert retcode == 0
 
-    @unittest.skip("Most folks don't have k_aug installed")
+    @unittest.skipUnless(k_aug_available, "k_aug solver not found")
     def test_theta_k_aug_for_Hessian(self):
         # this will fail if k_aug is not installed
         objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug")
         self.assertAlmostEqual(objval, 4.4675, places=2)
 
-    @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
-    @unittest.skipIf(
-        not parmest.inverse_reduced_hessian_available,
-        "Cannot test covariance matrix: required ASL dependency is missing",
-    )
+    @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6)
 
@@ -1485,11 +1469,7 @@ def SSE(model, data):
             },
         }
 
-    @unittest.skipIf(not pynumero_ASL_available, "pynumero ASL is not available")
-    @unittest.skipIf(
-        not parmest.inverse_reduced_hessian_available,
-        "Cannot test covariance matrix: required ASL dependency is missing",
-    )
+    @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def test_parmest_basics(self):
         for model_type, parmest_input in self.input.items():
             pest = parmest.Estimator(
@@ -1518,6 +1498,7 @@ def test_parmest_basics(self):
             obj_at_theta = pest.objective_at_theta(parmest_input["theta_vals"])
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_initialize_parmest_model_option(self):
         for model_type, parmest_input in self.input.items():
             pest = parmest.Estimator(
@@ -1549,6 +1530,7 @@ def test_parmest_basics_with_initialize_parmest_model_option(self):
 
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_square_problem_solve(self):
         for model_type, parmest_input in self.input.items():
             pest = parmest.Estimator(
@@ -1580,6 +1562,7 @@ def test_parmest_basics_with_square_problem_solve(self):
 
             self.assertAlmostEqual(obj_at_theta["obj"][0], 16.531953, places=2)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_parmest_basics_with_square_problem_solve_no_theta_vals(self):
         for model_type, parmest_input in self.input.items():
             pest = parmest.Estimator(
@@ -1923,6 +1906,7 @@ def test_return_continuous_set_multiple_datasets(self):
         self.assertAlmostEqual(return_vals1["time"].loc[1][18], 2.368, places=3)
         self.assertAlmostEqual(return_vals2["time"].loc[1][18], 2.368, places=3)
 
+    @unittest.skipUnless(pynumero_ASL_available, 'pynumero_ASL is not available')
     def test_covariance(self):
         from pyomo.contrib.interior_point.inverse_reduced_hessian import (
             inv_reduced_hessian_barrier,

From f4aeeedd902d0cd1c33108983ecced08a44d461b Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 21 Jun 2024 10:46:11 -0600
Subject: [PATCH 21/69] Relax baseline comparisons

---
 pyomo/contrib/trustregion/tests/test_interface.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/trustregion/tests/test_interface.py b/pyomo/contrib/trustregion/tests/test_interface.py
index 0922ccf950b..d241576f3ba 100644
--- a/pyomo/contrib/trustregion/tests/test_interface.py
+++ b/pyomo/contrib/trustregion/tests/test_interface.py
@@ -234,7 +234,7 @@ def test_updateSurrogateModel(self):
         for key, val in self.interface.data.grad_basis_model_output.items():
             self.assertEqual(value(val), 0)
         for key, val in self.interface.data.truth_model_output.items():
-            self.assertEqual(value(val), 0.8414709848078965)
+            self.assertAlmostEqual(value(val), 0.8414709848078965)
         # The truth gradients should equal the output of [cos(2-1), -cos(2-1)]
         truth_grads = []
         for key, val in self.interface.data.grad_truth_model_output.items():
@@ -332,7 +332,7 @@ def test_calculateFeasibility(self):
         # Check after a solve is completed
         self.interface.data.basis_constraint.activate()
         objective, step_norm, feasibility = self.interface.solveModel()
-        self.assertEqual(feasibility, 0.09569982275514467)
+        self.assertAlmostEqual(feasibility, 0.09569982275514467)
         self.interface.data.basis_constraint.deactivate()
 
     @unittest.skipIf(
@@ -361,7 +361,7 @@ def test_calculateStepSizeInfNorm(self):
         # Check after a solve is completed
         self.interface.data.basis_constraint.activate()
         objective, step_norm, feasibility = self.interface.solveModel()
-        self.assertEqual(step_norm, 3.393437471478297)
+        self.assertAlmostEqual(step_norm, 3.393437471478297)
         self.interface.data.basis_constraint.deactivate()
 
     @unittest.skipIf(

From a1a75b855e151f2d03cef9b2343c260607d28c24 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 21 Jun 2024 10:46:48 -0600
Subject: [PATCH 22/69] Update compare_baseline to use both ABS an REL
 tolerance

---
 pyomo/common/unittest.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/common/unittest.py b/pyomo/common/unittest.py
index 84d962eb784..f1fa01ac486 100644
--- a/pyomo/common/unittest.py
+++ b/pyomo/common/unittest.py
@@ -837,7 +837,7 @@ def filter_file_contents(self, lines, abstol=None):
 
         return filtered
 
-    def compare_baseline(self, test_output, baseline, abstol=1e-6, reltol=None):
+    def compare_baseline(self, test_output, baseline, abstol=1e-6, reltol=1e-8):
         # Filter files independently and then compare filtered contents
         out_filtered = self.filter_file_contents(
             test_output.strip().split('\n'), abstol

From 390463936dbcc2496b694059cad499095a344fea Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 14:21:14 -0600
Subject: [PATCH 23/69] Rework Set initialization to improve performance on
 large sets

---
 pyomo/core/base/set.py            | 423 +++++++++++++++---------------
 pyomo/core/tests/unit/test_set.py |  44 ++--
 2 files changed, 233 insertions(+), 234 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index 8b7c2a246d6..fcabf6674a4 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -16,15 +16,18 @@
 import math
 import sys
 import weakref
-from pyomo.common.pyomo_typing import overload
-from typing import Union, Type, Any as typingAny
+
 from collections.abc import Iterator
+from functools import partial
+from typing import Union, Type, Any as typingAny
 
+from pyomo.common.autoslots import AutoSlots
 from pyomo.common.collections import ComponentSet
 from pyomo.common.deprecation import deprecated, deprecation_warning, RenamedClass
 from pyomo.common.errors import DeveloperError, PyomoException
 from pyomo.common.log import is_debug_set
 from pyomo.common.modeling import NOTSET
+from pyomo.common.pyomo_typing import overload
 from pyomo.common.sorting import sorted_robust
 from pyomo.common.timing import ConstructionTimer
 
@@ -478,9 +481,7 @@ def __call__(self, parent, index):
 
         if not isinstance(_val, Sequence):
             _val = tuple(_val)
-        if len(_val) == 0:
-            return _val
-        if isinstance(_val[0], tuple):
+        if not _val or isinstance(_val[0], tuple):
             return _val
         return self._tuplize(_val, parent, index)
 
@@ -501,7 +502,7 @@ def _tuplize(self, _val, parent, index):
                 "length %s is not a multiple of dimen=%s" % (len(_val), d)
             )
 
-        return list(tuple(_val[d * i : d * (i + 1)]) for i in range(len(_val) // d))
+        return (tuple(_val[i : i + d]) for i in range(0, len(_val), d))
 
 
 class _NotFound(object):
@@ -1364,87 +1365,10 @@ def filter(self):
         return self._filter
 
     def add(self, *values):
-        count = 0
-        _block = self.parent_block()
-        for value in values:
-            if normalize_index.flatten:
-                _value = normalize_index(value)
-                if _value.__class__ is tuple:
-                    _d = len(_value)
-                else:
-                    _d = 1
-            else:
-                # If we are not normalizing indices, then we cannot reliably
-                # infer the set dimen
-                _d = 1
-                if isinstance(value, Sequence) and self.dimen != 1:
-                    _d = len(value)
-                _value = value
-            if _value not in self._domain:
-                raise ValueError(
-                    "Cannot add value %s to Set %s.\n"
-                    "\tThe value is not in the domain %s"
-                    % (value, self.name, self._domain)
-                )
-
-            # We wrap this check in a try-except because some values
-            #  (like lists) are not hashable and can raise exceptions.
-            try:
-                if _value in self:
-                    logger.warning(
-                        "Element %s already exists in Set %s; no action taken"
-                        % (value, self.name)
-                    )
-                    continue
-            except:
-                exc = sys.exc_info()
-                raise TypeError(
-                    "Unable to insert '%s' into Set %s:\n\t%s: %s"
-                    % (value, self.name, exc[0].__name__, exc[1])
-                )
-
-            if self._filter is not None:
-                if not self._filter(_block, _value):
-                    continue
-
-            if self._validate is not None:
-                try:
-                    flag = self._validate(_block, _value)
-                except:
-                    logger.error(
-                        "Exception raised while validating element '%s' "
-                        "for Set %s" % (value, self.name)
-                    )
-                    raise
-                if not flag:
-                    raise ValueError(
-                        "The value=%s violates the validation rule of Set %s"
-                        % (value, self.name)
-                    )
-
-            # If the Set has a fixed dimension, check that this element is
-            # compatible.
-            if self._dimen is not None:
-                if _d != self._dimen:
-                    if self._dimen is UnknownSetDimen:
-                        # The first thing added to a Set with unknown
-                        # dimension sets its dimension
-                        self._dimen = _d
-                    else:
-                        raise ValueError(
-                            "The value=%s has dimension %s and is not "
-                            "valid for Set %s which has dimen=%s"
-                            % (value, _d, self.name, self._dimen)
-                        )
-
-            # Add the value to this object (this last redirection allows
-            # derived classes to implement a different storage mechanism)
-            self._add_impl(_value)
-            count += 1
-        return count
+        return self.update(values)
 
-    def _add_impl(self, value):
-        self._values.add(value)
+    def _update_impl(self, values):
+        self._values.update(values)
 
     def remove(self, val):
         self._values.remove(val)
@@ -1457,17 +1381,144 @@ def clear(self):
 
     def set_value(self, val):
         self.clear()
-        for x in val:
-            self.add(x)
+        self.update(val)
 
     def update(self, values):
-        for v in values:
-            if v not in self:
-                self.add(v)
+        # _values was initialized above...
+        #
+        # Special case: set operations that are not first attached
+        # to the model must be constructed.
+        if isinstance(values, SetOperator):
+            values.construct()
+        try:
+            val_iter = iter(values)
+        except TypeError:
+            logger.error(
+                "Initializer for Set %s%s returned non-iterable object "
+                "of type %s."
+                % (
+                    self.name,
+                    ("[%s]" % (index,) if self.is_indexed() else ""),
+                    (values if values.__class__ is type else type(values).__name__),
+                )
+            )
+            raise
+
+        if self._dimen is not None:
+            if normalize_index.flatten:
+                val_iter = self._cb_normalized_dimen_verifier(self._dimen, val_iter)
+            else:
+                val_iter = self._cb_raw_dimen_verifier(self._dimen, val_iter)
+        elif normalize_index.flatten:
+            val_iter = map(normalize_index, val_iter)
+        else:
+            val_iter = self._cb_check_set_end(val_iter)
+
+        if self._domain is not Any:
+            val_iter = self._cb_domain_verifier(self._domain, val_iter)
+
+        if self._filter is not None:
+            val_iter = filter(partial(self._filter, self.parent_block()), val_iter)
+
+        if self._validate is not None:
+            val_iter = self._cb_validate(self._validate, self.parent_block(), val_iter)
+
+        nOrig = len(self)
+        # We wrap this check in a try-except because some values
+        #  (like lists) are not hashable and can raise exceptions.
+        try:
+            self._update_impl(val_iter)
+        except Set.End:
+            pass
+        return len(self) - nOrig
 
     def pop(self):
         return self._values.pop()
 
+    def _cb_domain_verifier(self, domain, val_iter):
+        for value in val_iter:
+            if value not in domain:
+                raise ValueError(
+                    "Cannot add value %s to Set %s.\n"
+                    "\tThe value is not in the domain %s"
+                    % (value, self.name, self._domain)
+                )
+            yield value
+
+    def _cb_check_set_end(self, val_iter):
+        for value in val_iter:
+            if value is Set.End:
+                return
+            yield value
+
+    def _cb_validate(self, validate, block, val_iter):
+        for value in val_iter:
+            try:
+                flag = validate(block, value)
+            except:
+                logger.error(
+                    "Exception raised while validating element '%s' "
+                    "for Set %s" % (value, self.name)
+                )
+                raise
+            if not flag:
+                raise ValueError(
+                    "The value=%s violates the validation rule of Set %s"
+                    % (value, self.name)
+                )
+            yield value
+
+    def _cb_normalized_dimen_verifier(self, dimen, val_iter):
+        # It is important that the iterator is an actual iterator
+        val_iter = iter(val_iter)
+        for value in val_iter:
+            if value.__class__ is tuple:
+                if dimen == len(value):
+                    yield value[0] if dimen == 1 else value
+                    continue
+            elif dimen == 1 and value.__class__ in native_types:
+                yield value
+                continue
+
+            # Note: normalize_index() will never return a 1-tuple
+            normalized_value = normalize_index(value)
+            _d = len(normalized_value) if normalized_value.__class__ is tuple else 1
+            if _d == dimen:
+                yield normalized_value
+            elif dimen is UnknownSetDimen:
+                # The first thing added to a Set with unknown dimension
+                # sets its dimension
+                self._dimen = dimen = _d
+                yield normalized_value
+            else:
+                raise ValueError(
+                    "The value=%s has dimension %s and is not "
+                    "valid for Set %s which has dimen=%s"
+                    % (value, _d, self.name, self._dimen)
+                )
+
+    def _cb_raw_dimen_verifier(self, dimen, val_iter):
+        for value in val_iter:
+            if isinstance(value, Sequence):
+                if dimen == len(value):
+                    yield value
+                    continue
+            elif dimen == 1:
+                yield value
+                continue
+            _d = len(value) if isinstance(value, Sequence) else 1
+            if dimen is UnknownSetDimen:
+                # The first thing added to a Set with unknown dimension
+                # sets its dimension
+                self._dimen = dimen = _d
+                yield value
+            else:
+                raise ValueError(
+                    "The value=%s has dimension %s and is not "
+                    "valid for Set %s which has dimen=%s"
+                    % (value, _d, self.name, self._dimen)
+                )
+
 
 class _FiniteSetData(metaclass=RenamedClass):
     __renamed__new_class__ = FiniteSetData
@@ -1655,27 +1706,30 @@ class OrderedSetData(_OrderedSetMixin, FiniteSetData):
 
     def __init__(self, component):
         self._values = {}
-        self._ordered_values = []
+        self._ordered_values = None
         FiniteSetData.__init__(self, component=component)
 
     def _iter_impl(self):
         """
         Return an iterator for the set.
         """
-        return iter(self._ordered_values)
+        return iter(self._values)
 
     def __reversed__(self):
-        return reversed(self._ordered_values)
+        return reversed(self._values)
 
-    def _add_impl(self, value):
-        self._values[value] = len(self._values)
-        self._ordered_values.append(value)
+    def _update_impl(self, values):
+        for val in values:
+            # Note that we reset _ordered_values within the loop because
+            # of an old example where the initializer rule makes
+            # reference to values previously inserted into the Set
+            # (which triggered the creation of the _ordered_values)
+            self._ordered_values = None
+            self._values[val] = None
 
     def remove(self, val):
-        idx = self._values.pop(val)
-        self._ordered_values.pop(idx)
-        for i in range(idx, len(self._ordered_values)):
-            self._values[self._ordered_values[i]] -= 1
+        self._values.pop(val)
+        self._ordered_values = None
 
     def discard(self, val):
         try:
@@ -1685,7 +1739,7 @@ def discard(self, val):
 
     def clear(self):
         self._values.clear()
-        self._ordered_values = []
+        self._ordered_values = None
 
     def pop(self):
         try:
@@ -1704,6 +1758,8 @@ def at(self, index):
         The public Set API is 1-based, even though the
         internal _lookup and _values are (pythonically) 0-based.
         """
+        if self._ordered_values is None:
+            self._rebuild_ordered_values()
         i = self._to_0_based_index(index)
         try:
             return self._ordered_values[i]
@@ -1723,6 +1779,8 @@ def ord(self, item):
         # when they are actually put as Set members.  So, we will look
         # for the exact thing that the user sent us and then fall back
         # on the scalar.
+        if self._ordered_values is None:
+            self._rebuild_ordered_values()
         try:
             return self._values[item] + 1
         except KeyError:
@@ -1733,6 +1791,12 @@ def ord(self, item):
         except KeyError:
             raise ValueError("%s.ord(x): x not in %s" % (self.name, self.name))
 
+    def _rebuild_ordered_values(self):
+        _set = self._values
+        self._ordered_values = list(_set)
+        for i, v in enumerate(self._ordered_values):
+            _set[v] = i
+
 
 class _OrderedSetData(metaclass=RenamedClass):
     __renamed__new_class__ = OrderedSetData
@@ -1760,7 +1824,8 @@ def set_value(self, val):
                 "This WILL potentially lead to nondeterministic behavior "
                 "in Pyomo" % (type(val).__name__,)
             )
-        super(InsertionOrderSetData, self).set_value(val)
+        self.clear()
+        super().update(val)
 
     def update(self, values):
         if type(values) in Set._UnorderedInitializers:
@@ -1770,7 +1835,7 @@ def update(self, values):
                 "This WILL potentially lead to nondeterministic behavior "
                 "in Pyomo" % (type(values).__name__,)
             )
-        super(InsertionOrderSetData, self).update(values)
+        return super().update(values)
 
 
 class _InsertionOrderSetData(metaclass=RenamedClass):
@@ -1800,73 +1865,38 @@ class SortedSetData(_SortedSetMixin, OrderedSetData):
     Public Class Attributes:
     """
 
-    __slots__ = ('_is_sorted',)
-
-    def __init__(self, component):
-        # An empty set is sorted...
-        self._is_sorted = True
-        OrderedSetData.__init__(self, component=component)
+    __slots__ = ()
 
     def _iter_impl(self):
         """
         Return an iterator for the set.
         """
-        if not self._is_sorted:
-            self._sort()
-        return super(SortedSetData, self)._iter_impl()
+        if self._ordered_values is None:
+            self._rebuild_ordered_values()
+        return iter(self._ordered_values)
 
     def __reversed__(self):
-        if not self._is_sorted:
-            self._sort()
-        return super(SortedSetData, self).__reversed__()
+        if self._ordered_values is None:
+            self._rebuild_ordered_values()
+        return reversed(self._ordered_values)
 
-    def _add_impl(self, value):
-        # Note that the sorted status has no bearing on insertion,
-        # so there is no reason to check if the data is correctly sorted
-        self._values[value] = len(self._values)
-        self._ordered_values.append(value)
-        self._is_sorted = False
+    def _update_impl(self, values):
+        for val in values:
+            self._values[val] = None
+        self._ordered_values = None
 
     # Note: removing data does not affect the sorted flag
     # def remove(self, val):
     # def discard(self, val):
 
-    def clear(self):
-        super(SortedSetData, self).clear()
-        self._is_sorted = True
-
-    def at(self, index):
-        """
-        Return the specified member of the set.
-
-        The public Set API is 1-based, even though the
-        internal _lookup and _values are (pythonically) 0-based.
-        """
-        if not self._is_sorted:
-            self._sort()
-        return super(SortedSetData, self).at(index)
-
-    def ord(self, item):
-        """
-        Return the position index of the input value.
-
-        Note that Pyomo Set objects have positions starting at 1 (not 0).
-
-        If the search item is not in the Set, then an IndexError is raised.
-        """
-        if not self._is_sorted:
-            self._sort()
-        return super(SortedSetData, self).ord(item)
-
     def sorted_data(self):
         return self.data()
 
-    def _sort(self):
-        self._ordered_values = list(
-            self.parent_component()._sort_fcn(self._ordered_values)
-        )
-        self._values = {j: i for i, j in enumerate(self._ordered_values)}
-        self._is_sorted = True
+    def _rebuild_ordered_values(self):
+        _set = self._values
+        self._ordered_values = list(self.parent_component()._sort_fcn(_set))
+        for i, v in enumerate(self._ordered_values):
+            _set[v] = i
 
 
 class _SortedSetData(metaclass=RenamedClass):
@@ -1974,10 +2004,11 @@ class Set(IndexedComponent):
 
     """
 
-    class End(object):
-        pass
+    class SetEndType(type):
+        def __hash__(self):
+            raise Set.End()
 
-    class Skip(object):
+    class End(Exception, metaclass=SetEndType):
         pass
 
     class InsertionOrder(object):
@@ -1988,6 +2019,7 @@ class SortedOrder(object):
 
     _ValidOrderedAuguments = {True, False, InsertionOrder, SortedOrder}
     _UnorderedInitializers = {set}
+    _SetEndEncountered = False
 
     @overload
     def __new__(cls: Type[Set], *args, **kwds) -> Union[SetData, IndexedSet]: ...
@@ -2222,21 +2254,6 @@ def _getitem_when_not_present(self, index):
         if _d is UnknownSetDimen and domain is not None and domain.dimen is not None:
             _d = domain.dimen
 
-        if self._init_values is not None:
-            self._init_values._dimen = _d
-            try:
-                _values = self._init_values(_block, index)
-            except TuplizeError as e:
-                raise ValueError(
-                    str(e) % (self._name, "[%s]" % index if self.is_indexed() else "")
-                )
-
-            if _values is Set.Skip:
-                return
-            elif _values is None:
-                raise ValueError(
-                    "Set rule or initializer returned None instead of Set.Skip"
-                )
         if index is None and not self.is_indexed():
             obj = self._data[index] = self
         else:
@@ -2258,55 +2275,35 @@ def _getitem_when_not_present(self, index):
                 obj._validate = self._init_validate
         if self._init_filter is not None:
             try:
-                _filter = Initializer(self._init_filter(_block, index))
-                if _filter.constant():
+                obj._filter = Initializer(self._init_filter(_block, index))
+                if obj._filter.constant():
                     # _init_filter was the actual filter function; use it.
-                    _filter = self._init_filter
+                    obj._filter = self._init_filter
             except:
                 # We will assume any exceptions raised when getting the
                 # filter for this index indicate that the function
                 # should have been passed directly to the underlying sets.
-                _filter = self._init_filter
+                obj._filter = self._init_filter
         else:
-            _filter = None
+            obj._filter = None
         if self._init_values is not None:
-            # _values was initialized above...
-            if obj.isordered() and type(_values) in Set._UnorderedInitializers:
-                logger.warning(
-                    "Initializing ordered Set %s with a fundamentally "
-                    "unordered data source (type: %s).  This WILL potentially "
-                    "lead to nondeterministic behavior in Pyomo"
-                    % (self.name, type(_values).__name__)
-                )
-            # Special case: set operations that are not first attached
-            # to the model must be constructed.
-            if isinstance(_values, SetOperator):
-                _values.construct()
+            # record the user-provided dimen in the initializer
+            self._init_values._dimen = _d
             try:
-                val_iter = iter(_values)
-            except TypeError:
-                logger.error(
-                    "Initializer for Set %s%s returned non-iterable object "
-                    "of type %s."
-                    % (
-                        self.name,
-                        ("[%s]" % (index,) if self.is_indexed() else ""),
-                        (
-                            _values
-                            if _values.__class__ is type
-                            else type(_values).__name__
-                        ),
-                    )
+                _values = self._init_values(_block, index)
+            except TuplizeError as e:
+                raise ValueError(
+                    str(e) % (self._name, "[%s]" % index if self.is_indexed() else "")
                 )
-                raise
-            for val in val_iter:
-                if val is Set.End:
-                    break
-                if _filter is None or _filter(_block, val):
-                    obj.add(val)
-        # We defer adding the filter until now so that add() doesn't
-        # call it a second time.
-        obj._filter = _filter
+            if _values is Set.Skip:
+                del self._data[index]
+                return
+            elif _values is None:
+                raise ValueError(
+                    "Set rule or initializer returned None instead of Set.Skip"
+                )
+
+            obj.set_value(_values)
         return obj
 
     @staticmethod
diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index f62589a6873..b5e47975c00 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -3763,8 +3763,8 @@ def I_init(m):
             m = ConcreteModel()
             m.I = Set(initialize={1, 3, 2, 4})
             ref = (
-                "Initializing ordered Set I with a "
-                "fundamentally unordered data source (type: set)."
+                'Calling set_value() on an insertion order Set with a fundamentally '
+                'unordered data source (type: set).'
             )
             self.assertIn(ref, output.getvalue())
         self.assertEqual(m.I.sorted_data(), (1, 2, 3, 4))
@@ -3877,12 +3877,13 @@ def _verify(_s, _l):
         m.I.add(4)
         _verify(m.I, [1, 3, 2, 4])
 
+        N = len(m.I)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             m.I.add(3)
-        self.assertEqual(
-            output.getvalue(), "Element 3 already exists in Set I; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
+        self.assertEqual(N, len(m.I))
         _verify(m.I, [1, 3, 2, 4])
 
         m.I.remove(3)
@@ -3959,12 +3960,13 @@ def _verify(_s, _l):
         m.I.add(4)
         _verify(m.I, [1, 2, 3, 4])
 
+        N = len(m.I)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             m.I.add(3)
-        self.assertEqual(
-            output.getvalue(), "Element 3 already exists in Set I; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
+        self.assertEqual(N, len(m.I))
         _verify(m.I, [1, 2, 3, 4])
 
         m.I.remove(3)
@@ -4052,12 +4054,13 @@ def _verify(_s, _l):
         m.I.add(4)
         _verify(m.I, [1, 2, 3, 4])
 
+        N = len(m.I)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             m.I.add(3)
-        self.assertEqual(
-            output.getvalue(), "Element 3 already exists in Set I; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
+        self.assertEqual(N, len(m.I))
         _verify(m.I, [1, 2, 3, 4])
 
         m.I.remove(3)
@@ -4248,24 +4251,23 @@ def test_add_filter_validate(self):
         self.assertIn(1, m.I)
         self.assertIn(1.0, m.I)
 
+        N = len(m.I)
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             self.assertFalse(m.I.add(1))
-        self.assertEqual(
-            output.getvalue(), "Element 1 already exists in Set I; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
+        self.assertEqual(N, len(m.I))
 
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             self.assertFalse(m.I.add((1,)))
-        self.assertEqual(
-            output.getvalue(), "Element (1,) already exists in Set I; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
 
         m.J = Set()
         # Note that pypy raises a different exception from cpython
         err = (
-            r"Unable to insert '{}' into Set J:\n\tTypeError: "
             r"((unhashable type: 'dict')|('dict' objects are unhashable))"
         )
         with self.assertRaisesRegex(TypeError, err):
@@ -4275,9 +4277,9 @@ def test_add_filter_validate(self):
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             self.assertFalse(m.J.add(1))
-        self.assertEqual(
-            output.getvalue(), "Element 1 already exists in Set J; no action taken\n"
-        )
+        # In Pyomo <= 6.7.3 duplicate values logged a warning.
+        self.assertEqual(output.getvalue(), "")
+        self.assertEqual(N, len(m.I))
 
         def _l_tri(model, i, j):
             self.assertIs(model, m)

From 5eb9875fe9d0e940fb914d867dc146c2103e8f64 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 14:28:20 -0600
Subject: [PATCH 24/69] implement Set.first() and .last() without forcing
 _ordered_values construction

---
 pyomo/core/base/set.py | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index fcabf6674a4..e6444c48452 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1596,10 +1596,16 @@ def ordered_iter(self):
         return iter(self)
 
     def first(self):
-        return self.at(1)
+        try:
+            return next(iter(self))
+        except StopIteration:
+            raise IndexError(f"{self.name} index out of range") from None
 
     def last(self):
-        return self.at(len(self))
+        try:
+            return next(reversed(self))
+        except StopIteration:
+            raise IndexError(f"{self.name} index out of range") from None
 
     def next(self, item, step=1):
         """
@@ -1745,9 +1751,9 @@ def pop(self):
         try:
             ans = self.last()
         except IndexError:
-            # Map the index error to a KeyError for consistency with
-            # set().pop()
-            raise KeyError('pop from an empty set')
+            # Map the exception for iterating over an empty dict to a
+            # KeyError for consistency with set().pop()
+            raise KeyError('pop from an empty set') from None
         self.discard(ans)
         return ans
 

From 4716f9d36bb6c54b3b7b448e7e9976f66f351ec8 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 15:42:19 -0600
Subject: [PATCH 25/69] NFC: apply black

---
 pyomo/core/tests/unit/test_set.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index b5e47975c00..6fb95a31a59 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -4267,9 +4267,7 @@ def test_add_filter_validate(self):
 
         m.J = Set()
         # Note that pypy raises a different exception from cpython
-        err = (
-            r"((unhashable type: 'dict')|('dict' objects are unhashable))"
-        )
+        err = r"((unhashable type: 'dict')|('dict' objects are unhashable))"
         with self.assertRaisesRegex(TypeError, err):
             m.J.add({})
 

From b6cb09cde03440c599946c376d63abaf4f144b1a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 15:48:26 -0600
Subject: [PATCH 26/69] Restore original ordered set initialization warning

---
 pyomo/core/base/set.py            | 15 ++++++++++++++-
 pyomo/core/tests/unit/test_set.py |  2 +-
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index e6444c48452..040c66ffd3a 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1383,6 +1383,9 @@ def set_value(self, val):
         self.clear()
         self.update(val)
 
+    def _initialize(self, val):
+        self.update(val)
+
     def update(self, values):
         # _values was initialized above...
         #
@@ -1822,6 +1825,16 @@ class InsertionOrderSetData(OrderedSetData):
 
     __slots__ = ()
 
+    def _initialize(self, val):
+        if type(val) in Set._UnorderedInitializers:
+            logger.warning(
+                "Initializing ordered Set %s with "
+                "a fundamentally unordered data source (type: %s).  "
+                "This WILL potentially lead to nondeterministic behavior "
+                "in Pyomo" % (self.name, type(val).__name__)
+            )
+        super().update(val)
+
     def set_value(self, val):
         if type(val) in Set._UnorderedInitializers:
             logger.warning(
@@ -2309,7 +2322,7 @@ def _getitem_when_not_present(self, index):
                     "Set rule or initializer returned None instead of Set.Skip"
                 )
 
-            obj.set_value(_values)
+            obj._initialize(_values)
         return obj
 
     @staticmethod
diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index 6fb95a31a59..ff0aaa5600b 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -3763,7 +3763,7 @@ def I_init(m):
             m = ConcreteModel()
             m.I = Set(initialize={1, 3, 2, 4})
             ref = (
-                'Calling set_value() on an insertion order Set with a fundamentally '
+                'Initializing ordered Set I with a fundamentally '
                 'unordered data source (type: set).'
             )
             self.assertIn(ref, output.getvalue())

From d5e4f05f58d4c25f831130040b4676b6f48555fb Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 16:54:25 -0600
Subject: [PATCH 27/69] Restore previous Set add()/update() return values

---
 pyomo/core/base/set.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index 040c66ffd3a..f183d906607 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1365,7 +1365,9 @@ def filter(self):
         return self._filter
 
     def add(self, *values):
-        return self.update(values)
+        N = len(self)
+        self.update(values)
+        return len(self) - N
 
     def _update_impl(self, values):
         self._values.update(values)
@@ -1426,14 +1428,12 @@ def update(self, values):
         if self._validate is not None:
             val_iter = self._cb_validate(self._validate, self.parent_block(), val_iter)
 
-        nOrig = len(self)
         # We wrap this check in a try-except because some values
         #  (like lists) are not hashable and can raise exceptions.
         try:
             self._update_impl(val_iter)
         except Set.End:
             pass
-        return len(self) - nOrig
 
     def pop(self):
         return self._values.pop()
@@ -1854,7 +1854,7 @@ def update(self, values):
                 "This WILL potentially lead to nondeterministic behavior "
                 "in Pyomo" % (type(values).__name__,)
             )
-        return super().update(values)
+        super().update(values)
 
 
 class _InsertionOrderSetData(metaclass=RenamedClass):

From d0338bc2e960aa12fbd4877cba87d4b7081de1b4 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 17:16:42 -0600
Subject: [PATCH 28/69] Rework Set.End to avoid conflict with
 dask.multiprocessing

(which attempts to hash all known exception types)
---
 pyomo/core/base/set.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index f183d906607..a691973a1a7 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1432,7 +1432,7 @@ def update(self, values):
         #  (like lists) are not hashable and can raise exceptions.
         try:
             self._update_impl(val_iter)
-        except Set.End:
+        except Set._SetEndException:
             pass
 
     def pop(self):
@@ -2023,11 +2023,14 @@ class Set(IndexedComponent):
 
     """
 
-    class SetEndType(type):
+    class _SetEndException(Exception):
+        pass
+
+    class _SetEndType(type):
         def __hash__(self):
-            raise Set.End()
+            raise Set._SetEndException()
 
-    class End(Exception, metaclass=SetEndType):
+    class End(metaclass=_SetEndType):
         pass
 
     class InsertionOrder(object):

From 1ad9d4eebdf34180e1045bbf1c20d6e097958217 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 25 Jun 2024 17:18:06 -0600
Subject: [PATCH 29/69] Remove unused class attribute

---
 pyomo/core/base/set.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index a691973a1a7..1ae866e84d7 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -2041,7 +2041,6 @@ class SortedOrder(object):
 
     _ValidOrderedAuguments = {True, False, InsertionOrder, SortedOrder}
     _UnorderedInitializers = {set}
-    _SetEndEncountered = False
 
     @overload
     def __new__(cls: Type[Set], *args, **kwds) -> Union[SetData, IndexedSet]: ...

From f658981cee62da884f60347ab1341d709d387d57 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 9 Jul 2024 12:54:52 -0600
Subject: [PATCH 30/69] NFC: update spaces in message

---
 pyomo/core/base/constraint.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index 3ee5a82ef58..defaea99dff 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -347,8 +347,8 @@ def set_value(self, expr):
             if getattr(expr, 'strict', False) in _strict_relational_exprs:
                 raise ValueError(
                     "Constraint '%s' encountered a strict "
-                    "inequality expression ('>' or '< '). All"
-                    " constraints must be formulated using "
+                    "inequality expression ('>' or '<').  All "
+                    "constraints must be formulated using "
                     "using '<=', '>=', or '=='." % (self.name,)
                 )
             self._expr = expr

From 533f6165bd7b9245ae3f69e1368dd6e54e101d1a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Wed, 10 Jul 2024 10:59:00 -0600
Subject: [PATCH 31/69] Fix (and test) FBBT error with RangedExpressions

---
 pyomo/contrib/fbbt/fbbt.py            |  2 +-
 pyomo/contrib/fbbt/tests/test_fbbt.py | 23 +++++++++++++++++++++++
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 60ac0603388..1a0b2992b07 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -1389,7 +1389,7 @@ def _fbbt_block(m, config):
     for c in m.component_data_objects(
         ctype=Constraint, active=True, descend_into=config.descend_into, sort=True
     ):
-        for v in identify_variables(c.body):
+        for v in identify_variables(c.expr):
             if v not in var_to_con_map:
                 var_to_con_map[v] = list()
             if v.lb is None:
diff --git a/pyomo/contrib/fbbt/tests/test_fbbt.py b/pyomo/contrib/fbbt/tests/test_fbbt.py
index f7d08d11215..033d3ac7a78 100644
--- a/pyomo/contrib/fbbt/tests/test_fbbt.py
+++ b/pyomo/contrib/fbbt/tests/test_fbbt.py
@@ -1331,6 +1331,29 @@ def test_named_expr(self):
         self.assertAlmostEqual(m.x.lb, 2)
         self.assertAlmostEqual(m.x.ub, 3)
 
+    def test_ranged_expression(self):
+        m = pyo.ConcreteModel()
+        m.l = pyo.Var(bounds=(2, None))
+        m.x = pyo.Var()
+        m.u = pyo.Var(bounds=(None, 8))
+        m.c = pyo.Constraint(expr=pyo.inequality(m.l, m.x, m.u))
+        self.tightener(m)
+        self.tightener(m)
+        self.assertEqual(m.l.bounds, (2, 8))
+        self.assertEqual(m.x.bounds, (2, 8))
+        self.assertEqual(m.u.bounds, (2, 8))
+
+        m = pyo.ConcreteModel()
+        m.l = pyo.Var(bounds=(2, None))
+        m.x = pyo.Var(bounds=(3, 7))
+        m.u = pyo.Var(bounds=(None, 8))
+        m.c = pyo.Constraint(expr=pyo.inequality(m.l, m.x, m.u))
+        self.tightener(m)
+        self.tightener(m)
+        self.assertEqual(m.l.bounds, (2, 7))
+        self.assertEqual(m.x.bounds, (3, 7))
+        self.assertEqual(m.u.bounds, (3, 8))
+
 
 class TestFBBT(FbbtTestBase, unittest.TestCase):
     def setUp(self) -> None:

From 399c7311afaebcfbf1fb8753a477ee6f82717f9e Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Wed, 10 Jul 2024 12:11:33 -0600
Subject: [PATCH 32/69] Only test variable RangedExpressions in Python FBBT

---
 pyomo/contrib/fbbt/tests/test_fbbt.py | 15 ++++++++++-----
 1 file changed, 10 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/fbbt/tests/test_fbbt.py b/pyomo/contrib/fbbt/tests/test_fbbt.py
index 033d3ac7a78..ff1cc8a5cfb 100644
--- a/pyomo/contrib/fbbt/tests/test_fbbt.py
+++ b/pyomo/contrib/fbbt/tests/test_fbbt.py
@@ -1331,7 +1331,17 @@ def test_named_expr(self):
         self.assertAlmostEqual(m.x.lb, 2)
         self.assertAlmostEqual(m.x.ub, 3)
 
+
+class TestFBBT(FbbtTestBase, unittest.TestCase):
+    def setUp(self) -> None:
+        self.tightener = fbbt
+
     def test_ranged_expression(self):
+        # The python version of FBBT is slightly more flexible than
+        # APPSI's cmodel (it allows - and correctly handles -
+        # RangedExpressions with variable lower / upper bounds.  If we
+        # ever port that functionality into APPSI, then this test can be
+        # moved into the base class.
         m = pyo.ConcreteModel()
         m.l = pyo.Var(bounds=(2, None))
         m.x = pyo.Var()
@@ -1353,8 +1363,3 @@ def test_ranged_expression(self):
         self.assertEqual(m.l.bounds, (2, 7))
         self.assertEqual(m.x.bounds, (3, 7))
         self.assertEqual(m.u.bounds, (3, 8))
-
-
-class TestFBBT(FbbtTestBase, unittest.TestCase):
-    def setUp(self) -> None:
-        self.tightener = fbbt

From 7f84d3e9b5e23a0273dc8c55527c74c538a9f436 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Wed, 10 Jul 2024 12:12:28 -0600
Subject: [PATCH 33/69] Move to using expr (not body) when gathering variabels
 from constraints

---
 pyomo/contrib/community_detection/community_graph.py | 2 +-
 pyomo/contrib/fbbt/fbbt.py                           | 6 +++---
 pyomo/util/infeasible.py                             | 2 +-
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/pyomo/contrib/community_detection/community_graph.py b/pyomo/contrib/community_detection/community_graph.py
index 889940b5996..c67a8cd6690 100644
--- a/pyomo/contrib/community_detection/community_graph.py
+++ b/pyomo/contrib/community_detection/community_graph.py
@@ -123,7 +123,7 @@ def generate_model_graph(
         # Create a list of the variable numbers that occur in the given constraint equation
         numbered_variables_in_constraint_equation = [
             component_number_map[constraint_variable]
-            for constraint_variable in identify_variables(model_constraint.body)
+            for constraint_variable in identify_variables(model_constraint.expr)
         ]
 
         # Update constraint_variable_map
diff --git a/pyomo/contrib/fbbt/fbbt.py b/pyomo/contrib/fbbt/fbbt.py
index 1a0b2992b07..4bd0e4552a1 100644
--- a/pyomo/contrib/fbbt/fbbt.py
+++ b/pyomo/contrib/fbbt/fbbt.py
@@ -1576,14 +1576,14 @@ def __init__(self, comp):
         if comp.ctype == Constraint:
             if comp.is_indexed():
                 for c in comp.values():
-                    self._vars.update(identify_variables(c.body))
+                    self._vars.update(identify_variables(c.expr))
             else:
-                self._vars.update(identify_variables(comp.body))
+                self._vars.update(identify_variables(comp.expr))
         else:
             for c in comp.component_data_objects(
                 Constraint, descend_into=True, active=True, sort=True
             ):
-                self._vars.update(identify_variables(c.body))
+                self._vars.update(identify_variables(c.expr))
 
     def save_bounds(self):
         bnds = ComponentMap()
diff --git a/pyomo/util/infeasible.py b/pyomo/util/infeasible.py
index 961d5b35036..6a90a4c3773 100644
--- a/pyomo/util/infeasible.py
+++ b/pyomo/util/infeasible.py
@@ -159,7 +159,7 @@ def log_infeasible_constraints(
         if log_variables:
             line += ''.join(
                 f"\n  - VAR {v.name}: {v.value}"
-                for v in identify_variables(constr.body, include_fixed=True)
+                for v in identify_variables(constr.expr, include_fixed=True)
             )
 
         logger.info(line)

From b2f728d694d9f1ab2ae09ef06a6c47315cc64563 Mon Sep 17 00:00:00 2001
From: Emma Johnson <emma58@users.noreply.github.com>
Date: Wed, 10 Jul 2024 13:35:49 -0600
Subject: [PATCH 34/69] Adding a test for the bug--we don't like setting
 something equivalent to 'True'

---
 .../cp/tests/test_logical_to_disjunctive.py   | 39 ++++++++++++++++++-
 1 file changed, 37 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
index 3f66aa57726..01ef37ada5f 100755
--- a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
+++ b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
@@ -200,6 +200,43 @@ def test_equivalence(self):
 
         assertExpressionsEqual(self, m.cons[10].expr, m.z[5] >= 1)
 
+    def test_equivalent_to_True(self):
+        m = self.make_model()
+        e = m.a.equivalent_to(True)
+
+        visitor = LogicalToDisjunctiveVisitor()
+        m.cons = visitor.constraints
+        m.z = visitor.z_vars
+
+        visitor.walk_expression(e)
+
+        self.assertIs(m.a.get_associated_binary(), m.z[1])
+        self.assertEqual(len(m.z), 4)
+        self.assertEqual(len(m.cons), 10)
+
+        # z[2] == !a v True
+        assertExpressionsEqual(
+            self, m.cons[1].expr, (1 - m.z[2]) + (1 - m.z[1]) + 1 >= 1
+        )
+        assertExpressionsEqual(self, m.cons[2].expr, 1 - (1 - m.z[1]) + m.z[2] >= 1)
+        assertExpressionsEqual(self, m.cons[3].expr, m.z[2] + (1 - 1) >= 1)
+
+        # z[3] == a v ! c
+        assertExpressionsEqual(
+            self, m.cons[4].expr, (1 - m.z[3]) + (1 - 1) + m.z[1] >= 1
+        )
+        assertExpressionsEqual(self, m.cons[5].expr, m.z[3] + (1 - m.z[1]) >= 1)
+        assertExpressionsEqual(self, m.cons[6].expr, 1 - (1 - 1) + m.z[3] >= 1)
+
+        # z[4] == z[2] ^ z[3]
+        assertExpressionsEqual(self, m.cons[7].expr, m.z[4] <= m.z[2])
+        assertExpressionsEqual(self, m.cons[8].expr, m.z[4] <= m.z[3])
+        assertExpressionsEqual(
+            self, m.cons[9].expr, 1 - m.z[4] <= 2 - (m.z[2] + m.z[3])
+        )
+
+        assertExpressionsEqual(self, m.cons[10].expr, m.z[4] >= 1)
+
     def test_xor(self):
         m = self.make_model()
         e = m.a.xor(m.b)
@@ -263,8 +300,6 @@ def test_at_most(self):
         # z3 = a ^ b
         assertExpressionsEqual(self, m.cons[1].expr, m.z[3] <= a)
         assertExpressionsEqual(self, m.cons[2].expr, m.z[3] <= b)
-        m.cons.pprint()
-        print(m.cons[3].expr)
         assertExpressionsEqual(self, m.cons[3].expr, 1 - m.z[3] <= 2 - sum([a, b]))
 
         # atmost in disjunctive form

From 531f758ee4235eeaedae12af8f44d9f8d083e308 Mon Sep 17 00:00:00 2001
From: Emma Johnson <emma58@users.noreply.github.com>
Date: Wed, 10 Jul 2024 13:51:07 -0600
Subject: [PATCH 35/69] Adding handling in beforeChild to cast bools to ints so
 that they will be binary and play nice with algebraic expressions as we build
 them

---
 pyomo/contrib/cp/tests/test_logical_to_disjunctive.py      | 4 ++--
 .../contrib/cp/transform/logical_to_disjunctive_walker.py  | 7 +++++++
 2 files changed, 9 insertions(+), 2 deletions(-)

diff --git a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
index 01ef37ada5f..9c79364c367 100755
--- a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
+++ b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
@@ -223,10 +223,10 @@ def test_equivalent_to_True(self):
 
         # z[3] == a v ! c
         assertExpressionsEqual(
-            self, m.cons[4].expr, (1 - m.z[3]) + (1 - 1) + m.z[1] >= 1
+            self, m.cons[4].expr, (1 - m.z[3]) + m.z[1] >= 1
         )
         assertExpressionsEqual(self, m.cons[5].expr, m.z[3] + (1 - m.z[1]) >= 1)
-        assertExpressionsEqual(self, m.cons[6].expr, 1 - (1 - 1) + m.z[3] >= 1)
+        assertExpressionsEqual(self, m.cons[6].expr, m.z[3] + 1 >= 1)
 
         # z[4] == z[2] ^ z[3]
         assertExpressionsEqual(self, m.cons[7].expr, m.z[4] <= m.z[2])
diff --git a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
index 09894b47090..47a85076ce1 100644
--- a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
+++ b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
@@ -243,6 +243,13 @@ def initializeWalker(self, expr):
         return True, expr
 
     def beforeChild(self, node, child, child_idx):
+        if child.__class__ is bool:
+            # If we encounter a bool, we are going to need to treat it as
+            # binary explicitly because we are finally pedantic enough in the
+            # expression system to not allow some of the mixing we will need
+            # (like summing a LinearExpression with a bool)
+            return False, int(child)
+
         if child.__class__ in EXPR.native_types:
             return False, child
 

From 18e5f48ecfd420ea2bdf3571e8038941fc0ae43b Mon Sep 17 00:00:00 2001
From: Emma Johnson <emma58@users.noreply.github.com>
Date: Wed, 10 Jul 2024 13:52:03 -0600
Subject: [PATCH 36/69] black

---
 pyomo/contrib/cp/tests/test_logical_to_disjunctive.py | 4 +---
 1 file changed, 1 insertion(+), 3 deletions(-)

diff --git a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
index 9c79364c367..d940468900b 100755
--- a/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
+++ b/pyomo/contrib/cp/tests/test_logical_to_disjunctive.py
@@ -222,9 +222,7 @@ def test_equivalent_to_True(self):
         assertExpressionsEqual(self, m.cons[3].expr, m.z[2] + (1 - 1) >= 1)
 
         # z[3] == a v ! c
-        assertExpressionsEqual(
-            self, m.cons[4].expr, (1 - m.z[3]) + m.z[1] >= 1
-        )
+        assertExpressionsEqual(self, m.cons[4].expr, (1 - m.z[3]) + m.z[1] >= 1)
         assertExpressionsEqual(self, m.cons[5].expr, m.z[3] + (1 - m.z[1]) >= 1)
         assertExpressionsEqual(self, m.cons[6].expr, m.z[3] + 1 >= 1)
 

From 89be970cad522d9a35911f5b38c90ef0198f9cdb Mon Sep 17 00:00:00 2001
From: Emma Johnson <emma58@users.noreply.github.com>
Date: Tue, 23 Jul 2024 16:07:07 -0600
Subject: [PATCH 37/69] Checking bool inside of native types

---
 .../cp/transform/logical_to_disjunctive_walker.py   | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
index 47a85076ce1..4bce9c7d6af 100644
--- a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
+++ b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
@@ -243,14 +243,13 @@ def initializeWalker(self, expr):
         return True, expr
 
     def beforeChild(self, node, child, child_idx):
-        if child.__class__ is bool:
-            # If we encounter a bool, we are going to need to treat it as
-            # binary explicitly because we are finally pedantic enough in the
-            # expression system to not allow some of the mixing we will need
-            # (like summing a LinearExpression with a bool)
-            return False, int(child)
-
         if child.__class__ in EXPR.native_types:
+            if child.__class__ is bool:
+                # If we encounter a bool, we are going to need to treat it as
+                # binary explicitly because we are finally pedantic enough in the
+                # expression system to not allow some of the mixing we will need
+                # (like summing a LinearExpression with a bool)
+                return False, int(child)
             return False, child
 
         if child.is_numeric_type():

From 507d828acc24f5b3b8b5edfb466f112655d4c227 Mon Sep 17 00:00:00 2001
From: Emma Johnson <emma58@users.noreply.github.com>
Date: Tue, 23 Jul 2024 16:11:00 -0600
Subject: [PATCH 38/69] Removing the integrality check for Param because it
 isn't really future-proof and it is redundant with the same check in the
 expression nodes where it actually is required

---
 .../cp/transform/logical_to_disjunctive_walker.py        | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
index 4bce9c7d6af..95cbaf57fa5 100644
--- a/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
+++ b/pyomo/contrib/cp/transform/logical_to_disjunctive_walker.py
@@ -51,14 +51,7 @@ def _dispatch_var(visitor, node):
 
 
 def _dispatch_param(visitor, node):
-    if int(value(node)) == value(node):
-        return False, node
-    else:
-        raise ValueError(
-            "Found non-integer valued Param '%s' in a logical "
-            "expression. This cannot be written to a disjunctive "
-            "form." % node.name
-        )
+    return False, node
 
 
 def _dispatch_expression(visitor, node):

From 08ccf404447c59d532fca8150f5b9233e0652abd Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Wed, 24 Jul 2024 10:52:34 -0600
Subject: [PATCH 39/69] Rename normalize_constraint -> to_bounded_expression

---
 pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp |  2 +-
 pyomo/contrib/appsi/cmodel/src/lp_writer.cpp  |  2 +-
 pyomo/contrib/appsi/cmodel/src/nl_writer.cpp  |  2 +-
 pyomo/core/base/constraint.py                 | 45 +++++++++++++++----
 pyomo/core/kernel/constraint.py               |  2 +-
 pyomo/gdp/plugins/bilinear.py                 |  2 +-
 pyomo/gdp/plugins/cuttingplane.py             |  4 +-
 .../plugins/solvers/persistent_solver.py      |  2 +-
 8 files changed, 44 insertions(+), 17 deletions(-)

diff --git a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
index 708cfd9e073..ca865d429e2 100644
--- a/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/fbbt_model.cpp
@@ -205,7 +205,7 @@ void process_fbbt_constraints(FBBTModel *model, PyomoExprTypes &expr_types,
   py::handle con_body;
 
   for (py::handle c : cons) {
-    lower_body_upper = c.attr("normalize_constraint")();
+    lower_body_upper = c.attr("to_bounded_expression")();
     con_lb = lower_body_upper[0];
     con_body = lower_body_upper[1];
     con_ub = lower_body_upper[2];
diff --git a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
index 996bb34f564..f33060ee523 100644
--- a/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/lp_writer.cpp
@@ -289,7 +289,7 @@ void process_lp_constraints(py::list cons, py::object writer) {
   py::object nonlinear_expr;
   PyomoExprTypes expr_types = PyomoExprTypes();
   for (py::handle c : cons) {
-    lower_body_upper = c.attr("normalize_constraint")();
+    lower_body_upper = c.attr("to_bounded_expression")();
     cname = getSymbol(c, labeler);
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = true);
diff --git a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
index 477bdd87aee..854262496ea 100644
--- a/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
+++ b/pyomo/contrib/appsi/cmodel/src/nl_writer.cpp
@@ -527,7 +527,7 @@ void process_nl_constraints(NLWriter *nl_writer, PyomoExprTypes &expr_types,
   py::handle repn_nonlinear_expr;
 
   for (py::handle c : cons) {
-    lower_body_upper = c.attr("normalize_constraint")();
+    lower_body_upper = c.attr("to_bounded_expression")();
     repn = generate_standard_repn(
         lower_body_upper[1], "compute_values"_a = false, "quadratic"_a = false);
     _const = appsi_expr_from_pyomo_expr(repn.attr("constant"), var_map,
diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index defaea99dff..5a9d1da5af1 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -174,12 +174,35 @@ def __init__(self, expr=None, component=None):
 
     def __call__(self, exception=True):
         """Compute the value of the body of this constraint."""
-        body = self.normalize_constraint()[1]
+        body = self.to_bounded_expression()[1]
         if body.__class__ not in native_numeric_types:
             body = value(self.body, exception=exception)
         return body
 
-    def normalize_constraint(self):
+    def to_bounded_expression(self):
+        """Convert this constraint to a tuple of 3 expressions (lb, body, ub)
+
+        This method "standardizes" the expression into a 3-tuple of
+        expressions: (`lower_bound`, `body`, `upper_bound`).  Upon
+        conversion, `lower_bound` and `upper_bound` are guaranteed to be
+        `None`, numeric constants, or fixed (not necessarily constant)
+        expressions.
+
+        Note
+        ----
+        As this method operates on the *current state* of the
+        expression, the any required expression manipulations (and by
+        extension, the result) can change after fixing / unfixing
+        :py:class:`Var` objects.
+
+        Raises
+        ------
+
+        ValueError: Raised if the expression cannot be mapped to this
+            form (i.e., :py:class:`RangedExpression` constraints with
+            variable lower of upper bounds.
+
+        """
         expr = self._expr
         if expr.__class__ is RangedExpression:
             lb, body, ub = ans = expr.args
@@ -217,8 +240,12 @@ def normalize_constraint(self):
     def body(self):
         """Access the body of a constraint expression."""
         try:
-            ans = self.normalize_constraint()[1]
+            ans = self.to_bounded_expression()[1]
         except ValueError:
+            # It is possible that the expression is not currently valid
+            # (i.e., a ranged expression with a non-fixed bound).  We
+            # will catch that exception here and - if this actually *is*
+            # a RangedExpression - return the body.
             if self._expr.__class__ is RangedExpression:
                 _, ans, _ = self._expr.args
             else:
@@ -229,14 +256,14 @@ def body(self):
             #
             # [JDS 6/2024: it would be nice to remove this behavior,
             # although possibly unnecessary, as people should use
-            # normalize_constraint() instead]
+            # to_bounded_expression() instead]
             return as_numeric(ans)
         return ans
 
     @property
     def lower(self):
         """Access the lower bound of a constraint expression."""
-        ans = self.normalize_constraint()[0]
+        ans = self.to_bounded_expression()[0]
         if ans.__class__ in native_types and ans is not None:
             # Historically, constraint.lower was guaranteed to return a type
             # derived from Pyomo NumericValue (or None).  Replicate that
@@ -250,7 +277,7 @@ def lower(self):
     @property
     def upper(self):
         """Access the upper bound of a constraint expression."""
-        ans = self.normalize_constraint()[2]
+        ans = self.to_bounded_expression()[2]
         if ans.__class__ in native_types and ans is not None:
             # Historically, constraint.upper was guaranteed to return a type
             # derived from Pyomo NumericValue (or None).  Replicate that
@@ -264,7 +291,7 @@ def upper(self):
     @property
     def lb(self):
         """Access the value of the lower bound of a constraint expression."""
-        bound = self.normalize_constraint()[0]
+        bound = self.to_bounded_expression()[0]
         if bound is None:
             return None
         if bound.__class__ not in native_numeric_types:
@@ -282,7 +309,7 @@ def lb(self):
     @property
     def ub(self):
         """Access the value of the upper bound of a constraint expression."""
-        bound = self.normalize_constraint()[2]
+        bound = self.to_bounded_expression()[2]
         if bound is None:
             return None
         if bound.__class__ not in native_numeric_types:
@@ -824,7 +851,7 @@ class SimpleConstraint(metaclass=RenamedClass):
     {
         'add',
         'set_value',
-        'normalize_constraint',
+        'to_bounded_expression',
         'body',
         'lower',
         'upper',
diff --git a/pyomo/core/kernel/constraint.py b/pyomo/core/kernel/constraint.py
index 6b8c4c619f5..fe8eb8b2c1f 100644
--- a/pyomo/core/kernel/constraint.py
+++ b/pyomo/core/kernel/constraint.py
@@ -177,7 +177,7 @@ class _MutableBoundsConstraintMixin(object):
     # Define some of the IConstraint abstract methods
     #
 
-    def normalize_constraint(self):
+    def to_bounded_expression(self):
         return self.lower, self.body, self.upper
 
     @property
diff --git a/pyomo/gdp/plugins/bilinear.py b/pyomo/gdp/plugins/bilinear.py
index 70b6e83b52f..bc91836ea9c 100644
--- a/pyomo/gdp/plugins/bilinear.py
+++ b/pyomo/gdp/plugins/bilinear.py
@@ -77,7 +77,7 @@ def _transformBlock(self, block, instance):
         for component in block.component_data_objects(
             Constraint, active=True, descend_into=False
         ):
-            lb, body, ub = component.normalize_constraint()
+            lb, body, ub = component.to_bounded_expression()
             expr = self._transformExpression(body, instance)
             instance.bilinear_data_.c_body[id(component)] = body
             component.set_value((lb, expr, ub))
diff --git a/pyomo/gdp/plugins/cuttingplane.py b/pyomo/gdp/plugins/cuttingplane.py
index a757f23c826..4cef098eba9 100644
--- a/pyomo/gdp/plugins/cuttingplane.py
+++ b/pyomo/gdp/plugins/cuttingplane.py
@@ -400,7 +400,7 @@ def back_off_constraint_with_calculated_cut_violation(
     val = value(transBlock_rHull.infeasibility_objective) - TOL
     if val <= 0:
         logger.info("\tBacking off cut by %s" % val)
-        lb, body, ub = cut.normalize_constraint()
+        lb, body, ub = cut.to_bounded_expression()
         cut.set_value((lb, body + abs(val), ub))
     # else there is nothing to do: restore the objective
     transBlock_rHull.del_component(transBlock_rHull.infeasibility_objective)
@@ -425,7 +425,7 @@ def back_off_constraint_by_fixed_tolerance(
                    this callback
     TOL: An absolute tolerance to be added to make cut more conservative.
     """
-    lb, body, ub = cut.normalize_constraint()
+    lb, body, ub = cut.to_bounded_expression()
     cut.set_value((lb, body + TOL, ub))
 
 
diff --git a/pyomo/solvers/plugins/solvers/persistent_solver.py b/pyomo/solvers/plugins/solvers/persistent_solver.py
index ef96bfa339f..ef883fe5496 100644
--- a/pyomo/solvers/plugins/solvers/persistent_solver.py
+++ b/pyomo/solvers/plugins/solvers/persistent_solver.py
@@ -262,7 +262,7 @@ def _add_and_collect_column_data(self, var, obj_coef, constraints, coefficients)
         coeff_list = list()
         constr_list = list()
         for val, c in zip(coefficients, constraints):
-            lb, body, ub = c.normalize_constraint()
+            lb, body, ub = c.to_bounded_expression()
             body += val * var
             c.set_value((lb, body, ub))
             self._vars_referenced_by_con[c].add(var)

From 2529557649b39abd0e963632f0cb608a2a2a6627 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Wed, 24 Jul 2024 11:23:51 -0600
Subject: [PATCH 40/69] NFC: fix comment typo

---
 pyomo/contrib/fbbt/tests/test_fbbt.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/fbbt/tests/test_fbbt.py b/pyomo/contrib/fbbt/tests/test_fbbt.py
index ff1cc8a5cfb..83e69233bb5 100644
--- a/pyomo/contrib/fbbt/tests/test_fbbt.py
+++ b/pyomo/contrib/fbbt/tests/test_fbbt.py
@@ -1339,7 +1339,7 @@ def setUp(self) -> None:
     def test_ranged_expression(self):
         # The python version of FBBT is slightly more flexible than
         # APPSI's cmodel (it allows - and correctly handles -
-        # RangedExpressions with variable lower / upper bounds.  If we
+        # RangedExpressions with variable lower / upper bounds).  If we
         # ever port that functionality into APPSI, then this test can be
         # moved into the base class.
         m = pyo.ConcreteModel()

From ac479032cb52dc14708bcd8d499642b874e76332 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 05:36:58 -0600
Subject: [PATCH 41/69] Rework _initialize and avoid redundant call to iter()

---
 pyomo/core/base/set.py            | 36 ++++++++++++++-----------------
 pyomo/core/tests/unit/test_set.py |  9 ++++++++
 2 files changed, 25 insertions(+), 20 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index 87e9549ec8d..e8ac978c355 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1386,29 +1386,27 @@ def set_value(self, val):
         self.update(val)
 
     def _initialize(self, val):
-        self.update(val)
+        try:
+            self.update(val)
+        except TypeError as e:
+            if 'not iterable' in str(e):
+                logger.error(
+                    "Initializer for Set %s returned non-iterable object "
+                    "of type %s."
+                    % (
+                        self.name,
+                        (val if val.__class__ is type else type(val).__name__),
+                    )
+                )
+            raise
 
     def update(self, values):
-        # _values was initialized above...
-        #
         # Special case: set operations that are not first attached
         # to the model must be constructed.
         if isinstance(values, SetOperator):
             values.construct()
-        try:
-            val_iter = iter(values)
-        except TypeError:
-            logger.error(
-                "Initializer for Set %s%s returned non-iterable object "
-                "of type %s."
-                % (
-                    self.name,
-                    ("[%s]" % (index,) if self.is_indexed() else ""),
-                    (values if values.__class__ is type else type(values).__name__),
-                )
-            )
-            raise
-
+        # It is important that val_iter is an actual iterator
+        val_iter = iter(values)
         if self._dimen is not None:
             if normalize_index.flatten:
                 val_iter = self._cb_normalized_dimen_verifier(self._dimen, val_iter)
@@ -1472,8 +1470,6 @@ def _cb_validate(self, validate, block, val_iter):
             yield value
 
     def _cb_normalized_dimen_verifier(self, dimen, val_iter):
-        # It is important that the iterator is an actual iterator
-        val_iter = iter(val_iter)
         for value in val_iter:
             if value.__class__ is tuple:
                 if dimen == len(value):
@@ -1833,7 +1829,7 @@ def _initialize(self, val):
                 "This WILL potentially lead to nondeterministic behavior "
                 "in Pyomo" % (self.name, type(val).__name__)
             )
-        super().update(val)
+        super()._initialize(val)
 
     def set_value(self, val):
         if type(val) in Set._UnorderedInitializers:
diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index ff0aaa5600b..7fe7d1fd7ae 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -3811,6 +3811,7 @@ def I_init(m):
         self.assertEqual(m.I.data(), (4, 3, 2, 1))
         self.assertEqual(m.I.dimen, 1)
 
+    def test_initialize_with_noniterable(self):
         output = StringIO()
         with LoggingIntercept(output, 'pyomo.core'):
             with self.assertRaisesRegex(TypeError, "'int' object is not iterable"):
@@ -3819,6 +3820,14 @@ def I_init(m):
             ref = "Initializer for Set I returned non-iterable object of type int."
             self.assertIn(ref, output.getvalue())
 
+        output = StringIO()
+        with LoggingIntercept(output, 'pyomo.core'):
+            with self.assertRaisesRegex(TypeError, "'int' object is not iterable"):
+                m = ConcreteModel()
+                m.I = Set([1,2], initialize=5)
+            ref = "Initializer for Set I[1] returned non-iterable object of type int."
+            self.assertIn(ref, output.getvalue())
+
     def test_scalar_indexed_api(self):
         m = ConcreteModel()
         m.I = Set(initialize=range(3))

From 9aab8c2e487cefd8f0ed05ef19351cf19ed6e88b Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 05:37:35 -0600
Subject: [PATCH 42/69] Duplicate _update_impl workaround from ordered sets in
 sorted sets

---
 pyomo/core/base/set.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index e8ac978c355..b62f61f8557 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1897,8 +1897,12 @@ def __reversed__(self):
 
     def _update_impl(self, values):
         for val in values:
+            # Note that we reset _ordered_values within the loop because
+            # of an old example where the initializer rule makes
+            # reference to values previously inserted into the Set
+            # (which triggered the creation of the _ordered_values)
+            self._ordered_values = None
             self._values[val] = None
-        self._ordered_values = None
 
     # Note: removing data does not affect the sorted flag
     # def remove(self, val):

From 26488f38c49961d5b8c6ac5b847096c720971edb Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 07:26:58 -0600
Subject: [PATCH 43/69] Catch (and test) an edge case when flattening sets

---
 pyomo/core/base/set.py            | 19 +++++++++++--------
 pyomo/core/tests/unit/test_set.py | 15 +++++++++++++++
 2 files changed, 26 insertions(+), 8 deletions(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index b62f61f8557..045fc134f74 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1471,16 +1471,19 @@ def _cb_validate(self, validate, block, val_iter):
 
     def _cb_normalized_dimen_verifier(self, dimen, val_iter):
         for value in val_iter:
-            if value.__class__ is tuple:
-                if dimen == len(value):
-                    yield value[0] if dimen == 1 else value
+            if value.__class__ in native_types:
+                if dimen == 1:
+                    yield value
                     continue
-            elif dimen == 1 and value.__class__ in native_types:
-                yield value
-                continue
+                normalized_value = value
+            else:
+                normalized_value = normalize_index(value)
+                # Note: normalize_index() will never return a 1-tuple
+                if normalized_value.__class__ is tuple:
+                    if dimen == len(normalized_value):
+                        yield normalized_value[0] if dimen == 1 else normalized_value
+                        continue
 
-            # Note: normalize_index() will never return a 1-tuple
-            normalized_value = normalize_index(value)
             _d = len(normalized_value) if normalized_value.__class__ is tuple else 1
             if _d == dimen:
                 yield normalized_value
diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index 7fe7d1fd7ae..2a2651ca554 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -5256,6 +5256,21 @@ def Bindex(m):
         self.assertIs(m.K.index_set()._domain, Integers)
         self.assertEqual(m.K.index_set(), [0, 1, 2, 3, 4])
 
+    def test_normalize_index(self):
+        try:
+            _oldFlatten = normalize_index.flatten
+            normalize_index.flatten = True
+
+            m = ConcreteModel()
+            with self.assertRaisesRegex(
+                ValueError,
+                r"The value=\(\(2, 3\),\) has dimension 2 and is not "
+                "valid for Set I which has dimen=1",
+            ):
+                m.I = Set(initialize=[1, ((2, 3),)])
+        finally:
+            normalize_index.flatten = _oldFlatten
+
     def test_no_normalize_index(self):
         try:
             _oldFlatten = normalize_index.flatten

From 0fc78ad21925e647da1f58483e2ab37be48f716a Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 07:27:25 -0600
Subject: [PATCH 44/69] NFC: apply black

---
 pyomo/core/tests/unit/test_set.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/core/tests/unit/test_set.py b/pyomo/core/tests/unit/test_set.py
index 2a2651ca554..8c0360d618b 100644
--- a/pyomo/core/tests/unit/test_set.py
+++ b/pyomo/core/tests/unit/test_set.py
@@ -3824,7 +3824,7 @@ def test_initialize_with_noniterable(self):
         with LoggingIntercept(output, 'pyomo.core'):
             with self.assertRaisesRegex(TypeError, "'int' object is not iterable"):
                 m = ConcreteModel()
-                m.I = Set([1,2], initialize=5)
+                m.I = Set([1, 2], initialize=5)
             ref = "Initializer for Set I[1] returned non-iterable object of type int."
             self.assertIn(ref, output.getvalue())
 

From 49a53d89cef82445091afb352a59129914b73f0e Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 07:46:49 -0600
Subject: [PATCH 45/69] NFC: fix doc typo

Co-authored-by: Bethany Nicholson <blnicho@users.noreply.github.com>
---
 pyomo/core/base/constraint.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index 5a9d1da5af1..b79bc178e80 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -191,7 +191,7 @@ def to_bounded_expression(self):
         Note
         ----
         As this method operates on the *current state* of the
-        expression, the any required expression manipulations (and by
+        expression, any required expression manipulations (and by
         extension, the result) can change after fixing / unfixing
         :py:class:`Var` objects.
 

From 21502959923070e54076d943647fcaa723a224cc Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 07:47:04 -0600
Subject: [PATCH 46/69] NFC: fix doc typo

Co-authored-by: Bethany Nicholson <blnicho@users.noreply.github.com>
---
 pyomo/core/base/constraint.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/core/base/constraint.py b/pyomo/core/base/constraint.py
index b79bc178e80..bc9a32f5404 100644
--- a/pyomo/core/base/constraint.py
+++ b/pyomo/core/base/constraint.py
@@ -200,7 +200,7 @@ def to_bounded_expression(self):
 
         ValueError: Raised if the expression cannot be mapped to this
             form (i.e., :py:class:`RangedExpression` constraints with
-            variable lower of upper bounds.
+            variable lower or upper bounds.
 
         """
         expr = self._expr

From d8900b296b698b29d5b285dd115cd899014994c0 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 08:10:54 -0600
Subject: [PATCH 47/69] Remove redundant header, simplify imports

---
 pyomo/contrib/sensitivity_toolbox/sens.py | 13 +------------
 1 file changed, 1 insertion(+), 12 deletions(-)

diff --git a/pyomo/contrib/sensitivity_toolbox/sens.py b/pyomo/contrib/sensitivity_toolbox/sens.py
index a3d69b2c7b1..818f13cb789 100644
--- a/pyomo/contrib/sensitivity_toolbox/sens.py
+++ b/pyomo/contrib/sensitivity_toolbox/sens.py
@@ -9,16 +9,6 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
-# ______________________________________________________________________________
-#
-# Pyomo: Python Optimization Modeling Objects
-# Copyright (c) 2008-2024
-#  National Technology and Engineering Solutions of Sandia, LLC
-# Under the terms of Contract DE-NA0003525 with National Technology and
-# Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
-# rights in this software.
-# This software is distributed under the 3-clause BSD License
-# ______________________________________________________________________________
 from pyomo.environ import (
     Param,
     Var,
@@ -36,6 +26,7 @@
 from pyomo.core.expr import ExpressionReplacementVisitor
 
 from pyomo.common.modeling import unique_component_name
+from pyomo.common.dependencies import numpy as np, scipy
 from pyomo.common.deprecation import deprecated
 from pyomo.common.tempfiles import TempfileManager
 from pyomo.opt import SolverFactory, SolverStatus
@@ -44,8 +35,6 @@
 import os
 import io
 import shutil
-from pyomo.common.dependencies import numpy as np, numpy_available
-from pyomo.common.dependencies import scipy, scipy_available
 
 logger = logging.getLogger('pyomo.contrib.sensitivity_toolbox')
 

From bb17e9c9bc115238c34ff67fd9e3354eff63c782 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 08:11:41 -0600
Subject: [PATCH 48/69] Update constraint processing to leverage new Constraint
 internal storage

---
 pyomo/contrib/sensitivity_toolbox/sens.py | 35 +++++++++++++----------
 1 file changed, 20 insertions(+), 15 deletions(-)

diff --git a/pyomo/contrib/sensitivity_toolbox/sens.py b/pyomo/contrib/sensitivity_toolbox/sens.py
index 818f13cb789..34fbb92327a 100644
--- a/pyomo/contrib/sensitivity_toolbox/sens.py
+++ b/pyomo/contrib/sensitivity_toolbox/sens.py
@@ -24,6 +24,7 @@
 
 from pyomo.common.sorting import sorted_robust
 from pyomo.core.expr import ExpressionReplacementVisitor
+from pyomo.core.expr.numvalue import is_potentially_variable
 
 from pyomo.common.modeling import unique_component_name
 from pyomo.common.dependencies import numpy as np, scipy
@@ -673,25 +674,29 @@ def _replace_parameters_in_constraints(self, variableSubMap):
         )
         last_idx = 0
         for con in old_con_list:
-            if con.equality or con.lower is None or con.upper is None:
-                new_expr = param_replacer.walk_expression(con.expr)
-                block.constList.add(expr=new_expr)
+            new_expr = param_replacer.walk_expression(con.expr)
+            # TODO: We could only create new constraints for expressions
+            # where substitution actually happened, but that breaks some
+            # current tests:
+            #
+            # if new_expr is con.expr:
+            #     # No params were substituted.  We can ignore this constraint
+            #     continue
+            if new_expr.nargs() == 3 and (
+                is_potentially_variable(new_expr.arg(0))
+                or is_potentially_variable(new_expr.arg(2))
+            ):
+                # This is a potentially "invalid" range constraint: it
+                # may now have variables in the bounds.  For safety, we
+                # will split it into two simple inequalities.
+                block.constList.add(expr=(new_expr.arg(0) <= new_expr.arg(1)))
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
-            else:
-                # Constraint must be a ranged inequality, break into
-                # separate constraints
-                new_body = param_replacer.walk_expression(con.body)
-                new_lower = param_replacer.walk_expression(con.lower)
-                new_upper = param_replacer.walk_expression(con.upper)
-
-                # Add constraint for lower bound
-                block.constList.add(expr=(new_lower <= new_body))
+                block.constList.add(expr=(new_expr.arg(1) <= new_expr.arg(2)))
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
-
-                # Add constraint for upper bound
-                block.constList.add(expr=(new_body <= new_upper))
+            else:
+                block.constList.add(expr=new_expr)
                 last_idx += 1
                 new_old_comp_map[block.constList[last_idx]] = con
             con.deactivate()

From 7df3801d9679d88e4f195fd38bc6f0a6750a08c2 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Mon, 29 Jul 2024 08:28:11 -0600
Subject: [PATCH 49/69] Avoid duplicate logging of unordered Set data warning

---
 pyomo/core/base/set.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/pyomo/core/base/set.py b/pyomo/core/base/set.py
index fadd0c3e7d0..e4a6d13e96e 100644
--- a/pyomo/core/base/set.py
+++ b/pyomo/core/base/set.py
@@ -1387,7 +1387,10 @@ def set_value(self, val):
 
     def _initialize(self, val):
         try:
-            self.update(val)
+            # We want to explicitly call the update() on *this class* to
+            # bypass potential double logging of the use of unordered
+            # data with ordered Sets
+            FiniteSetData.update(self, val)
         except TypeError as e:
             if 'not iterable' in str(e):
                 logger.error(

From 017258dd593817f96ee9da4207e379807e5c8e52 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 12:33:02 -0600
Subject: [PATCH 50/69] Restore 'options' as an alias of the new
 config.solver_options

---
 pyomo/contrib/solver/base.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/solver/base.py b/pyomo/contrib/solver/base.py
index 98bf3836004..ac07bda4e67 100644
--- a/pyomo/contrib/solver/base.py
+++ b/pyomo/contrib/solver/base.py
@@ -353,15 +353,20 @@ def __init__(self, **kwargs):
             raise NotImplementedError('Still working on this')
         # There is no reason for a user to be trying to mix both old
         # and new options. That is silly. So we will yell at them.
-        self.options = kwargs.pop('options', None)
+        _options = kwargs.pop('options', None)
         if 'solver_options' in kwargs:
-            if self.options is not None:
+            if _options is not None:
                 raise ValueError(
                     "Both 'options' and 'solver_options' were requested. "
                     "Please use one or the other, not both."
                 )
-            self.options = kwargs.pop('solver_options')
+            _options = kwargs.pop('solver_options')
+        if _options is not None:
+            kwargs['solver_options'] = _options
         super().__init__(**kwargs)
+        # Make the legacy 'options' attribute an alias of the new
+        # config.solver_options
+        self.options = self.config.solver_options
 
     #
     # Support "with" statements

From 61bb290ee5a2b1c9600d9abc7b381202b52dc238 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 13:40:20 -0600
Subject: [PATCH 51/69] Update tests to not directly instantiate the
 LegacySolverWrapper mixin

---
 pyomo/contrib/solver/tests/unit/test_base.py | 35 ++++++++++++--------
 1 file changed, 21 insertions(+), 14 deletions(-)

diff --git a/pyomo/contrib/solver/tests/unit/test_base.py b/pyomo/contrib/solver/tests/unit/test_base.py
index b52f96ba903..b7937d16af3 100644
--- a/pyomo/contrib/solver/tests/unit/test_base.py
+++ b/pyomo/contrib/solver/tests/unit/test_base.py
@@ -16,6 +16,10 @@
 from pyomo.contrib.solver import base
 
 
+class _LegacyWrappedSolverBase(base.LegacySolverWrapper, base.SolverBase):
+    pass
+
+
 class TestSolverBase(unittest.TestCase):
     def test_abstract_member_list(self):
         expected_list = ['solve', 'available', 'version']
@@ -192,11 +196,13 @@ def test_class_method_list(self):
         ]
         self.assertEqual(sorted(expected_list), sorted(method_list))
 
+    @unittest.mock.patch.multiple(_LegacyWrappedSolverBase, __abstractmethods__=set())
     def test_context_manager(self):
-        with base.LegacySolverWrapper() as instance:
-            with self.assertRaises(AttributeError):
-                instance.available()
+        with _LegacyWrappedSolverBase() as instance:
+            self.assertIsInstance(instance, _LegacyWrappedSolverBase)
+            self.assertFalse(instance.available(False))
 
+    @unittest.mock.patch.multiple(_LegacyWrappedSolverBase, __abstractmethods__=set())
     def test_map_config(self):
         # Create a fake/empty config structure that can be added to an empty
         # instance of LegacySolverWrapper
@@ -205,7 +211,7 @@ def test_map_config(self):
             'solver_options',
             ConfigDict(implicit=True, description="Options to pass to the solver."),
         )
-        instance = base.LegacySolverWrapper()
+        instance = _LegacyWrappedSolverBase()
         instance.config = self.config
         instance._map_config(
             True, False, False, 20, True, False, None, None, None, False, None, None
@@ -272,20 +278,21 @@ def test_map_config(self):
         with self.assertRaises(AttributeError):
             print(instance.config.keepfiles)
 
+    @unittest.mock.patch.multiple(_LegacyWrappedSolverBase, __abstractmethods__=set())
     def test_solver_options_behavior(self):
         # options can work in multiple ways (set from instantiation, set
         # after instantiation, set during solve).
         # Test case 1: Set at instantiation
-        solver = base.LegacySolverWrapper(options={'max_iter': 6})
+        solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
         self.assertEqual(solver.options, {'max_iter': 6})
 
         # Test case 2: Set later
-        solver = base.LegacySolverWrapper()
+        solver = _LegacyWrappedSolverBase()
         solver.options = {'max_iter': 4, 'foo': 'bar'}
         self.assertEqual(solver.options, {'max_iter': 4, 'foo': 'bar'})
 
         # Test case 3: pass some options to the mapping (aka, 'solve' command)
-        solver = base.LegacySolverWrapper()
+        solver = _LegacyWrappedSolverBase()
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
@@ -296,7 +303,7 @@ def test_solver_options_behavior(self):
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
 
         # Test case 4: Set at instantiation and override during 'solve' call
-        solver = base.LegacySolverWrapper(options={'max_iter': 6})
+        solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
@@ -309,11 +316,11 @@ def test_solver_options_behavior(self):
 
         # solver_options are also supported
         # Test case 1: set at instantiation
-        solver = base.LegacySolverWrapper(solver_options={'max_iter': 6})
+        solver = _LegacyWrappedSolverBase(solver_options={'max_iter': 6})
         self.assertEqual(solver.options, {'max_iter': 6})
 
         # Test case 2: pass some solver_options to the mapping (aka, 'solve' command)
-        solver = base.LegacySolverWrapper()
+        solver = _LegacyWrappedSolverBase()
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
@@ -324,7 +331,7 @@ def test_solver_options_behavior(self):
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
 
         # Test case 3: Set at instantiation and override during 'solve' call
-        solver = base.LegacySolverWrapper(solver_options={'max_iter': 6})
+        solver = _LegacyWrappedSolverBase(solver_options={'max_iter': 6})
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
@@ -337,7 +344,7 @@ def test_solver_options_behavior(self):
 
         # users can mix... sort of
         # Test case 1: Initialize with options, solve with solver_options
-        solver = base.LegacySolverWrapper(options={'max_iter': 6})
+        solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
@@ -351,11 +358,11 @@ def test_solver_options_behavior(self):
         # do we know what to do with it then?
         # Test case 1: Class instance
         with self.assertRaises(ValueError):
-            solver = base.LegacySolverWrapper(
+            solver = _LegacyWrappedSolverBase(
                 options={'max_iter': 6}, solver_options={'max_iter': 4}
             )
         # Test case 2: Passing to `solve`
-        solver = base.LegacySolverWrapper()
+        solver = _LegacyWrappedSolverBase()
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',

From 896eae16203d9b3ab8865ce68ae5f13db6057b23 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 14:07:31 -0600
Subject: [PATCH 52/69] LegacySolverWrapper: 'options' and 'config' should be
 singleton attributes

---
 pyomo/contrib/solver/base.py | 11 ++++++++---
 1 file changed, 8 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/solver/base.py b/pyomo/contrib/solver/base.py
index ac07bda4e67..4fe8bee4e53 100644
--- a/pyomo/contrib/solver/base.py
+++ b/pyomo/contrib/solver/base.py
@@ -377,6 +377,14 @@ def __enter__(self):
     def __exit__(self, t, v, traceback):
         """Exit statement - enables `with` statements."""
 
+    def __setattr__(self, attr, value):
+        # 'options' and 'config' are really singleton attributes.  Map
+        # any assignment to set_value()
+        if attr in ('options', 'config') and attr in self.__dict__:
+            getattr(self, attr).set_value(value)
+        else:
+            super().__setattr__(attr, value)
+
     def _map_config(
         self,
         tee=NOTSET,
@@ -395,7 +403,6 @@ def _map_config(
         writer_config=NOTSET,
     ):
         """Map between legacy and new interface configuration options"""
-        self.config = self.config()
         if 'report_timing' not in self.config:
             self.config.declare(
                 'report_timing', ConfigValue(domain=bool, default=False)
@@ -410,8 +417,6 @@ def _map_config(
             self.config.time_limit = timelimit
         if report_timing is not NOTSET:
             self.config.report_timing = report_timing
-        if self.options is not None:
-            self.config.solver_options.set_value(self.options)
         if (options is not NOTSET) and (solver_options is not NOTSET):
             # There is no reason for a user to be trying to mix both old
             # and new options. That is silly. So we will yell at them.

From 6975cb6b068beda72f748396dd51644f03b895ce Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 14:08:02 -0600
Subject: [PATCH 53/69] Update 'options' tests to reflect new bas class

---
 pyomo/contrib/solver/tests/unit/test_base.py | 50 +++++++-------------
 1 file changed, 16 insertions(+), 34 deletions(-)

diff --git a/pyomo/contrib/solver/tests/unit/test_base.py b/pyomo/contrib/solver/tests/unit/test_base.py
index b7937d16af3..ba62f97542a 100644
--- a/pyomo/contrib/solver/tests/unit/test_base.py
+++ b/pyomo/contrib/solver/tests/unit/test_base.py
@@ -285,73 +285,49 @@ def test_solver_options_behavior(self):
         # Test case 1: Set at instantiation
         solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
         self.assertEqual(solver.options, {'max_iter': 6})
+        self.assertEqual(solver.config.solver_options, {'max_iter': 6})
 
         # Test case 2: Set later
         solver = _LegacyWrappedSolverBase()
         solver.options = {'max_iter': 4, 'foo': 'bar'}
         self.assertEqual(solver.options, {'max_iter': 4, 'foo': 'bar'})
+        self.assertEqual(solver.config.solver_options, {'max_iter': 4, 'foo': 'bar'})
 
         # Test case 3: pass some options to the mapping (aka, 'solve' command)
         solver = _LegacyWrappedSolverBase()
-        config = ConfigDict(implicit=True)
-        config.declare(
-            'solver_options',
-            ConfigDict(implicit=True, description="Options to pass to the solver."),
-        )
-        solver.config = config
         solver._map_config(options={'max_iter': 4})
+        self.assertEqual(solver.options, {'max_iter': 4})
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
 
         # Test case 4: Set at instantiation and override during 'solve' call
         solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
-        config = ConfigDict(implicit=True)
-        config.declare(
-            'solver_options',
-            ConfigDict(implicit=True, description="Options to pass to the solver."),
-        )
-        solver.config = config
         solver._map_config(options={'max_iter': 4})
+        self.assertEqual(solver.options, {'max_iter': 4})
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
-        self.assertEqual(solver.options, {'max_iter': 6})
 
         # solver_options are also supported
         # Test case 1: set at instantiation
         solver = _LegacyWrappedSolverBase(solver_options={'max_iter': 6})
         self.assertEqual(solver.options, {'max_iter': 6})
+        self.assertEqual(solver.config.solver_options, {'max_iter': 6})
 
         # Test case 2: pass some solver_options to the mapping (aka, 'solve' command)
         solver = _LegacyWrappedSolverBase()
-        config = ConfigDict(implicit=True)
-        config.declare(
-            'solver_options',
-            ConfigDict(implicit=True, description="Options to pass to the solver."),
-        )
-        solver.config = config
         solver._map_config(solver_options={'max_iter': 4})
+        self.assertEqual(solver.options, {'max_iter': 4})
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
 
         # Test case 3: Set at instantiation and override during 'solve' call
         solver = _LegacyWrappedSolverBase(solver_options={'max_iter': 6})
-        config = ConfigDict(implicit=True)
-        config.declare(
-            'solver_options',
-            ConfigDict(implicit=True, description="Options to pass to the solver."),
-        )
-        solver.config = config
         solver._map_config(solver_options={'max_iter': 4})
+        self.assertEqual(solver.options, {'max_iter': 4})
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
-        self.assertEqual(solver.options, {'max_iter': 6})
 
         # users can mix... sort of
         # Test case 1: Initialize with options, solve with solver_options
         solver = _LegacyWrappedSolverBase(options={'max_iter': 6})
-        config = ConfigDict(implicit=True)
-        config.declare(
-            'solver_options',
-            ConfigDict(implicit=True, description="Options to pass to the solver."),
-        )
-        solver.config = config
         solver._map_config(solver_options={'max_iter': 4})
+        self.assertEqual(solver.options, {'max_iter': 4})
         self.assertEqual(solver.config.solver_options, {'max_iter': 4})
 
         # users CANNOT initialize both values at the same time, because how
@@ -363,14 +339,20 @@ def test_solver_options_behavior(self):
             )
         # Test case 2: Passing to `solve`
         solver = _LegacyWrappedSolverBase()
+        with self.assertRaises(ValueError):
+            solver._map_config(solver_options={'max_iter': 4}, options={'max_iter': 6})
+
+        # Test that assignment to maps to set_vlaue:
+        solver = _LegacyWrappedSolverBase()
         config = ConfigDict(implicit=True)
         config.declare(
             'solver_options',
             ConfigDict(implicit=True, description="Options to pass to the solver."),
         )
         solver.config = config
-        with self.assertRaises(ValueError):
-            solver._map_config(solver_options={'max_iter': 4}, options={'max_iter': 6})
+        solver.config.solver_options.max_iter = 6
+        self.assertEqual(solver.options, {'max_iter': 6})
+        self.assertEqual(solver.config.solver_options, {'max_iter': 6})
 
     def test_map_results(self):
         # Unclear how to test this

From a61df0c68d98a60090bc884521769ced4caa5b46 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 14:12:52 -0600
Subject: [PATCH 54/69] Fix typo

---
 pyomo/contrib/solver/tests/unit/test_base.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/solver/tests/unit/test_base.py b/pyomo/contrib/solver/tests/unit/test_base.py
index ba62f97542a..e9ea717593f 100644
--- a/pyomo/contrib/solver/tests/unit/test_base.py
+++ b/pyomo/contrib/solver/tests/unit/test_base.py
@@ -342,7 +342,7 @@ def test_solver_options_behavior(self):
         with self.assertRaises(ValueError):
             solver._map_config(solver_options={'max_iter': 4}, options={'max_iter': 6})
 
-        # Test that assignment to maps to set_vlaue:
+        # Test that assignment to maps to set_value:
         solver = _LegacyWrappedSolverBase()
         config = ConfigDict(implicit=True)
         config.declare(

From 73ca0fcab8bdac103b8445a7d102969e63030259 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:09:15 -0600
Subject: [PATCH 55/69] Guard additional tests for pynumero availability

---
 pyomo/contrib/parmest/tests/test_examples.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index 450863a08a4..552c568cc7c 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -44,6 +44,7 @@ def test_model_with_constraint(self):
 
         rooney_biegler_with_constraint.main()
 
+    @unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
     @unittest.skipUnless(seaborn_available, "test requires seaborn")
     def test_parameter_estimation_example(self):
         from pyomo.contrib.parmest.examples.rooney_biegler import (
@@ -67,11 +68,12 @@ def test_likelihood_ratio_example(self):
         likelihood_ratio_example.main()
 
 
-@unittest.skipIf(
-    not parmest.parmest_available,
+@unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
+@unittest.skipUnless(ipopt_available, "The 'ipopt' solver is not available")
+@unittest.skipUnless(
+    parmest.parmest_available,
     "Cannot test parmest: required dependencies are missing",
 )
-@unittest.skipIf(not ipopt_available, "The 'ipopt' solver is not available")
 class TestReactionKineticsExamples(unittest.TestCase):
     @classmethod
     def setUpClass(self):
@@ -141,6 +143,7 @@ def test_model(self):
 
         reactor_design.main()
 
+    @unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
     def test_parameter_estimation_example(self):
         from pyomo.contrib.parmest.examples.reactor_design import (
             parameter_estimation_example,

From 31cb42014aeee6834da21945294165a3e529a612 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:09:53 -0600
Subject: [PATCH 56/69] Remove tests of k_aug functionality that is no longer
 supported

---
 pyomo/contrib/parmest/tests/test_parmest.py | 13 -------------
 1 file changed, 13 deletions(-)

diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index 1ff42a38d9e..fbe4290dd10 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -35,7 +35,6 @@
 
 is_osx = platform.mac_ver()[0] != ""
 ipopt_available = SolverFactory("ipopt").available()
-k_aug_available = SolverFactory("k_aug").available(exception_flag=False)
 pynumero_ASL_available = AmplInterface.available()
 testdir = this_file_dir()
 
@@ -199,12 +198,6 @@ def test_parallel_parmest(self):
             retcode = subprocess.call(rlist)
         self.assertEqual(retcode, 0)
 
-    @unittest.skipUnless(k_aug_available, "k_aug solver not found")
-    def test_theta_k_aug_for_Hessian(self):
-        # this will fail if k_aug is not installed
-        objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug")
-        self.assertAlmostEqual(objval, 4.4675, places=2)
-
     @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6)
@@ -1205,12 +1198,6 @@ def test_parallel_parmest(self):
             retcode = subprocess.call(rlist)
         assert retcode == 0
 
-    @unittest.skipUnless(k_aug_available, "k_aug solver not found")
-    def test_theta_k_aug_for_Hessian(self):
-        # this will fail if k_aug is not installed
-        objval, thetavals, Hessian = self.pest.theta_est(solver="k_aug")
-        self.assertAlmostEqual(objval, 4.4675, places=2)
-
     @unittest.skipIf(not pynumero_ASL_available, "pynumero_ASL is not available")
     def test_theta_est_cov(self):
         objval, thetavals, cov = self.pest.theta_est(calc_cov=True, cov_n=6)

From 3ba93e555fbcb1ea3221cc028c0e1642bb488e3c Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:11:17 -0600
Subject: [PATCH 57/69] Fix GJH version() to return the expected tuple

---
 pyomo/solvers/plugins/solvers/ASL.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/solvers/plugins/solvers/ASL.py b/pyomo/solvers/plugins/solvers/ASL.py
index 7acd59936b1..c912f2a30ee 100644
--- a/pyomo/solvers/plugins/solvers/ASL.py
+++ b/pyomo/solvers/plugins/solvers/ASL.py
@@ -108,7 +108,7 @@ def _get_version(self):
             if ver is None:
                 # Some ASL solvers do not export a version number
                 if results.stdout.strip().split()[-1].startswith('ASL('):
-                    return '0.0.0'
+                    return (0, 0, 0)
             return ver
         except OSError:
             pass

From d1c7952fa40cc340bb54bd87471e59a6d68a80ea Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:11:30 -0600
Subject: [PATCH 58/69] Update "pyomo help -s"

- Silence all output generated when gathering solver availability / version information
- Improve handling of version() not returning a tuple
- Fix a bug where an exception raised when creating a solver caused an
  unhandled UnknownSolver error
---
 pyomo/scripting/driver_help.py | 64 ++++++++++++++++++++--------------
 1 file changed, 37 insertions(+), 27 deletions(-)

diff --git a/pyomo/scripting/driver_help.py b/pyomo/scripting/driver_help.py
index 38d1a4c16bf..ce7c36932ec 100644
--- a/pyomo/scripting/driver_help.py
+++ b/pyomo/scripting/driver_help.py
@@ -20,6 +20,7 @@
 
 import pyomo.common
 from pyomo.common.collections import Bunch
+from pyomo.common.tee import capture_output
 import pyomo.scripting.pyomo_parser
 
 logger = logging.getLogger('pyomo.solvers')
@@ -235,33 +236,42 @@ def help_solvers():
     try:
         # Disable warnings
         logging.disable(logging.WARNING)
-        for s in solver_list:
-            # Create a solver, and see if it is available
-            with pyomo.opt.SolverFactory(s) as opt:
-                ver = ''
-                if opt.available(False):
-                    avail = '-'
-                    if opt.license_is_valid():
-                        avail = '+'
-                    try:
-                        ver = opt.version()
-                        if ver:
-                            while len(ver) > 2 and ver[-1] == 0:
-                                ver = ver[:-1]
-                            ver = '.'.join(str(v) for v in ver)
-                        else:
-                            ver = ''
-                    except (AttributeError, NameError):
-                        pass
-                elif s == 'py' or (hasattr(opt, "_metasolver") and opt._metasolver):
-                    # py is a metasolver, but since we don't specify a subsolver
-                    # for this test, opt is actually an UnknownSolver, so we
-                    # can't try to get the _metasolver attribute from it.
-                    # Also, default to False if the attribute isn't implemented
-                    avail = '*'
-                else:
-                    avail = ''
-                _data.append((avail, s, ver, pyomo.opt.SolverFactory.doc(s)))
+        # suppress ALL output
+        with capture_output(capture_fd=True):
+            for s in solver_list:
+                # Create a solver, and see if it is available
+                with pyomo.opt.SolverFactory(s) as opt:
+                    ver = ''
+                    if opt.available(False):
+                        avail = '-'
+                        if opt.license_is_valid():
+                            avail = '+'
+                        try:
+                            ver = opt.version()
+                            if isinstance(ver, str):
+                                pass
+                            elif ver:
+                                while len(ver) > 2 and ver[-1] == 0:
+                                    ver = ver[:-1]
+                                ver = '.'.join(str(v) for v in ver)
+                            else:
+                                ver = ''
+                        except (AttributeError, NameError):
+                            pass
+                    elif isinstance(s, UnknownSolver):
+                        # We can get here if creating a registered
+                        # solver failed (i.e., an exception was raised
+                        # in __init__)
+                        avail = ''
+                    elif s == 'py' or (hasattr(opt, "_metasolver") and opt._metasolver):
+                        # py is a metasolver, but since we don't specify a subsolver
+                        # for this test, opt is actually an UnknownSolver, so we
+                        # can't try to get the _metasolver attribute from it.
+                        # Also, default to False if the attribute isn't implemented
+                        avail = '*'
+                    else:
+                        avail = ''
+                    _data.append((avail, s, ver, pyomo.opt.SolverFactory.doc(s)))
     finally:
         # Reset logging level
         logging.disable(logging.NOTSET)

From 8a4dcf1d1b32969681f877de3f1b3a8bb3d464d3 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:17:01 -0600
Subject: [PATCH 59/69] Guard against platforms missing lsb_release

---
 pyomo/common/tests/test_download.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/pyomo/common/tests/test_download.py b/pyomo/common/tests/test_download.py
index 4fde029f1b1..4ee781d5738 100644
--- a/pyomo/common/tests/test_download.py
+++ b/pyomo/common/tests/test_download.py
@@ -22,7 +22,7 @@
 import pyomo.common.envvar as envvar
 
 from pyomo.common import DeveloperError
-from pyomo.common.fileutils import this_file
+from pyomo.common.fileutils import this_file, Executable
 from pyomo.common.download import FileDownloader, distro_available
 from pyomo.common.log import LoggingIntercept
 from pyomo.common.tee import capture_output
@@ -173,7 +173,8 @@ def test_get_os_version(self):
                 self.assertTrue(v.replace('.', '').startswith(dist_ver))
 
             if (
-                subprocess.run(
+                Executable('lsb_release').available()
+                and subprocess.run(
                     ['lsb_release'],
                     stdout=subprocess.DEVNULL,
                     stderr=subprocess.DEVNULL,

From 0a87a535ef1c07ca00077b904a1d631d7092cf00 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:18:59 -0600
Subject: [PATCH 60/69] Add missing import

---
 pyomo/repn/plugins/__init__.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/pyomo/repn/plugins/__init__.py b/pyomo/repn/plugins/__init__.py
index d3804c55106..ffe131b9b8b 100644
--- a/pyomo/repn/plugins/__init__.py
+++ b/pyomo/repn/plugins/__init__.py
@@ -37,6 +37,8 @@ def load():
 
 def activate_writer_version(name, ver):
     """DEBUGGING TOOL to switch the "default" writer implementation"""
+    from pyomo.opt import WriterFactory
+
     doc = WriterFactory.doc(name)
     WriterFactory.unregister(name)
     WriterFactory.register(name, doc)(WriterFactory.get_class(f'{name}_v{ver}'))

From 353fb6c66cce5bfdb2bbded0801d5166f309ebb6 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:19:41 -0600
Subject: [PATCH 61/69] NFC: apply black

---
 pyomo/contrib/appsi/solvers/ipopt.py         | 3 ++-
 pyomo/contrib/parmest/tests/test_examples.py | 4 ++--
 pyomo/contrib/parmest/tests/test_parmest.py  | 8 ++------
 3 files changed, 6 insertions(+), 9 deletions(-)

diff --git a/pyomo/contrib/appsi/solvers/ipopt.py b/pyomo/contrib/appsi/solvers/ipopt.py
index 97e8122fe78..af40d2e88d2 100644
--- a/pyomo/contrib/appsi/solvers/ipopt.py
+++ b/pyomo/contrib/appsi/solvers/ipopt.py
@@ -571,9 +571,10 @@ def get_reduced_costs(
     def has_linear_solver(self, linear_solver):
         import pyomo.core as AML
         from pyomo.common.tee import capture_output
+
         m = AML.ConcreteModel()
         m.x = AML.Var()
-        m.o = AML.Objective(expr=(m.x-2)**2)
+        m.o = AML.Objective(expr=(m.x - 2) ** 2)
         with capture_output() as OUT:
             solver = self.__class__()
             solver.config.stream_solver = True
diff --git a/pyomo/contrib/parmest/tests/test_examples.py b/pyomo/contrib/parmest/tests/test_examples.py
index 552c568cc7c..3b0c869affa 100644
--- a/pyomo/contrib/parmest/tests/test_examples.py
+++ b/pyomo/contrib/parmest/tests/test_examples.py
@@ -18,6 +18,7 @@
 ipopt_available = SolverFactory("ipopt").available()
 pynumero_ASL_available = AmplInterface.available()
 
+
 @unittest.skipIf(
     not parmest.parmest_available,
     "Cannot test parmest: required dependencies are missing",
@@ -71,8 +72,7 @@ def test_likelihood_ratio_example(self):
 @unittest.skipUnless(pynumero_ASL_available, "test requires libpynumero_ASL")
 @unittest.skipUnless(ipopt_available, "The 'ipopt' solver is not available")
 @unittest.skipUnless(
-    parmest.parmest_available,
-    "Cannot test parmest: required dependencies are missing",
+    parmest.parmest_available, "Cannot test parmest: required dependencies are missing"
 )
 class TestReactionKineticsExamples(unittest.TestCase):
     @classmethod
diff --git a/pyomo/contrib/parmest/tests/test_parmest.py b/pyomo/contrib/parmest/tests/test_parmest.py
index fbe4290dd10..52b7cd390e8 100644
--- a/pyomo/contrib/parmest/tests/test_parmest.py
+++ b/pyomo/contrib/parmest/tests/test_parmest.py
@@ -22,12 +22,7 @@
 import pyomo.environ as pyo
 import pyomo.dae as dae
 
-from pyomo.common.dependencies import (
-    numpy as np,
-    pandas as pd,
-    scipy,
-    matplotlib,
-)
+from pyomo.common.dependencies import numpy as np, pandas as pd, scipy, matplotlib
 from pyomo.common.fileutils import this_file_dir
 from pyomo.contrib.parmest.experiment import Experiment
 from pyomo.contrib.pynumero.asl import AmplInterface
@@ -38,6 +33,7 @@
 pynumero_ASL_available = AmplInterface.available()
 testdir = this_file_dir()
 
+
 @unittest.skipIf(
     not parmest.parmest_available,
     "Cannot test parmest: required dependencies are missing",

From 7ae66e33d29236caf37413106632d0f954fcfe05 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 17:36:25 -0600
Subject: [PATCH 62/69] fix UnknownSolver logic, missing import

---
 pyomo/scripting/driver_help.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/pyomo/scripting/driver_help.py b/pyomo/scripting/driver_help.py
index ce7c36932ec..45fdc711137 100644
--- a/pyomo/scripting/driver_help.py
+++ b/pyomo/scripting/driver_help.py
@@ -258,16 +258,18 @@ def help_solvers():
                                 ver = ''
                         except (AttributeError, NameError):
                             pass
-                    elif isinstance(s, UnknownSolver):
+                    elif s == 'py':
+                        # py is a metasolver, but since we don't specify a subsolver
+                        # for this test, opt is actually an UnknownSolver, so we
+                        # can't try to get the _metasolver attribute from it.
+                        avail = '*'
+                    elif isinstance(s, pyomo.opt.solvers.UnknownSolver):
                         # We can get here if creating a registered
                         # solver failed (i.e., an exception was raised
                         # in __init__)
                         avail = ''
-                    elif s == 'py' or (hasattr(opt, "_metasolver") and opt._metasolver):
-                        # py is a metasolver, but since we don't specify a subsolver
-                        # for this test, opt is actually an UnknownSolver, so we
-                        # can't try to get the _metasolver attribute from it.
-                        # Also, default to False if the attribute isn't implemented
+                    elif getattr(opt, "_metasolver", False):
+                        # Note: default to False if the attribute isn't implemented
                         avail = '*'
                     else:
                         avail = ''

From 0f79ae044118e3e7607bf09a89ea9bacd42bb05c Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 22:27:46 -0600
Subject: [PATCH 63/69] Updating baseline due to new ipopt has_linear_solver
 method

---
 pyomo/contrib/solver/tests/unit/test_ipopt.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyomo/contrib/solver/tests/unit/test_ipopt.py b/pyomo/contrib/solver/tests/unit/test_ipopt.py
index cc459245506..9769eadecae 100644
--- a/pyomo/contrib/solver/tests/unit/test_ipopt.py
+++ b/pyomo/contrib/solver/tests/unit/test_ipopt.py
@@ -84,6 +84,7 @@ def test_class_member_list(self):
             'CONFIG',
             'config',
             'available',
+            'has_linear_solver',
             'is_persistent',
             'solve',
             'version',

From 2629cfdcdef022e11fa9458db883865e065ff38c Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Tue, 6 Aug 2024 23:41:22 -0600
Subject: [PATCH 64/69] Additional test guard for ipopt availability

---
 pyomo/contrib/doe/tests/test_fim_doe.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/pyomo/contrib/doe/tests/test_fim_doe.py b/pyomo/contrib/doe/tests/test_fim_doe.py
index d9a8d60fdb4..8891bd072a4 100644
--- a/pyomo/contrib/doe/tests/test_fim_doe.py
+++ b/pyomo/contrib/doe/tests/test_fim_doe.py
@@ -35,6 +35,9 @@
     VariablesWithIndices,
 )
 from pyomo.contrib.doe.examples.reactor_kinetics import create_model, disc_for_measure
+from pyomo.environ import SolverFactory
+
+ipopt_available = SolverFactory("ipopt").available()
 
 
 class TestMeasurementError(unittest.TestCase):
@@ -196,6 +199,7 @@ def test(self):
 
 
 @unittest.skipIf(not numpy_available, "Numpy is not available")
+@unittest.skipIf(not ipopt_available, "Numpy is not available")
 class TestPriorFIMError(unittest.TestCase):
     def test(self):
         # Control time set [h]

From 12930fe61af1f48ac5500db4fda2e207b4c88a32 Mon Sep 17 00:00:00 2001
From: Miranda Mundt <mmundt@sandia.gov>
Date: Thu, 8 Aug 2024 06:48:01 -0600
Subject: [PATCH 65/69] Bug fix: Issue 3336

---
 pyomo/opt/base/solvers.py           | 4 ++--
 pyomo/opt/tests/base/test_solver.py | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyomo/opt/base/solvers.py b/pyomo/opt/base/solvers.py
index c0698165603..4ffef7e7cac 100644
--- a/pyomo/opt/base/solvers.py
+++ b/pyomo/opt/base/solvers.py
@@ -470,8 +470,8 @@ def set_results_format(self, format):
         Set the current results format (if it's valid for the current
         problem format).
         """
-        if (self._problem_format in self._valid_results_formats) and (
-            format in self._valid_results_formats[self._problem_format]
+        if (self._problem_format in self._valid_result_formats) and (
+            format in self._valid_result_formats[self._problem_format]
         ):
             self._results_format = format
         else:
diff --git a/pyomo/opt/tests/base/test_solver.py b/pyomo/opt/tests/base/test_solver.py
index 8ffc647804d..919e9375f60 100644
--- a/pyomo/opt/tests/base/test_solver.py
+++ b/pyomo/opt/tests/base/test_solver.py
@@ -109,7 +109,7 @@ def test_set_problem_format(self):
     def test_set_results_format(self):
         opt = pyomo.opt.SolverFactory("stest1")
         opt._valid_problem_formats = ['a']
-        opt._valid_results_formats = {'a': 'b'}
+        opt._valid_result_formats = {'a': 'b'}
         self.assertEqual(opt.problem_format(), None)
         try:
             opt.set_results_format('b')

From a4946df48905ce94f3d8fc840374a03d30f5a599 Mon Sep 17 00:00:00 2001
From: Miranda Mundt <55767766+mrmundt@users.noreply.github.com>
Date: Thu, 8 Aug 2024 12:53:45 -0600
Subject: [PATCH 66/69] Update pyomo/contrib/doe/tests/test_fim_doe.py

---
 pyomo/contrib/doe/tests/test_fim_doe.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/contrib/doe/tests/test_fim_doe.py b/pyomo/contrib/doe/tests/test_fim_doe.py
index 8891bd072a4..9cae2fe6278 100644
--- a/pyomo/contrib/doe/tests/test_fim_doe.py
+++ b/pyomo/contrib/doe/tests/test_fim_doe.py
@@ -199,7 +199,7 @@ def test(self):
 
 
 @unittest.skipIf(not numpy_available, "Numpy is not available")
-@unittest.skipIf(not ipopt_available, "Numpy is not available")
+@unittest.skipIf(not ipopt_available, "ipopt is not available")
 class TestPriorFIMError(unittest.TestCase):
     def test(self):
         # Control time set [h]

From 0855619fb68261e87ab3f03791032fa9b610e931 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 9 Aug 2024 12:59:05 -0600
Subject: [PATCH 67/69] Add testing of ipopt.has_linear_solver()

---
 pyomo/contrib/appsi/tests/test_ipopt.py       | 42 +++++++++++++++++++
 pyomo/contrib/solver/tests/unit/test_ipopt.py | 23 ++++++++++
 pyomo/solvers/plugins/solvers/IPOPT.py        |  8 +++-
 pyomo/solvers/tests/checks/test_ipopt.py      | 42 +++++++++++++++++++
 4 files changed, 113 insertions(+), 2 deletions(-)
 create mode 100644 pyomo/contrib/appsi/tests/test_ipopt.py
 create mode 100644 pyomo/solvers/tests/checks/test_ipopt.py

diff --git a/pyomo/contrib/appsi/tests/test_ipopt.py b/pyomo/contrib/appsi/tests/test_ipopt.py
new file mode 100644
index 00000000000..b3697b9b233
--- /dev/null
+++ b/pyomo/contrib/appsi/tests/test_ipopt.py
@@ -0,0 +1,42 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2024
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+from pyomo.common import unittest
+from pyomo.contrib.appsi.solvers import ipopt
+
+
+ipopt_available = ipopt.Ipopt().available()
+
+
+@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
+class TestIpoptInterface(unittest.TestCase):
+    def test_has_linear_solver(self):
+        opt = ipopt.Ipopt()
+        self.assertTrue(
+            any(
+                map(
+                    opt.has_linear_solver,
+                    [
+                        'mumps',
+                        'ma27',
+                        'ma57',
+                        'ma77',
+                        'ma86',
+                        'ma97',
+                        'pardiso',
+                        'pardisomkl',
+                        'spral',
+                        'wsmp',
+                    ],
+                )
+            )
+        )
+        self.assertFalse(opt.has_linear_solver('bogus_linear_solver'))
diff --git a/pyomo/contrib/solver/tests/unit/test_ipopt.py b/pyomo/contrib/solver/tests/unit/test_ipopt.py
index 9769eadecae..27a80feede0 100644
--- a/pyomo/contrib/solver/tests/unit/test_ipopt.py
+++ b/pyomo/contrib/solver/tests/unit/test_ipopt.py
@@ -168,6 +168,29 @@ def test_write_options_file(self):
                 data = f.readlines()
                 self.assertEqual(len(data), len(list(opt.config.solver_options.keys())))
 
+    def test_has_linear_solver(self):
+        opt = ipopt.Ipopt()
+        self.assertTrue(
+            any(
+                map(
+                    opt.has_linear_solver,
+                    [
+                        'mumps',
+                        'ma27',
+                        'ma57',
+                        'ma77',
+                        'ma86',
+                        'ma97',
+                        'pardiso',
+                        'pardisomkl',
+                        'spral',
+                        'wsmp',
+                    ],
+                )
+            )
+        )
+        self.assertFalse(opt.has_linear_solver('bogus_linear_solver'))
+
     def test_create_command_line(self):
         opt = ipopt.Ipopt()
         # No custom options, no file created. Plain and simple.
diff --git a/pyomo/solvers/plugins/solvers/IPOPT.py b/pyomo/solvers/plugins/solvers/IPOPT.py
index be0f143ea46..82dcfdb75a0 100644
--- a/pyomo/solvers/plugins/solvers/IPOPT.py
+++ b/pyomo/solvers/plugins/solvers/IPOPT.py
@@ -14,6 +14,7 @@
 
 from pyomo.common import Executable
 from pyomo.common.collections import Bunch
+from pyomo.common.errors import ApplicationError
 from pyomo.common.tee import capture_output
 from pyomo.common.tempfiles import TempfileManager
 
@@ -215,6 +216,9 @@ def has_linear_solver(self, linear_solver):
         m = AML.ConcreteModel()
         m.x = AML.Var()
         m.o = AML.Objective(expr=(m.x - 2) ** 2)
-        with capture_output() as OUT:
-            self.solve(m, tee=True, options={'linear_solver': linear_solver})
+        try:
+            with capture_output() as OUT:
+                self.solve(m, tee=True, options={'linear_solver': linear_solver})
+        except ApplicationError:
+            return False
         return 'running with linear solver' in OUT.getvalue()
diff --git a/pyomo/solvers/tests/checks/test_ipopt.py b/pyomo/solvers/tests/checks/test_ipopt.py
new file mode 100644
index 00000000000..b7d00c35a6f
--- /dev/null
+++ b/pyomo/solvers/tests/checks/test_ipopt.py
@@ -0,0 +1,42 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2024
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+from pyomo.common import unittest
+from pyomo.solvers.plugins.solvers import IPOPT
+import pyomo.environ
+
+ipopt_available = IPOPT.IPOPT().available()
+
+
+@unittest.skipIf(not ipopt_available, "The 'ipopt' command is not available")
+class TestIpoptInterface(unittest.TestCase):
+    def test_has_linear_solver(self):
+        opt = IPOPT.IPOPT()
+        self.assertTrue(
+            any(
+                map(
+                    opt.has_linear_solver,
+                    [
+                        'mumps',
+                        'ma27',
+                        'ma57',
+                        'ma77',
+                        'ma86',
+                        'ma97',
+                        'pardiso',
+                        'pardisomkl',
+                        'spral',
+                        'wsmp',
+                    ],
+                )
+            )
+        )
+        self.assertFalse(opt.has_linear_solver('bogus_linear_solver'))

From cbd702dadd80e218a8e8597e511c57b2d0b449a9 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 9 Aug 2024 13:14:13 -0600
Subject: [PATCH 68/69] Add active_writer_version, test writer activation
 functions

---
 pyomo/repn/plugins/__init__.py   | 15 ++++++++++
 pyomo/repn/tests/test_plugins.py | 50 ++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+)
 create mode 100644 pyomo/repn/tests/test_plugins.py

diff --git a/pyomo/repn/plugins/__init__.py b/pyomo/repn/plugins/__init__.py
index ffe131b9b8b..4029f44a03d 100644
--- a/pyomo/repn/plugins/__init__.py
+++ b/pyomo/repn/plugins/__init__.py
@@ -42,3 +42,18 @@ def activate_writer_version(name, ver):
     doc = WriterFactory.doc(name)
     WriterFactory.unregister(name)
     WriterFactory.register(name, doc)(WriterFactory.get_class(f'{name}_v{ver}'))
+
+
+def active_writer_version(name):
+    """DEBUGGING TOOL to switch the "default" writer implementation"""
+    from pyomo.opt import WriterFactory
+
+    ref = WriterFactory.get_class(name)
+    ver = 1
+    try:
+        while 1:
+            if WriterFactory.get_class(f'{name}_v{ver}') is ref:
+                return ver
+            ver += 1
+    except KeyError:
+        return None
diff --git a/pyomo/repn/tests/test_plugins.py b/pyomo/repn/tests/test_plugins.py
new file mode 100644
index 00000000000..fa6026ea74e
--- /dev/null
+++ b/pyomo/repn/tests/test_plugins.py
@@ -0,0 +1,50 @@
+#  ___________________________________________________________________________
+#
+#  Pyomo: Python Optimization Modeling Objects
+#  Copyright (c) 2008-2024
+#  National Technology and Engineering Solutions of Sandia, LLC
+#  Under the terms of Contract DE-NA0003525 with National Technology and
+#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
+#  rights in this software.
+#  This software is distributed under the 3-clause BSD License.
+#  ___________________________________________________________________________
+
+from pyomo.common import unittest
+
+from pyomo.opt import WriterFactory
+from pyomo.repn.plugins import activate_writer_version, active_writer_version
+
+import pyomo.environ
+
+
+class TestPlugins(unittest.TestCase):
+    def test_active(self):
+        with self.assertRaises(KeyError):
+            active_writer_version('nonexistant_writer')
+        ver = active_writer_version('lp')
+        self.assertIs(
+            WriterFactory.get_class('lp'), WriterFactory.get_class(f'lp_v{ver}')
+        )
+
+        class TMP(object):
+            pass
+
+        WriterFactory.register('test_writer')(TMP)
+        try:
+            self.assertIsNone(active_writer_version('test_writer'))
+        finally:
+            WriterFactory.unregister('test_writer')
+
+    def test_activate(self):
+        ver = active_writer_version('lp')
+        try:
+            activate_writer_version('lp', 2)
+            self.assertIs(
+                WriterFactory.get_class('lp'), WriterFactory.get_class(f'lp_v2')
+            )
+            activate_writer_version('lp', 1)
+            self.assertIs(
+                WriterFactory.get_class('lp'), WriterFactory.get_class(f'lp_v1')
+            )
+        finally:
+            activate_writer_version('lp', ver)

From 658dab36b72f3bd18759dad286f5e20b4ec936e7 Mon Sep 17 00:00:00 2001
From: John Siirola <jsiirola@users.noreply.github.com>
Date: Fri, 9 Aug 2024 13:29:54 -0600
Subject: [PATCH 69/69] NFC: fix typo

---
 pyomo/repn/tests/test_plugins.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyomo/repn/tests/test_plugins.py b/pyomo/repn/tests/test_plugins.py
index fa6026ea74e..1152131f6b6 100644
--- a/pyomo/repn/tests/test_plugins.py
+++ b/pyomo/repn/tests/test_plugins.py
@@ -20,7 +20,7 @@
 class TestPlugins(unittest.TestCase):
     def test_active(self):
         with self.assertRaises(KeyError):
-            active_writer_version('nonexistant_writer')
+            active_writer_version('nonexistent_writer')
         ver = active_writer_version('lp')
         self.assertIs(
             WriterFactory.get_class('lp'), WriterFactory.get_class(f'lp_v{ver}')