get_constraints: Fetch slack and dual values

miplearn/features.py

@@ -43,12 +43,17 @@ class VariableFeatures:
 
 @dataclass
 class Constraint:
-    rhs: float = 0.0
+    basis_status: Optional[str] = None
+    category: Optional[Hashable] = None
+    dual_value: Optional[float] = None
+    lazy: bool = False
     lhs: Dict[str, float] = lambda: {}  # type: ignore
+    rhs: float = 0.0
+    sa_rhs_down: Optional[float] = None
+    sa_rhs_up: Optional[float] = None
     sense: str = "<"
+    slack: Optional[float] = None
     user_features: Optional[List[float]] = None
-    lazy: bool = False
-    category: Hashable = None
 
 
 @dataclass

miplearn/solvers/gurobi.py

@@ -5,6 +5,7 @@ import logging
 import re
 import sys
 from dataclasses import dataclass
+from enum import Enum
 from io import StringIO
 from random import randint
 from typing import List, Any, Dict, Optional, Hashable
@@ -68,6 +69,9 @@ class GurobiSolver(InternalSolver):
         self.lazy_cb_frequency = lazy_cb_frequency
         self._bin_vars: List["gurobipy.Var"] = []
         self._varname_to_var: Dict[str, "gurobipy.Var"] = {}
+        self._dirty = True
+        self._has_lp_solution = False
+        self._has_mip_solution = False
 
         if self.lazy_cb_frequency == 1:
             self.lazy_cb_where = [self.gp.GRB.Callback.MIPSOL]
@@ -136,9 +140,12 @@ class GurobiSolver(InternalSolver):
             var.ub = 1.0
         with _RedirectOutput(streams):
             self.model.optimize()
+            self._dirty = False
         for var in self._bin_vars:
             var.vtype = self.gp.GRB.BINARY
         log = streams[0].getvalue()
+        self._has_lp_solution = self.model.solCount > 0
+        self._has_mip_solution = False
         opt_value = None
         if not self.is_infeasible():
             opt_value = self.model.objVal
@@ -191,6 +198,7 @@ class GurobiSolver(InternalSolver):
         while True:
             with _RedirectOutput(streams):
                 self.model.optimize(cb_wrapper)
+                self._dirty = False
             if len(callback_exceptions) > 0:
                 raise callback_exceptions[0]
             total_wallclock_time += self.model.runtime
@@ -198,6 +206,8 @@ class GurobiSolver(InternalSolver):
             should_repeat = iteration_cb()
             if not should_repeat:
                 break
+        self._has_lp_solution = False
+        self._has_mip_solution = self.model.solCount > 0
 
         # Fetch results and stats
         log = streams[0].getvalue()
@@ -306,6 +316,9 @@ class GurobiSolver(InternalSolver):
             self.model.addConstr(lhs <= constr.rhs, name=name)
         else:
             self.model.addConstr(lhs >= constr.rhs, name=name)
+        self._dirty = True
+        self._has_lp_solution = False
+        self._has_mip_solution = False
 
     @overrides
     def remove_constraint(self, name: str) -> None:
@@ -341,12 +354,6 @@ class GurobiSolver(InternalSolver):
             var.lb = value
             var.ub = value
 
-    @overrides
-    def get_inequality_slacks(self) -> Dict[str, float]:
-        assert self.model is not None
-        ineqs = [c for c in self.model.getConstrs() if c.sense != "="]
-        return {c.ConstrName: c.Slack for c in ineqs}
-
     @overrides
     def relax(self) -> None:
         assert self.model is not None
@@ -414,7 +421,9 @@ class GurobiSolver(InternalSolver):
     def get_constraints(self) -> Dict[str, Constraint]:
         assert self.model is not None
         self._raise_if_callback()
+        if self._dirty:
             self.model.update()
+            self._dirty = False
         constraints: Dict[str, Constraint] = {}
         for c in self.model.getConstrs():
             constr = self._parse_gurobi_constraint(c)
@@ -422,13 +431,22 @@ class GurobiSolver(InternalSolver):
             constraints[c.constrName] = constr
         return constraints
 
-    def _parse_gurobi_constraint(self, c: Any) -> Constraint:
+    def _parse_gurobi_constraint(self, gp_constr: Any) -> Constraint:
         assert self.model is not None
-        expr = self.model.getRow(c)
+        expr = self.model.getRow(gp_constr)
         lhs: Dict[str, float] = {}
         for i in range(expr.size()):
             lhs[expr.getVar(i).varName] = expr.getCoeff(i)
-        return Constraint(rhs=c.rhs, lhs=lhs, sense=c.sense)
+        constr = Constraint(
+            rhs=gp_constr.rhs,
+            lhs=lhs,
+            sense=gp_constr.sense,
+        )
+        if self._has_lp_solution:
+            constr.dual_value = gp_constr.pi
+        if self._has_lp_solution or self._has_mip_solution:
+            constr.slack = gp_constr.slack
+        return constr
 
     @overrides
     def are_callbacks_supported(self) -> bool:

miplearn/solvers/internal.py

@@ -184,14 +184,6 @@ class InternalSolver(ABC, EnforceOverrides):
         """
         pass
 
-    @abstractmethod
-    def get_inequality_slacks(self) -> Dict[str, float]:
-        """
-        Returns a dictionary mapping constraint name to the constraint slack
-        in the current solution.
-        """
-        pass
-
     @abstractmethod
     def is_infeasible(self) -> bool:
         """

miplearn/solvers/pyomo/base.py

@@ -12,7 +12,7 @@ import numpy as np
 import pyomo
 from overrides import overrides
 from pyomo import environ as pe
-from pyomo.core import Var
+from pyomo.core import Var, Suffix
 from pyomo.core.base import _GeneralVarData
 from pyomo.core.base.constraint import ConstraintList
 from pyomo.core.expr.numeric_expr import SumExpression, MonomialTermExpression
@@ -61,6 +61,8 @@ class BasePyomoSolver(InternalSolver):
         self._varname_to_var: Dict[str, pe.Var] = {}
         self._cname_to_constr: Dict[str, pe.Constraint] = {}
         self._termination_condition: str = ""
+        self._has_lp_solution = False
+        self._has_mip_solution = False
 
         for (key, value) in params.items():
             self._pyomo_solver.options[key] = value
@@ -76,10 +78,14 @@ class BasePyomoSolver(InternalSolver):
             streams += [sys.stdout]
         with _RedirectOutput(streams):
             results = self._pyomo_solver.solve(tee=True)
+        self._termination_condition = results["Solver"][0]["Termination condition"]
         self._restore_integrality()
         opt_value = None
+        self._has_lp_solution = False
+        self._has_mip_solution = False
         if not self.is_infeasible():
             opt_value = results["Problem"][0]["Lower bound"]
+            self._has_lp_solution = True
         return {
             "LP value": opt_value,
             "LP log": streams[0].getvalue(),
@@ -122,7 +128,10 @@ class BasePyomoSolver(InternalSolver):
         ws_value = self._extract_warm_start_value(log)
         self._termination_condition = results["Solver"][0]["Termination condition"]
         lb, ub = None, None
+        self._has_mip_solution = False
+        self._has_lp_solution = False
         if not self.is_infeasible():
+            self._has_mip_solution = True
             lb = results["Problem"][0]["Lower bound"]
             ub = results["Problem"][0]["Upper bound"]
         stats: MIPSolveStats = {
@@ -185,6 +194,7 @@ class BasePyomoSolver(InternalSolver):
         self.instance = instance
         self.model = model
         self.model.extra_constraints = ConstraintList()
+        self.model.dual = Suffix(direction=Suffix.IMPORT)
         self._pyomo_solver.set_instance(model)
         self._update_obj()
         self._update_vars()
@@ -256,6 +266,9 @@ class BasePyomoSolver(InternalSolver):
             self._cname_to_constr[name] = cl
         else:
             self._pyomo_solver.add_constraint(constr)
+        self._termination_condition = ""
+        self._has_lp_solution = False
+        self._has_mip_solution = False
 
     @overrides
     def remove_constraint(self, name: str) -> None:
@@ -321,22 +334,14 @@ class BasePyomoSolver(InternalSolver):
             var.domain = pyomo.core.base.set_types.Reals
             self._pyomo_solver.update_var(var)
 
-    @overrides
-    def get_inequality_slacks(self) -> Dict[str, float]:
-        result: Dict[str, float] = {}
-        for (cname, cobj) in self._cname_to_constr.items():
-            if cobj.equality:
-                continue
-            result[cname] = cobj.slack()
-        return result
-
     @overrides
     def is_infeasible(self) -> bool:
         return self._termination_condition == TerminationCondition.infeasible
 
     @overrides
     def get_dual(self, cid: str) -> float:
-        raise NotImplementedError()
+        constr = self._cname_to_constr[cid]
+        return self._pyomo_solver.dual[constr]
 
     @overrides
     def get_sense(self) -> str:
@@ -376,45 +381,55 @@ class BasePyomoSolver(InternalSolver):
 
         return constraints
 
-    @staticmethod
+    def _parse_pyomo_constraint(
-    def _parse_pyomo_constraint(c: pyomo.core.Constraint) -> Constraint:
+        self,
+        pyomo_constr: pyomo.core.Constraint,
+    ) -> Constraint:
+        constr = Constraint()
+
         # Extract RHS and sense
-        has_ub = c.has_ub()
+        has_ub = pyomo_constr.has_ub()
-        has_lb = c.has_lb()
+        has_lb = pyomo_constr.has_lb()
         assert (
-            (not has_lb) or (not has_ub) or c.upper() == c.lower()
+            (not has_lb) or (not has_ub) or pyomo_constr.upper() == pyomo_constr.lower()
         ), "range constraints not supported"
         if has_lb:
-            sense = ">"
+            constr.sense = ">"
-            rhs = c.lower()
+            constr.rhs = pyomo_constr.lower()
         elif has_ub:
-            sense = "<"
+            constr.sense = "<"
-            rhs = c.upper()
+            constr.rhs = pyomo_constr.upper()
         else:
-            sense = "="
+            constr.sense = "="
-            rhs = c.upper()
+            constr.rhs = pyomo_constr.upper()
 
         # Extract LHS
         lhs = {}
-        if isinstance(c.body, SumExpression):
+        if isinstance(pyomo_constr.body, SumExpression):
-            for term in c.body._args_:
+            for term in pyomo_constr.body._args_:
                 if isinstance(term, MonomialTermExpression):
                     lhs[term._args_[1].name] = term._args_[0]
                 elif isinstance(term, _GeneralVarData):
                     lhs[term.name] = 1.0
                 else:
                     raise Exception(f"Unknown term type: {term.__class__.__name__}")
-        elif isinstance(c.body, _GeneralVarData):
+        elif isinstance(pyomo_constr.body, _GeneralVarData):
-            lhs[c.body.name] = 1.0
+            lhs[pyomo_constr.body.name] = 1.0
         else:
-            raise Exception(f"Unknown expression type: {c.body.__class__.__name__}")
+            raise Exception(
+                f"Unknown expression type: {pyomo_constr.body.__class__.__name__}"
+            )
+        constr.lhs = lhs
+
+        # Extract solution attributes
+        if self._has_lp_solution:
+            constr.dual_value = self.model.dual[pyomo_constr]
+
+        if self._has_mip_solution or self._has_lp_solution:
+            constr.slack = pyomo_constr.slack()
 
         # Build constraint
-        return Constraint(
+        return constr
-            lhs=lhs,
-            rhs=rhs,
-            sense=sense,
-        )
 
     @overrides
     def are_callbacks_supported(self) -> bool:

miplearn/solvers/tests/__init__.py

@@ -12,6 +12,14 @@ from miplearn.solvers.internal import InternalSolver
 # This file is in the main source folder, so that it can be called from Julia.
 
 
+def _round_constraints(constraints):
+    for (cname, c) in constraints.items():
+        for attr in ["slack", "dual_value"]:
+            if getattr(c, attr) is not None:
+                setattr(c, attr, round(getattr(c, attr), 6))
+    return constraints
+
+
 def run_internal_solver_tests(solver: InternalSolver) -> None:
     run_basic_usage_tests(solver.clone())
     run_warm_start_tests(solver.clone())
@@ -22,20 +30,38 @@ def run_internal_solver_tests(solver: InternalSolver) -> None:
 
 
 def run_basic_usage_tests(solver: InternalSolver) -> None:
+    # Create and set instance
     instance = solver.build_test_instance_knapsack()
     model = instance.to_model()
     solver.set_instance(instance, model)
+
+    # Fetch variables (after-load)
     assert_equals(
         solver.get_variable_names(),
         ["x[0]", "x[1]", "x[2]", "x[3]"],
     )
 
+    # Fetch constraints (after-load)
+    assert_equals(
+        _round_constraints(solver.get_constraints()),
+        {
+            "eq_capacity": Constraint(
+                lazy=False,
+                lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0},
+                rhs=67.0,
+                sense="<",
+            )
+        },
+    )
+
+    # Solve linear programming relaxation
     lp_stats = solver.solve_lp()
     assert not solver.is_infeasible()
     assert lp_stats["LP value"] is not None
     assert_equals(round(lp_stats["LP value"], 3), 1287.923)
     assert len(lp_stats["LP log"]) > 100
 
+    # Fetch variables (after-lp)
     solution = solver.get_solution()
     assert solution is not None
     assert solution["x[0]"] is not None
@@ -47,6 +73,22 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     assert_equals(round(solution["x[2]"], 3), 1.000)
     assert_equals(round(solution["x[3]"], 3), 0.000)
 
+    # Fetch constraints (after-lp)
+    assert_equals(
+        _round_constraints(solver.get_constraints()),
+        {
+            "eq_capacity": Constraint(
+                lazy=False,
+                lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0},
+                rhs=67.0,
+                sense="<",
+                slack=0.0,
+                dual_value=13.538462,
+            )
+        },
+    )
+
+    # Solve MIP
     mip_stats = solver.solve(
         tee=True,
         iteration_cb=None,
@@ -60,6 +102,7 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     assert_equals(mip_stats["Sense"], "max")
     assert isinstance(mip_stats["Wallclock time"], float)
 
+    # Fetch variables (after-mip)
     solution = solver.get_solution()
     assert solution is not None
     assert solution["x[0]"] is not None
@@ -71,13 +114,15 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     assert_equals(solution["x[2]"], 1.0)
     assert_equals(solution["x[3]"], 1.0)
 
+    # Fetch constraints (after-mip)
     assert_equals(
-        solver.get_constraints(),
+        _round_constraints(solver.get_constraints()),
         {
             "eq_capacity": Constraint(
                 lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0},
                 rhs=67.0,
                 sense="<",
+                slack=6.0,
             ),
         },
     )
@@ -86,10 +131,11 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     cut = Constraint(lhs={"x[0]": 1.0}, sense="<", rhs=0.0)
     assert not solver.is_constraint_satisfied(cut)
 
-    # Add new constraint and verify that it is listed
+    # Add new constraint and verify that it is listed. Modifying the model should
+    # also clear the current solution.
     solver.add_constraint(cut, "cut")
     assert_equals(
-        solver.get_constraints(),
+        _round_constraints(solver.get_constraints()),
         {
             "eq_capacity": Constraint(
                 lhs={"x[0]": 23.0, "x[1]": 26.0, "x[2]": 20.0, "x[3]": 18.0},
@@ -104,17 +150,11 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
         },
     )
 
-    # New constraint should affect the solution
+    # Re-solve MIP and verify that constraint affects the solution
     stats = solver.solve()
     assert_equals(stats["Lower bound"], 1030.0)
     assert solver.is_constraint_satisfied(cut)
 
-    # Verify slacks
-    assert_equals(
-        solver.get_inequality_slacks(),
-        {"cut": 0.0, "eq_capacity": 3.0},
-    )
-
     # Remove the new constraint
     solver.remove_constraint("cut")
 
@@ -122,9 +162,6 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     stats = solver.solve()
     assert_equals(stats["Lower bound"], 1183.0)
 
-    # Constraint should not be satisfied by current solution
-    assert not solver.is_constraint_satisfied(cut)
-
 
 def run_warm_start_tests(solver: InternalSolver) -> None:
     instance = solver.build_test_instance_knapsack()