Merge pull request #3672 from emma58/linear-dual-fix

Bugfixes for `core.lp_dual` transformation

Author: blnicho

pyomo/core/plugins/transform/lp_dual.py

@@ -24,8 +24,10 @@
     NonPositiveReals,
     maximize,
     minimize,
+    RangeSet,
     Reals,
 )
+from pyomo.core.expr.numvalue import native_numeric_types
 from pyomo.opt import WriterFactory
 from pyomo.repn.standard_repn import isclose_const
 from pyomo.util.config_domains import ComponentDataSet
@@ -110,6 +112,11 @@ def _take_dual(self, model, std_form):
                 "Model '%s' has no objective or multiple active objectives. Can "
                 "only take dual with exactly one active objective!" % model.name
             )
+        if len(std_form.columns) == 0 and std_form.c.shape[1] == 0:
+            raise ValueError(
+                f"Model '{model.name}' has no variables in the active Constraints "
+                f"or Objective."
+            )
         primal_sense = std_form.objectives[0].sense
 
         dual = ConcreteModel(name="%s dual" % model.name)
@@ -121,7 +128,29 @@ def _take_dual(self, model, std_form):
         dual_cols = range(A.shape[0])
         dual.x = Var(dual_cols, domain=NonNegativeReals)
         trans_info = dual.private_data()
+        A_csr = A.tocsr()
         for j, (primal_cons, ineq) in enumerate(std_form.rows):
+            # We need to check this constraint isn't trivial due to the
+            # parameterization, which we can detect if the row is all 0's.
+            if A_csr.indptr[j] == A_csr.indptr[j + 1]:
+                # All 0's in the coefficient matrix: check what's on the RHS
+                b = std_form.rhs[j]
+                if type(b) not in native_numeric_types:
+                    # The parameterization made this trivial. I'm not sure what's
+                    # really expected here, so maybe we just scream? Or we leave
+                    # the constraint in the model as it is written...
+                    raise ValueError(
+                        f"The primal model contains a constraint that the "
+                        f"parameterization makes trivial: '{primal_cons.name}'"
+                        f"\nPlease deactivate it or declare it on another Block "
+                        f"before taking the dual."
+                    )
+                else:
+                    # The whole constraint is trivial--it will already have been
+                    # checked compiling the standard form, so we can safely ignore
+                    # it.
+                    pass
+
             # maximize is -1 and minimize is +1 and ineq is +1 for <= and -1 for
             # >=, so we need to change domain to NonPositiveReals if the product
             # of these is +1.
@@ -141,17 +170,28 @@ def _take_dual(self, model, std_form):
                 primal_row = A.indices[j]
                 lhs += coef * dual.x[primal_row]
 
-            if primal.domain is Reals:
+            domain = primal.domain
+            lb, ub = domain.bounds()
+            # Note: the following checks the domain for continuity and compactness:
+            if not domain == RangeSet(*domain.bounds(), 0):
+                raise ValueError(
+                    f"The domain of the primal variable '{primal.name}' "
+                    f"is not continuous."
+                )
+            unrestricted = (lb is None or lb < 0) and (ub is None or ub > 0)
+            nonneg = (lb is not None) and lb >= 0
+
+            if unrestricted:
                 dual.constraints[i] = lhs == c[i]
             elif primal_sense is minimize:
-                if primal.domain is NonNegativeReals:
+                if nonneg:
                     dual.constraints[i] = lhs <= c[i]
-                else:  # primal.domain is NonPositiveReals
+                else:  # primal domain is nonpositive
                     dual.constraints[i] = lhs >= c[i]
             else:
-                if primal.domain is NonNegativeReals:
+                if nonneg:
                     dual.constraints[i] = lhs >= c[i]
-                else:  # primal.domain is NonPositiveReals
+                else:  # primal domain is nonpositive
                     dual.constraints[i] = lhs <= c[i]
             trans_info.dual_constraint[primal] = dual.constraints[i]
             trans_info.primal_var[dual.constraints[i]] = primal

pyomo/core/tests/unit/test_lp_dual.py

@@ -17,9 +17,11 @@
     Constraint,
     maximize,
     minimize,
+    NonNegativeIntegers,
     NonNegativeReals,
     NonPositiveReals,
     Objective,
+    Param,
     Reals,
     Suffix,
     TerminationCondition,
@@ -398,3 +400,45 @@ def test_dual_var_map_error(self):
             "on model 'primal'",
         ):
             thing = lp_dual.get_dual_var(m, m.c_new)
+
+    def test_parameterization_makes_constraint_trivial(self):
+        m = self.get_bilevel_model()
+        m.budgetish = Constraint(expr=m.outer[2] + m.outer[3] == 1)
+
+        lp_dual = TransformationFactory('core.lp_dual')
+        with self.assertRaisesRegex(
+            ValueError,
+            "The primal model contains a constraint that the "
+            "parameterization makes trivial: 'budgetish'"
+            "\nPlease deactivate it or declare it on another Block "
+            "before taking the dual.",
+        ):
+            dual = lp_dual.create_using(m, parameterize_wrt=m.outer)
+
+    def test_normal_trivial_constraint_error(self):
+        m = ConcreteModel()
+        m.p = Param(initialize=3, mutable=True)
+        m.x = Var(bounds=(0, 9))
+        m.c = Constraint(expr=m.x * m.p <= 8)
+        m.x.fix(2)
+
+        m.obj = Objective(expr=m.x)
+
+        lp_dual = TransformationFactory('core.lp_dual')
+        with self.assertRaisesRegex(
+            ValueError,
+            "Model 'unknown' has no variables in the active Constraints "
+            "or Objective.",
+        ):
+            dual = lp_dual.create_using(m)
+
+    def test_discrete_primal_var_error(self):
+        m = self.get_bilevel_model()
+        m.x.domain = NonNegativeIntegers
+
+        with self.assertRaisesRegex(
+            ValueError, "The domain of the primal variable 'x' is not continuous"
+        ):
+            dual = TransformationFactory('core.lp_dual').create_using(
+                m, parameterize_wrt=[m.outer, m.outer1]
+            )

pyomo/repn/parameterized.py

@@ -83,12 +83,12 @@ def _before_general_expression(visitor, child):
     def _before_var(visitor, child):
         _id = id(child)
         if _id not in visitor.var_map:
-            if child.fixed:
-                return False, (_CONSTANT, visitor.check_constant(child.value, child))
             if child in visitor.wrt:
-                # pseudo-constant
+                # pseudo-constant, and we need to leave it in even if it is fixed!
                 # We aren't treating this Var as a Var for the purposes of this walker
                 return False, (_FIXED, child)
+            if child.fixed:
+                return False, (_CONSTANT, visitor.check_constant(child.value, child))
             # This is a normal situation
             visitor.var_recorder.add(child)
         ans = visitor.Result()

pyomo/repn/plugins/parameterized_standard_form.py

@@ -224,6 +224,64 @@ def todense(self):
 
         return dense
 
+    def tocsr(self):
+        """Implements the same algorithm as scipy's csr_tocsc function from
+        sparsetools.
+        """
+        csc_data = self.data
+        row_index = self.indices
+        col_index_ptr = self.indptr
+        ncols = self.shape[1]
+
+        num_nonzeros = len(csc_data)
+        csr_data = np.empty(csc_data.shape[0], dtype=object)
+        col_index = np.empty(num_nonzeros, dtype=int)
+        # tally the nonzeros in each column
+        row_index_ptr = np.zeros(self.shape[0], dtype=int)
+        for i in row_index:
+            row_index_ptr[int(i)] += 1
+
+        # cumulative sum the tally to get the column index pointer
+        cum_sum = 0
+        for i, tally in enumerate(row_index_ptr):
+            row_index_ptr[i] = cum_sum
+            cum_sum += tally
+        # We have now initialized the row_index_ptr to the *starting* position
+        # of each column in the data vector.  Note that row_index_ptr is only
+        # num_rows long: we have ignored the last entry in the standard CSR
+        # row_index_ptr (the total number of nonzeros).  This will get resolved
+        # below when we shift this vector by one position.
+
+        # Now we are actually going to mess up what we just did while we
+        # construct the row index: We can imagine that row_index_ptr holds the
+        # position of the *next* nonzero in each column, so each time we move a
+        # data element into a column we will increment that row_index_ptr by
+        # one.  This is beautiful because by "messing up" the row_index_pointer,
+        # we are just transforming the vector of *starting* indices for each
+        # column to a vector of *ending* indices (actually 1 past the last
+        # index) of each column. Thank you, scipy.
+        for col in range(ncols):
+            for j in range(col_index_ptr[col], col_index_ptr[col + 1]):
+                row = row_index[j]
+                dest = row_index_ptr[row]
+                col_index[dest] = col
+                # Note that the data changes order because now we are looking
+                # for nonzeros through the columns rather than through the rows.
+                csr_data[dest] = csc_data[j]
+
+                row_index_ptr[row] += 1
+
+        # Fix the row index pointer by inserting 0 at the beginning. The
+        # row_index_ptr currently holds pointers to 1 past the last element of
+        # each row, which is really the starting index for the next
+        # row. Inserting the 0 (the starting index for the first column)
+        # shifts everything by one column, "converting" the vector to the
+        # starting indices of each row, and extending the vector length to
+        # num_rows + 1 (as is expected by the CSR matrix).
+        row_index_ptr = np.insert(row_index_ptr, 0, 0)
+
+        return _CSRMatrix((csr_data, col_index, row_index_ptr), self.shape)
+
     def sum_duplicates(self):
         """Implements the algorithm from scipy's csr_sum_duplicates function
         in sparsetools.

pyomo/repn/tests/test_parameterized_linear.py

@@ -293,15 +293,11 @@ def test_ANY_over_constant_division(self):
         m.x = Var()
         m.z = Var()
         m.y = Var()
-        # We will use the fixed value regardless of the fact that we aren't
-        # treating this as a Var.
         m.y.fix(1)
 
         expr = m.y + m.x + m.z + ((3 * m.z * m.x) / m.p) / m.y
         cfg = VisitorConfig()
-        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.y, m.z]).walk_expression(
-            expr
-        )
+        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.z]).walk_expression(expr)
 
         self.assertEqual(repn.multiplier, 1)
         assertExpressionsEqual(self, repn.constant, 1 + m.z)
@@ -330,17 +326,15 @@ def test_errors_propagate_nan(self):
             "\texpression: 3*z*x/p\n",
         )
         self.assertEqual(repn.multiplier, 1)
-        assertExpressionsEqual(self, repn.constant, 1 + m.z)
+        assertExpressionsEqual(self, repn.constant, m.y + m.z)
         self.assertEqual(len(repn.linear), 1)
         self.assertIsInstance(repn.linear[id(m.x)], InvalidNumber)
-        assertExpressionsEqual(self, repn.linear[id(m.x)].value, 1 + float('nan'))
+        assertExpressionsEqual(self, repn.linear[id(m.x)].value, 1 + float('nan') / m.y)
         self.assertEqual(repn.nonlinear, None)
 
         m.y.fix(None)
         expr = m.z * log(m.y) + 3
-        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.y, m.z]).walk_expression(
-            expr
-        )
+        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.z]).walk_expression(expr)
         self.assertEqual(repn.multiplier, 1)
         self.assertIsInstance(repn.constant, InvalidNumber)
         assertExpressionsEqual(self, repn.constant.value, float('nan') * m.z + 3)
@@ -523,7 +517,7 @@ def test_0_mult_nan_param(self):
         e = m.p * (m.y**2)
 
         cfg = VisitorConfig()
-        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.y]).walk_expression(e)
+        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[]).walk_expression(e)
 
         self.assertEqual(len(repn.linear), 0)
         self.assertEqual(repn.multiplier, 1)
@@ -539,11 +533,11 @@ def test_0_mult_linear_with_nan(self):
         e = m.p * (3 * m.x * m.y + m.z)
 
         cfg = VisitorConfig()
-        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.x]).walk_expression(e)
+        repn = ParameterizedLinearRepnVisitor(**cfg, wrt=[m.z]).walk_expression(e)
         self.assertEqual(len(repn.linear), 1)
         self.assertIn(id(m.y), repn.linear)
         self.assertIsInstance(repn.linear[id(m.y)], InvalidNumber)
         assertExpressionsEqual(self, repn.linear[id(m.y)].value, 0 * 3 * float('nan'))
         self.assertEqual(repn.multiplier, 1)
         self.assertIsNone(repn.nonlinear)
-        self.assertEqual(repn.constant, 0)
+        assertExpressionsEqual(self, repn.constant, 0 * m.z)

pyomo/repn/tests/test_parameterized_quadratic.py

@@ -1059,9 +1059,7 @@ def test_negation_nonlinear_wrt_y_fix_z(self):
         )
 
         cfg = VisitorConfig()
-        # note: variable fixing takes precedence over inclusion in
-        # the `wrt` list; that is tested here
-        visitor = ParameterizedQuadraticRepnVisitor(**cfg, wrt=[m.y, m.z])
+        visitor = ParameterizedQuadraticRepnVisitor(**cfg, wrt=[m.y])
         repn = visitor.walk_expression(expr)
 
         self.assertEqual(cfg.subexpr, {})
@@ -1391,7 +1389,7 @@ def test_pow_integer_fixed_var(self):
         expr = (1 + 2 * m.x + 3 * m.y) ** (m.z)
 
         cfg = VisitorConfig()
-        visitor = ParameterizedQuadraticRepnVisitor(**cfg, wrt=[m.y, m.z])
+        visitor = ParameterizedQuadraticRepnVisitor(**cfg, wrt=[m.y])
         repn = visitor.walk_expression(expr)
 
         self.assertEqual(cfg.subexpr, {})

pyomo/repn/tests/test_parameterized_standard_form.py

@@ -51,6 +51,11 @@ def test_csr_to_csc_only_data(self):
         self.assertTrue(np.all(thing.indices == np.array([0, 1, 3, 2])))
         self.assertTrue(np.all(thing.indptr == np.array([0, 1, 3, 4, 4])))
 
+        should_be_A = thing.tocsr()
+        self.assertTrue(np.all(should_be_A.data == A.data))
+        self.assertTrue(np.all(should_be_A.indices == A.indices))
+        self.assertTrue(np.all(should_be_A.indptr == A.indptr))
+
     def test_csr_to_csc_pyomo_exprs(self):
         m = ConcreteModel()
         m.x = Var()
@@ -70,6 +75,15 @@ def test_csr_to_csc_pyomo_exprs(self):
         self.assertTrue(np.all(thing.indices == np.array([0, 1, 3, 2])))
         self.assertTrue(np.all(thing.indptr == np.array([0, 1, 3, 4, 4])))
 
+        should_be_A = thing.tocsr()
+        self.assertEqual(should_be_A.data[0], 5)
+        assertExpressionsEqual(self, should_be_A.data[1], 8 * m.x)
+        assertExpressionsEqual(self, should_be_A.data[2], 3 * m.x * m.y**2)
+        self.assertEqual(should_be_A.data[3], 6)
+        self.assertEqual(should_be_A.data.shape, (4,))
+        self.assertTrue(np.all(should_be_A.indices == np.array([0, 1, 2, 1])))
+        self.assertTrue(np.all(should_be_A.indptr == np.array([0, 1, 2, 3, 4])))
+
     def test_csr_to_csc_empty_matrix(self):
         A = _CSRMatrix(([], [], [0]), [0, 4])
         thing = A.tocsc()
@@ -79,6 +93,12 @@ def test_csr_to_csc_empty_matrix(self):
         self.assertEqual(thing.shape, (0, 4))
         self.assertTrue(np.all(thing.indptr == np.zeros(5)))
 
+        should_be_A = thing.tocsr()
+        self.assertEqual(should_be_A.data.size, 0)
+        self.assertEqual(should_be_A.indices.size, 0)
+        self.assertEqual(should_be_A.shape, (0, 4))
+        self.assertTrue(np.all(should_be_A.indptr == np.zeros(5)))
+
     def test_todense(self):
         A = _CSRMatrix(([5, 8, 3, 6], [0, 1, 2, 1], [0, 1, 2, 3, 4]), [4, 4])
         dense = np.array([[5, 0, 0, 0], [0, 8, 0, 0], [0, 0, 3, 0], [0, 6, 0, 0]])