Merge branch 'master' into release-automation

Author: Joao-Dionisio

CHANGELOG.md

@@ -2,9 +2,11 @@
 
 ## Unreleased
 ### Added
-- `Expr` and `GenExpr` support NumPy unary functions (`np.sin`, `np.cos`, `np.sqrt`, `np.exp`, `np.log`, `np.absolute`)
+- `Expr` and `GenExpr` support NumPy unary functions (`np.sin`, `np.cos`, `np.sqrt`, `np.exp`, `np.log`, `np.absolute`, `np.negative`)
+- `Expr` and `GenExpr` support NumPy binary functions (`np.add`, `np.subtract`, `np.multiply`, `np.divide`, `np.true_divide`, `np.power`, `np.less_equal`, `np.greater_equal`, `np.equal`)
 - Added `getBase()` and `setBase()` methods to `LP` class for getting/setting basis status
 - Added `getMemUsed()`, `getMemTotal()`, and `getMemExternEstim()` methods
+- Added `isReoptEnabled()` and raising error if not enabled upon calling `reoptSolve()`
 ### Fixed
 - Removed `Py_INCREF`/`Py_DECREF` on `Model` in `catchEvent`/`dropEvent` that caused memory leak for imbalanced usage
 - Used `getIndex()` instead of `ptr()` for sorting nonlinear expression terms to avoid nondeterministic behavior

src/pyscipopt/expr.pxi

@@ -246,7 +246,32 @@ cdef class ExprLike:
             )
 
         if method == "__call__":
-            if ufunc is np.absolute:
+            if arrays := [a for a in args if type(a) is np.ndarray]:
+                if any(a.dtype.kind not in "fiub" for a in arrays):
+                    return NotImplemented
+                # If the np.ndarray is of numeric type, all arguments are converted to
+                # MatrixExpr or MatrixGenExpr and then the ufunc is applied.
+                return ufunc(*[_ensure_matrix(a) for a in args], **kwargs)
+
+            if ufunc is np.add:
+                return args[0] + args[1]
+            elif ufunc is np.subtract:
+                return args[0] - args[1]
+            elif ufunc is np.multiply:
+                return args[0] * args[1]
+            elif ufunc in {np.divide, np.true_divide}:
+                return args[0] / args[1]
+            elif ufunc is np.power:
+                return args[0] ** args[1]
+            elif ufunc is np.negative:
+                return -args[0]
+            elif ufunc is np.less_equal:
+                return args[0] <= args[1]
+            elif ufunc is np.greater_equal:
+                return args[0] >= args[1]
+            elif ufunc is np.equal:
+                return args[0] == args[1]
+            elif ufunc is np.absolute:
                 return args[0].__abs__()
             elif ufunc is np.exp:
                 return args[0].exp()
@@ -1031,6 +1056,12 @@ cdef inline object _wrap_ufunc(object x, object ufunc):
         return res.view(MatrixGenExpr) if isinstance(res, np.ndarray) else res
     return ufunc(_to_const(x))
 
+cdef inline object _ensure_matrix(object arg):
+    if type(arg) is np.ndarray:
+        return arg.view(MatrixExpr)
+    matrix = MatrixExpr if isinstance(arg, Expr) else MatrixGenExpr
+    return np.array(arg, dtype=object).view(matrix)
+
 
 def expr_to_nodes(expr):
     '''transforms tree to an array of nodes. each node is an operator and the position of the 

src/pyscipopt/scip.pxd

@@ -1550,6 +1550,7 @@ cdef extern from "scip/scip.h":
     SCIP_RETCODE SCIPfreeReoptSolve(SCIP* scip)
     SCIP_RETCODE SCIPchgReoptObjective(SCIP* scip, SCIP_OBJSENSE objsense, SCIP_VAR** vars, SCIP_Real* coefs, int nvars)
     SCIP_RETCODE SCIPenableReoptimization(SCIP* scip, SCIP_Bool enable)
+    SCIP_Bool SCIPisReoptEnabled(SCIP* scip)
 
     BMS_BLKMEM* SCIPblkmem(SCIP* scip)
 

src/pyscipopt/scip.pxi

@@ -3973,6 +3973,17 @@ cdef class Model:
         """
         PY_SCIP_CALL(SCIPenableReoptimization(self._scip, enable))
 
+    def isReoptEnabled(self):
+        """
+        Returns whether reoptimization is enabled.
+
+        Returns
+        -------
+        bool
+
+        """
+        return SCIPisReoptEnabled(self._scip)
+
     def lpiGetIterations(self):
         """
         Get the iteration count of the last solved LP.
@@ -11950,6 +11961,10 @@ cdef class Model:
     def freeReoptSolve(self):
         """Frees all solution process data and prepares for reoptimization."""
 
+        if not SCIPisReoptEnabled(self._scip):
+            raise ValueError("freeReoptSolve requires reoptimization to be enabled. "
+                             "Call enableReoptimization() before solving.")
+
         if self.getStage() not in [SCIP_STAGE_INIT,
                                  SCIP_STAGE_PROBLEM,
                                  SCIP_STAGE_TRANSFORMED,

src/pyscipopt/scip.pyi

@@ -1454,6 +1454,7 @@ class Model:
     def isObjChangedProbing(self) -> Incomplete: ...
     def isObjIntegral(self) -> Incomplete: ...
     def isPositive(self, val: Incomplete) -> Incomplete: ...
+    def isReoptEnabled(self) -> bool: ...
     def isZero(self, value: Incomplete) -> Incomplete: ...
     def lpiGetIterations(self) -> Incomplete: ...
     def markDoNotAggrVar(self, var: Incomplete) -> Incomplete: ...

tests/test_expr.py

@@ -222,7 +222,7 @@ def test_getVal_with_GenExpr():
         m.getVal(1 / z)
 
 
-def test_unary(model):
+def test_unary_ufunc(model):
     m, x, y, z = model
 
     res = "abs(sum(0.0,prod(1.0,x)))"
@@ -276,6 +276,57 @@ def test_unary(model):
         # forbid modifying Variable/Expr/GenExpr in-place via out parameter
         np.sin(x, out=np.array([0]))
 
+    # test np.negative
+    assert str(np.negative(x)) == "Expr({Term(x): -1.0})"
+
+
+def test_binary_ufunc(model):
+    m, x, y, z = model
+
+    # test np.add
+    assert str(np.add(x, 1)) == "Expr({Term(x): 1.0, Term(): 1.0})"
+    assert str(np.add(1, x)) == "Expr({Term(x): 1.0, Term(): 1.0})"
+    a = np.array([1])
+    assert str(np.add(x, a)) == "[Expr({Term(x): 1.0, Term(): 1.0})]"
+    assert str(np.add(a, x)) == "[Expr({Term(x): 1.0, Term(): 1.0})]"
+
+    # test np.subtract
+    assert str(np.subtract(x, 1)) == "Expr({Term(x): 1.0, Term(): -1.0})"
+    assert str(np.subtract(1, x)) == "Expr({Term(x): -1.0, Term(): 1.0})"
+    assert str(np.subtract(x, a)) == "[Expr({Term(x): 1.0, Term(): -1.0})]"
+    assert str(np.subtract(a, x)) == "[Expr({Term(x): -1.0, Term(): 1.0})]"
+
+    # test np.multiply
+    a = np.array([2])
+    assert str(np.multiply(x, 2)) == "Expr({Term(x): 2.0})"
+    assert str(np.multiply(2, x)) == "Expr({Term(x): 2.0})"
+    assert str(np.multiply(x, a)) == "[Expr({Term(x): 2.0})]"
+    assert str(np.multiply(a, x)) == "[Expr({Term(x): 2.0})]"
+
+    # test np.divide
+    assert str(np.divide(x, 2)) == "Expr({Term(x): 0.5})"
+    assert str(np.divide(2, x)) == "prod(2.0,**(sum(0.0,prod(1.0,x)),-1))"
+    assert str(np.divide(x, a)) == "[Expr({Term(x): 0.5})]"
+    assert str(np.divide(a, x)) == "[prod(2.0,**(sum(0.0,prod(1.0,x)),-1))]"
+
+    # test np.power
+    assert str(np.power(x, 2)) == "Expr({Term(x, x): 1.0})"
+    assert str(np.power(2, x)) == "exp(prod(1.0,sum(0.0,prod(1.0,x)),log(2.0)))"
+    assert str(np.power(x, a)) == "[Expr({Term(x, x): 1.0})]"
+    assert str(np.power(a, x)) == "[exp(prod(1.0,sum(0.0,prod(1.0,x)),log(2.0)))]"
+
+    # test np.less_equal
+    assert str(np.less_equal(x, a)) == "[ExprCons(Expr({Term(x): 1.0}), None, 2.0)]"
+    assert str(np.less_equal(a, x)) == "[ExprCons(Expr({Term(x): 1.0}), 2.0, None)]"
+
+    # test np.equal
+    assert str(np.equal(x, a)) == "[ExprCons(Expr({Term(x): 1.0}), 2.0, 2.0)]"
+    assert str(np.equal(a, x)) == "[ExprCons(Expr({Term(x): 1.0}), 2.0, 2.0)]"
+
+    # test np.greater_equal
+    assert str(np.greater_equal(x, a)) == "[ExprCons(Expr({Term(x): 1.0}), 2.0, None)]"
+    assert str(np.greater_equal(a, x)) == "[ExprCons(Expr({Term(x): 1.0}), None, 2.0)]"
+
 
 def test_mul():
     m = Model()

tests/test_reopt.py

@@ -73,6 +73,31 @@ def test_reopt_many_variables_sparse_objective(self):
         self.assertAlmostEqual(m.getVal(vars[50]), 0.0)
         self.assertAlmostEqual(m.getVal(vars[99]), 0.0)
 
+    def test_freeReoptSolve_requires_reoptimization_enabled(self):
+        """freeReoptSolve must raise ValueError instead of crashing when
+        reoptimization was not enabled (see issue #624)."""
+        m = Model()
+        m.hideOutput()
+        self.assertFalse(m.isReoptEnabled())
+
+        x = m.addVar(name="x", vtype="I", ub=10)
+        m.addCons(x >= 3)
+        m.setObjective(x)
+        m.optimize()
+
+        with self.assertRaises(ValueError):
+            m.freeReoptSolve()
+
+        m2 = Model()
+        m2.enableReoptimization()
+        self.assertTrue(m2.isReoptEnabled())
+        m2.hideOutput()
+        y = m2.addVar(name="y", vtype="I", ub=10)
+        m2.addCons(y >= 3)
+        m2.setObjective(y)
+        m2.optimize()
+        m2.freeReoptSolve()
+
     def test_reopt_zero_objective(self):
         """Test reoptimization with zero objective (no variables, all coefficients zero)."""
         m = Model()