WIP: properly handle cons_{get,set}_expr and slack

Author: rschwarz

mip/highs.py

@@ -638,30 +638,61 @@ def constr_get_expr(self: "SolverHighs", constr: "mip.Constr") -> "mip.LinExpr":
         )
 
         #  - second, to get the coefficients in pre-allocated arrays.
-        matrix_start = ffi.new("int[]", 1)
-        matrix_index = ffi.new("int[]", num_nz[0])
-        matrix_value = ffi.new("double[]", num_nz[0])
-        check(
-            self._lib.Highs_getRowsByRange(
-                self._model,
-                row,
-                row,
-                num_row,
-                lower,
-                upper,
-                num_nz,
-                matrix_start,
-                matrix_index,
-                matrix_value,
+        if num_nz[0] == 0:
+            # early exit for empty expressions
+            expr = mip.xsum([])
+        else:
+            matrix_start = ffi.new("int[]", 1)
+            matrix_index = ffi.new("int[]", num_nz[0])
+            matrix_value = ffi.new("double[]", num_nz[0])
+            check(
+                self._lib.Highs_getRowsByRange(
+                    self._model,
+                    row,
+                    row,
+                    num_row,
+                    lower,
+                    upper,
+                    num_nz,
+                    matrix_start,
+                    matrix_index,
+                    matrix_value,
+                )
+            )
+            expr = mip.xsum(
+                matrix_value[i] * self.model.vars[i] for i in range(num_nz[0])
             )
-        )
 
-        return mip.xsum(matrix_value[i] * self.model.vars[i] for i in range(num_nz))
+        # Also set sense and constant
+        lhs, rhs = lower[0], upper[0]
+        if rhs < mip.INF:
+            expr -= rhs
+            if lhs > -mip.INF:
+                assert lhs == rhs
+                expr.sense = mip.EQUAL
+            else:
+                expr.sense = mip.LESS_OR_EQUAL
+        else:
+            if lhs > -mip.INF:
+                expr -= lhs
+                expr.sense = mip.GREATER_OR_EQUAL
+            else:
+                raise ValueError("Unbounded constraint?!")
+        return expr
 
     def constr_set_expr(
         self: "SolverHighs", constr: "mip.Constr", value: "mip.LinExpr"
     ) -> "mip.LinExpr":
-        raise NotImplementedError()
+        # We also have to set to 0 all coefficients of the old row, so we
+        # fetch that first.
+        coeffs = {var: 0.0 for var in constr.expr}
+
+        # Then we fetch the new coefficients and overwrite.
+        coeffs.update(value.expr.items())
+
+        # Finally, we change the coeffs in HiGHS' matrix one-by-one.
+        for var, coef in coeffs.items():
+            self._change_coef(constr.idx, var.idx, coef)
 
     def constr_get_rhs(self: "SolverHighs", idx: int) -> numbers.Real:
         # fetch both lower and upper bound
@@ -731,7 +762,19 @@ def constr_get_pi(self: "SolverHighs", constr: "mip.Constr") -> numbers.Real:
             return self._pi[constr.idx]
 
     def constr_get_slack(self: "SolverHighs", constr: "mip.Constr") -> numbers.Real:
-        raise NotImplementedError()
+        expr = constr.expr
+        activity = sum(coef * var.x for var, coef in expr.expr.items())
+        rhs = -expr.const
+        slack = rhs - activity
+        assert False
+        if expr.sense == mip.LESS_OR_EQUAL:
+            return slack
+        elif expr.sense == mip.GREATER_OR_EQUAL:
+            return -slack
+        elif expr.sense == mip.EQUAL:
+            return -abs(slack)
+        else:
+            raise ValueError(f"Invalid constraint sense: {expr.sense}")
 
     def remove_constrs(self: "SolverHighs", constrsList: List[int]):
         set_ = ffi.new("int[]", constrsList)

test/test_model.py

@@ -1368,18 +1368,29 @@ def test_solve_relaxation(solver):
     y = m.add_var("y", var_type=INTEGER)
     z = m.add_var("z", var_type=BINARY)
 
-    m.add_constr(x <= 10 * z)
-    m.add_constr(x <= 9.5)
-    m.add_constr(x + y <= 20)
+    c1 = m.add_constr(x <= 10 * z)
+    c2 = m.add_constr(x <= 9.5)
+    c3 = m.add_constr(x + y <= 20)
     m.objective = mip.maximize(4*x + y - z)
 
+    # double-check constraint expressions
+    assert c1.idx == 0
+    expr1 = c1.expr
+    assert expr1.expr == pytest.approx({x: 1.0, z: -10.0})
+    assert expr1.const == pytest.approx(0.0)
+    assert expr1.sense == mip.LESS_OR_EQUAL
+
     # first solve proper MIP
     status = m.optimize()
     assert status == OptimizationStatus.OPTIMAL
     assert x.x == pytest.approx(9.5)
     assert y.x == pytest.approx(10.0)
     assert z.x == pytest.approx(1.0)
 
+    assert c1.slack == pytest.approx(0.5)
+    assert c2.slack == pytest.approx(0.0)
+    assert c3.slack == pytest.approx(0.0)
+
     # then compare LP relaxation
     # (seems to fail for CBC?!)
     if solver == HIGHS:
@@ -1388,3 +1399,7 @@ def test_solve_relaxation(solver):
         assert x.x == pytest.approx(9.5)
         assert y.x == pytest.approx(10.5)
         assert z.x == pytest.approx(0.95)
+
+        assert c1.slack == pytest.approx(0.0)
+        assert c2.slack == pytest.approx(0.0)
+        assert c3.slack == pytest.approx(0.0)