Merge branch 'nonlinear-to-pwl' into bashar_pwl

.github/workflows/test_pr_and_main.yml

@@ -29,7 +29,7 @@ defaults:
 env:
   PYTHONWARNINGS: ignore::UserWarning
   PYTHON_CORE_PKGS: wheel
-  PYPI_ONLY: z3-solver
+  PYPI_ONLY: z3-solver linear-tree
   PYPY_EXCLUDE: scipy numdifftools seaborn statsmodels
   CACHE_VER: v221013.1
   NEOS_EMAIL: [email protected]

pyomo/contrib/piecewise/tests/test_nonlinear_to_pwl.py

@@ -9,7 +9,11 @@
 #  This software is distributed under the 3-clause BSD License.
 #  ___________________________________________________________________________
 
+from io import StringIO
+import logging
+
 from pyomo.common.dependencies import attempt_import, scipy_available, numpy_available
+from pyomo.common.log import LoggingIntercept
 import pyomo.common.unittest as unittest
 from pyomo.contrib.piecewise import PiecewiseLinearFunction
 from pyomo.contrib.piecewise.transform.nonlinear_to_pwl import (
@@ -23,9 +27,11 @@
 )
 from pyomo.core.expr.numeric_expr import SumExpression
 from pyomo.environ import (
+    Binary,
     ConcreteModel,
     Var,
     Constraint,
+    Integers,
     TransformationFactory,
     log,
     Objective,
@@ -303,6 +309,102 @@ def test_do_not_additively_decompose_below_min_dimension(self):
         # This is only approximated by one pwlf:
         self.assertIsInstance(transformed_c.body, _ExpressionData)
 
+    @unittest.skipUnless(numpy_available, "Numpy is not available")
+    def test_uniform_sampling_discrete_vars(self):
+        m = ConcreteModel()
+        m.x = Var(['rocky', 'bullwinkle'], domain=Binary)
+        m.y = Var(domain=Integers, bounds=(0, 5))
+        m.c = Constraint(expr=m.x['rocky'] * m.x['bullwinkle'] + m.y <= 4)
+
+        n_to_pwl = TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
+        output = StringIO()
+        with LoggingIntercept(output, 'pyomo.core', logging.WARNING):
+            n_to_pwl.apply_to(
+                m,
+                num_points=3,
+                additively_decompose=False,
+                domain_partitioning_method=DomainPartitioningMethod.UNIFORM_GRID,
+            )
+        # No warnings (this is to check that we aren't emitting a bunch of
+        # warnings about setting variables outside of their domains)
+        self.assertEqual(output.getvalue().strip(), "")
+
+        transformed_c = n_to_pwl.get_transformed_component(m.c)
+        pwlf = transformed_c.body.expr.pw_linear_function
+
+        # should sample 0, 1 for th m.x's
+        # should sample 0, 2, 5 for m.y (because of half to even rounding (*sigh*))
+        points = set(pwlf._points)
+        self.assertEqual(len(points), 12)
+        for x in [0, 1]:
+            for y in [0, 1]:
+                for z in [0, 2, 5]:
+                    self.assertIn((x, y, z), points)
+
+    @unittest.skipUnless(numpy_available, "Numpy is not available")
+    def test_uniform_sampling_discrete_vars(self):
+        m = ConcreteModel()
+        m.x = Var(['rocky', 'bullwinkle'], domain=Binary)
+        m.y = Var(domain=Integers, bounds=(0, 5))
+        m.c = Constraint(expr=m.x['rocky'] * m.x['bullwinkle'] + m.y <= 4)
+
+        n_to_pwl = TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
+        output = StringIO()
+        with LoggingIntercept(output, 'pyomo.core', logging.WARNING):
+            n_to_pwl.apply_to(
+                m,
+                num_points=3,
+                additively_decompose=False,
+                domain_partitioning_method=DomainPartitioningMethod.UNIFORM_GRID,
+            )
+        # No warnings (this is to check that we aren't emitting a bunch of
+        # warnings about setting variables outside of their domains)
+        self.assertEqual(output.getvalue().strip(), "")
+
+        transformed_c = n_to_pwl.get_transformed_component(m.c)
+        pwlf = transformed_c.body.expr.pw_linear_function
+
+        # should sample 0, 1 for th m.x's
+        # should sample 0, 2, 5 for m.y (because of half to even rounding (*sigh*))
+        points = set(pwlf._points)
+        self.assertEqual(len(points), 12)
+        for x in [0, 1]:
+            for y in [0, 1]:
+                for z in [0, 2, 5]:
+                    self.assertIn((x, y, z), points)
+
+    @unittest.skipUnless(numpy_available, "Numpy is not available")
+    def test_random_sampling_discrete_vars(self):
+        m = ConcreteModel()
+        m.x = Var(['rocky', 'bullwinkle'], domain=Binary)
+        m.y = Var(domain=Integers, bounds=(0, 5))
+        m.c = Constraint(expr=m.x['rocky'] * m.x['bullwinkle'] + m.y <= 4)
+
+        n_to_pwl = TransformationFactory('contrib.piecewise.nonlinear_to_pwl')
+        output = StringIO()
+        with LoggingIntercept(output, 'pyomo.core', logging.WARNING):
+            n_to_pwl.apply_to(
+                m,
+                num_points=3,
+                additively_decompose=False,
+                domain_partitioning_method=DomainPartitioningMethod.RANDOM_GRID,
+            )
+        # No warnings (this is to check that we aren't emitting a bunch of
+        # warnings about setting variables outside of their domains)
+        self.assertEqual(output.getvalue().strip(), "")
+
+        transformed_c = n_to_pwl.get_transformed_component(m.c)
+        pwlf = transformed_c.body.expr.pw_linear_function
+
+        # should sample 0, 1 for th m.x's
+        # Happen to get 0, 1, 5 for m.y
+        points = set(pwlf._points)
+        self.assertEqual(len(points), 12)
+        for x in [0, 1]:
+            for y in [0, 1]:
+                for z in [0, 1, 5]:
+                    self.assertIn((x, y, z), points)
+
 
 class TestNonlinearToPWL_2D(unittest.TestCase):
     def make_paraboloid_model(self):

pyomo/contrib/piecewise/transform/nonlinear_to_pwl.py

@@ -93,19 +93,30 @@ def __init__(self):
 def _get_random_point_grid(bounds, n, func, config, seed=42):
     # Generate randomized grid of points
     linspaces = []
-    for lb, ub in bounds:
-        np.random.seed(seed)
-        linspaces.append(np.random.uniform(lb, ub, n))
+    np.random.seed(seed)
+    for (lb, ub), is_integer in bounds:
+        if not is_integer:
+            linspaces.append(np.random.uniform(lb, ub, n))
+        else:
+            size = min(n, ub - lb + 1)
+            linspaces.append(
+                np.random.choice(range(lb, ub + 1), size=size, replace=False)
+            )
     return list(itertools.product(*linspaces))
 
 
 def _get_uniform_point_grid(bounds, n, func, config):
     # Generate non-randomized grid of points
     linspaces = []
-    for lb, ub in bounds:
-        # Issues happen when exactly using the boundary
-        nudge = (ub - lb) * 1e-4
-        linspaces.append(np.linspace(lb + nudge, ub - nudge, n))
+    for (lb, ub), is_integer in bounds:
+        if not is_integer:
+            # Issues happen when exactly using the boundary
+            nudge = (ub - lb) * 1e-4
+            linspaces.append(np.linspace(lb + nudge, ub - nudge, n))
+        else:
+            size = min(n, ub - lb + 1)
+            pts = np.linspace(lb, ub, size)
+            linspaces.append(np.array([round(i) for i in pts]))
     return list(itertools.product(*linspaces))
 
 
@@ -159,8 +170,8 @@ def _get_pwl_function_approximation(func, config, bounds):
     func: function to approximate
     config: ConfigDict for transformation, specifying domain_partitioning_method,
        num_points, and max_depth (if using linear trees)
-    bounds: list of tuples giving upper and lower bounds for each of func's
-       arguments
+    bounds: list of tuples giving upper and lower bounds and a boolean indicating
+       if the variable's domain is discrete or not, for each of func's arguments
     """
     method = config.domain_partitioning_method
     n = config.num_points
@@ -195,8 +206,8 @@ def _generate_bound_points(leaves, bounds):
         for pt in [lower_corner_list, upper_corner_list]:
             for i in range(len(pt)):
                 # clamp within bounds range
-                pt[i] = max(pt[i], bounds[i][0])
-                pt[i] = min(pt[i], bounds[i][1])
+                pt[i] = max(pt[i], bounds[i][0][0])
+                pt[i] = min(pt[i], bounds[i][0][1])
 
         if tuple(lower_corner_list) not in bound_points:
             bound_points.append(tuple(lower_corner_list))
@@ -206,7 +217,7 @@ def _generate_bound_points(leaves, bounds):
     # This process should have gotten every interior bound point. However, all
     # but two of the corners of the overall bounding box should have been
     # missed. Let's fix that now.
-    for outer_corner in itertools.product(*bounds):
+    for outer_corner in itertools.product(*[b[0] for b in bounds]):
         if outer_corner not in bound_points:
             bound_points.append(outer_corner)
     return bound_points
@@ -296,9 +307,9 @@ def _reassign_none_bounds(leaves, input_bounds):
     for l in L:
         for f in features:
             if leaves[l]['bounds'][f][0] == None:
-                leaves[l]['bounds'][f][0] = input_bounds[f][0]
+                leaves[l]['bounds'][f][0] = input_bounds[f][0][0]
             if leaves[l]['bounds'][f][1] == None:
-                leaves[l]['bounds'][f][1] = input_bounds[f][1]
+                leaves[l]['bounds'][f][1] = input_bounds[f][0][1]
     return leaves
 
 
@@ -615,7 +626,7 @@ def _get_bounds_list(self, var_list, obj):
                     "at least one bound" % (v.name, obj.name)
                 )
             else:
-                bounds.append(v.bounds)
+                bounds.append((v.bounds, v.is_integer()))
         return bounds
 
     def _needs_approximating(self, expr, approximate_quadratic):

pyomo/core/base/constraint.py

@@ -179,7 +179,7 @@ def __call__(self, exception=True):
             body = value(self.body, exception=exception)
         return body
 
-    def to_bounded_expression(self):
+    def to_bounded_expression(self, evaluate_bounds=False):
         """Convert this constraint to a tuple of 3 expressions (lb, body, ub)
 
         This method "standardizes" the expression into a 3-tuple of
@@ -195,6 +195,13 @@ def to_bounded_expression(self):
         extension, the result) can change after fixing / unfixing
         :py:class:`Var` objects.
 
+        Parameters
+        ----------
+        evaluate_bounds: bool
+
+            If True, then the lower and upper bounds will be evaluated
+            to a finite numeric constant or None.
+
         Raises
         ------
 
@@ -226,15 +233,36 @@ def to_bounded_expression(self):
                     "variable upper bound.  Cannot normalize the "
                     "constraint or send it to a solver."
                 )
-            return ans
-        elif expr is not None:
+        elif expr is None:
+            ans = None, None, None
+        else:
             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
+                ans = rhs if expr.__class__ is EqualityExpression else None, lhs, rhs
+            elif lhs.__class__ in native_types or not lhs.is_potentially_variable():
+                ans = lhs, rhs, lhs if expr.__class__ is EqualityExpression else None
+            else:
+                ans = 0 if expr.__class__ is EqualityExpression else None, lhs - rhs, 0
+
+        if evaluate_bounds:
+            lb, body, ub = ans
+            return self._evaluate_bound(lb, True), body, self._evaluate_bound(ub, False)
+        return ans
+
+    def _evaluate_bound(self, bound, is_lb):
+        if bound is None:
+            return None
+        if bound.__class__ not in native_numeric_types:
+            bound = float(value(bound))
+        # Note that "bound != bound" catches float('nan')
+        if bound in _nonfinite_values or bound != bound:
+            if bound == (-_inf if is_lb else _inf):
+                return None
+            raise ValueError(
+                f"Constraint '{self.name}' created with an invalid non-finite "
+                f"{'lower' if is_lb else 'upper'} bound ({bound})."
+            )
+        return bound
 
     @property
     def body(self):
@@ -291,38 +319,12 @@ def upper(self):
     @property
     def lb(self):
         """Access the value of the lower bound of a constraint expression."""
-        bound = self.to_bounded_expression()[0]
-        if bound is None:
-            return None
-        if bound.__class__ not in native_numeric_types:
-            bound = float(value(bound))
-        # Note that "bound != bound" catches float('nan')
-        if bound in _nonfinite_values or bound != bound:
-            if bound == -_inf:
-                return None
-            raise ValueError(
-                f"Constraint '{self.name}' created with an invalid non-finite "
-                f"lower bound ({bound})."
-            )
-        return bound
+        return self._evaluate_bound(self.to_bounded_expression()[0], True)
 
     @property
     def ub(self):
         """Access the value of the upper bound of a constraint expression."""
-        bound = self.to_bounded_expression()[2]
-        if bound is None:
-            return None
-        if bound.__class__ not in native_numeric_types:
-            bound = float(value(bound))
-        # Note that "bound != bound" catches float('nan')
-        if bound in _nonfinite_values or bound != bound:
-            if bound == _inf:
-                return None
-            raise ValueError(
-                f"Constraint '{self.name}' created with an invalid non-finite "
-                f"upper bound ({bound})."
-            )
-        return bound
+        return self._evaluate_bound(self.to_bounded_expression()[2], False)
 
     @property
     def equality(self):

pyomo/core/kernel/constraint.py

@@ -160,6 +160,17 @@ def has_ub(self):
         ub = self.ub
         return (ub is not None) and (value(ub) != float('inf'))
 
+    def to_bounded_expression(self, evaluate_bounds=False):
+        if evaluate_bounds:
+            lb = self.lb
+            if lb == -float('inf'):
+                lb = None
+            ub = self.ub
+            if ub == float('inf'):
+                ub = None
+            return lb, self.body, ub
+        return self.lower, self.body, self.upper
+
 
 class _MutableBoundsConstraintMixin(object):
     """
@@ -177,9 +188,6 @@ class _MutableBoundsConstraintMixin(object):
     # Define some of the IConstraint abstract methods
     #
 
-    def to_bounded_expression(self):
-        return self.lower, self.body, self.upper
-
     @property
     def lower(self):
         """The expression for the lower bound of the constraint"""

pyomo/repn/beta/matrix.py

@@ -587,6 +587,11 @@ def constant(self):
     # Abstract Interface (ConstraintData)
     #
 
+    def to_bounded_expression(self, evaluate_bounds=False):
+        """Access this constraint as a single expression."""
+        # Note that the bounds are always going to be floats...
+        return self.lower, self.body, self.upper
+
     @property
     def body(self):
         """Access the body of a constraint expression."""