Merge pull request scipy#20044 from andyfaff/matrix

TST: _ConstraintWrapper returns a violation of the correct shape

scipy/optimize/_differentialevolution.py

@@ -1799,6 +1799,12 @@ class _ConstraintWrapper:
         Contains lower and upper bounds for the constraints --- lb and ub.
         These are converted to ndarray and have a size equal to the number of
         the constraints.
+
+    Notes
+    -----
+    _ConstraintWrapper.fun and _ConstraintWrapper.violation can get sent
+    arrays of shape (N, S) or (N,), where S is the number of vectors of shape
+    (N,) to consider constraints for.
     """
     def __init__(self, constraint, x0):
         self.constraint = constraint
@@ -1813,7 +1819,19 @@ def fun(x):
                     A = constraint.A
                 else:
                     A = np.atleast_2d(constraint.A)
-                return A.dot(x)
+
+                res = A.dot(x)
+                # x either has shape (N, S) or (N)
+                # (M, N) x (N, S) --> (M, S)
+                # (M, N) x (N,)   --> (M,)
+                # However, if (M, N) is a matrix then:
+                # (M, N) * (N,)   --> (M, 1), we need this to be (M,)
+                if x.ndim == 1 and res.ndim == 2:
+                    # deal with case that constraint.A is an np.matrix
+                    # see gh20041
+                    res = np.asarray(res)[:, 0]
+
+                return res
         elif isinstance(constraint, Bounds):
             def fun(x):
                 return np.asarray(x)

scipy/optimize/tests/test__differential_evolution.py

@@ -986,6 +986,26 @@ def cons_f(x):
         assert pc.num_constr == 2
         assert pc.parameter_count == 2
 
+    def test_matrix_linear_constraint(self):
+        # gh20041 supplying an np.matrix to construct a LinearConstraint caused
+        # _ConstraintWrapper to start returning constraint violations of the
+        # wrong shape.
+        with suppress_warnings() as sup:
+            sup.filter(PendingDeprecationWarning)
+            matrix = np.matrix([[1, 1, 1, 1.],
+                                [2, 2, 2, 2.]])
+        lc = LinearConstraint(matrix, 0, 1)
+        x0 = np.ones(4)
+        cw = _ConstraintWrapper(lc, x0)
+        # the shape of the constraint violation should be the same as the number
+        # of constraints applied.
+        assert cw.violation(x0).shape == (2,)
+
+        # let's try a vectorised violation call.
+        xtrial = np.arange(4 * 5).reshape(4, 5)
+        assert cw.violation(xtrial).shape == (2, 5)
+
+
     def test_L1(self):
         # Lampinen ([5]) test problem 1