Extract LHS as a sparse matrix

miplearn/features/extractor.py

@@ -39,7 +39,7 @@ class FeaturesExtractor:
         sample.put_array("static_var_types", variables.types)
         sample.put_array("static_var_upper_bounds", variables.upper_bounds)
         sample.put_array("static_constr_names", constraints.names)
-        # sample.put("static_constr_lhs", constraints.lhs)
+        sample.put_sparse("static_constr_lhs", constraints.lhs)
         sample.put_array("static_constr_rhs", constraints.rhs)
         sample.put_array("static_constr_senses", constraints.senses)
 

miplearn/features/sample.py

@@ -76,7 +76,9 @@ class Sample(ABC):
         assert value.dtype.kind in "biufS", f"Unsupported dtype: {value.dtype}"
 
     def _assert_is_sparse(self, value: Any) -> None:
-        assert isinstance(value, coo_matrix)
+        assert isinstance(
+            value, coo_matrix
+        ), f"coo_matrix expected; found: {value.__class__}"
         self._assert_is_array(value.data)
 
 

miplearn/solvers/gurobi.py

@@ -10,6 +10,7 @@ from typing import List, Any, Dict, Optional, TYPE_CHECKING
 
 import numpy as np
 from overrides import overrides
+from scipy.sparse import coo_matrix, lil_matrix
 
 from miplearn.instance.base import Instance
 from miplearn.solvers import _RedirectOutput
@@ -99,17 +100,19 @@ class GurobiSolver(InternalSolver):
         assert cf.lhs is not None
         assert cf.rhs is not None
         assert self.model is not None
+        lhs = cf.lhs.tocsr()
         for i in range(len(cf.names)):
             sense = cf.senses[i]
-            lhs = self.gp.quicksum(
-                self._varname_to_var[varname] * coeff for (varname, coeff) in cf.lhs[i]
+            row = lhs[i, :]
+            row_expr = self.gp.quicksum(
+                self._gp_vars[row.indices[j]] * row.data[j] for j in range(row.getnnz())
             )
             if sense == b"=":
-                self.model.addConstr(lhs == cf.rhs[i], name=cf.names[i])
+                self.model.addConstr(row_expr == cf.rhs[i], name=cf.names[i])
             elif sense == b"<":
-                self.model.addConstr(lhs <= cf.rhs[i], name=cf.names[i])
+                self.model.addConstr(row_expr <= cf.rhs[i], name=cf.names[i])
             elif sense == b">":
-                self.model.addConstr(lhs >= cf.rhs[i], name=cf.names[i])
+                self.model.addConstr(row_expr >= cf.rhs[i], name=cf.names[i])
             else:
                 raise Exception(f"Unknown sense: {sense}")
         self.model.update()
@@ -133,18 +136,18 @@ class GurobiSolver(InternalSolver):
         assert cf.rhs is not None
         assert self.model is not None
         result = []
+        x = np.array(self.model.getAttr("x", self.model.getVars()))
+        lhs = cf.lhs.tocsr() * x
         for i in range(len(cf.names)):
             sense = cf.senses[i]
-            lhs = sum(
-                self._varname_to_var[varname].x * coeff
-                for (varname, coeff) in cf.lhs[i]
-            )
-            if sense == "<":
-                result.append(lhs <= cf.rhs[i] + tol)
-            elif sense == ">":
-                result.append(lhs >= cf.rhs[i] - tol)
+            if sense == b"<":
+                result.append(lhs[i] <= cf.rhs[i] + tol)
+            elif sense == b">":
+                result.append(lhs[i] >= cf.rhs[i] - tol)
+            elif sense == b"<":
+                result.append(abs(cf.rhs[i] - lhs[i]) <= tol)
             else:
-                result.append(abs(cf.rhs[i] - lhs) <= tol)
+                raise Exception(f"unknown sense: {sense}")
         return result
 
     @overrides
@@ -214,7 +217,7 @@ class GurobiSolver(InternalSolver):
 
         gp_constrs = model.getConstrs()
         constr_names = np.array(model.getAttr("constrName", gp_constrs), dtype="S")
-        lhs: Optional[List] = None
+        lhs: Optional[coo_matrix] = None
         rhs, senses, slacks, basis_status = None, None, None, None
         dual_value, basis_status, sa_rhs_up, sa_rhs_down = None, None, None, None
 
@@ -222,13 +225,14 @@ class GurobiSolver(InternalSolver):
             rhs = np.array(model.getAttr("rhs", gp_constrs), dtype=float)
             senses = np.array(model.getAttr("sense", gp_constrs), dtype="S")
             if with_lhs:
-                lhs = [None for _ in gp_constrs]
+                nrows = len(gp_constrs)
+                ncols = len(self._var_names)
+                tmp = lil_matrix((nrows, ncols), dtype=float)
                 for (i, gp_constr) in enumerate(gp_constrs):
                     expr = model.getRow(gp_constr)
-                    lhs[i] = [
-                        (self._var_names[expr.getVar(j).index], expr.getCoeff(j))
-                        for j in range(expr.size())
-                    ]
+                    for j in range(expr.size()):
+                        tmp[i, j] = expr.getCoeff(j)
+                lhs = tmp.tocoo()
 
         if self._has_lp_solution:
             dual_value = np.array(model.getAttr("pi", gp_constrs), dtype=float)

miplearn/solvers/internal.py

@@ -5,9 +5,10 @@
 import logging
 from abc import ABC, abstractmethod
 from dataclasses import dataclass
-from typing import Any, Optional, List, Tuple, TYPE_CHECKING
+from typing import Any, Optional, List, TYPE_CHECKING
 
 import numpy as np
+from scipy.sparse import coo_matrix
 
 from miplearn.instance.base import Instance
 from miplearn.types import (
@@ -71,7 +72,7 @@ class Constraints:
     basis_status: Optional[np.ndarray] = None
     dual_values: Optional[np.ndarray] = None
     lazy: Optional[np.ndarray] = None
-    lhs: Optional[List[List[Tuple[bytes, float]]]] = None
+    lhs: Optional[coo_matrix] = None
     names: Optional[np.ndarray] = None
     rhs: Optional[np.ndarray] = None
     sa_rhs_down: Optional[np.ndarray] = None
@@ -104,7 +105,7 @@ class Constraints:
             ),
             names=(None if self.names is None else self.names[selected]),
             lazy=(None if self.lazy is None else self.lazy[selected]),
-            lhs=self._filter(self.lhs, selected),
+            lhs=(None if self.lhs is None else self.lhs.tocsr()[selected].tocoo()),
             rhs=(None if self.rhs is None else self.rhs[selected]),
             sa_rhs_down=(
                 None if self.sa_rhs_down is None else self.sa_rhs_down[selected]
@@ -114,15 +115,6 @@ class Constraints:
             slacks=(None if self.slacks is None else self.slacks[selected]),
         )
 
-    def _filter(
-        self,
-        obj: Optional[List],
-        selected: List[bool],
-    ) -> Optional[List]:
-        if obj is None:
-            return None
-        return [obj[i] for (i, selected_i) in enumerate(selected) if selected_i]
-
 
 class InternalSolver(ABC):
     """

miplearn/solvers/pyomo/base.py

@@ -6,7 +6,7 @@ import logging
 import re
 import sys
 from io import StringIO
-from typing import Any, List, Dict, Optional, Tuple
+from typing import Any, List, Dict, Optional
 
 import numpy as np
 import pyomo
@@ -18,6 +18,7 @@ from pyomo.core.base.constraint import ConstraintList
 from pyomo.core.expr.numeric_expr import SumExpression, MonomialTermExpression
 from pyomo.opt import TerminationCondition
 from pyomo.opt.base.solvers import SolverFactory
+from scipy.sparse import coo_matrix
 
 from miplearn.instance.base import Instance
 from miplearn.solvers import _RedirectOutput, _none_if_empty
@@ -58,6 +59,7 @@ class BasePyomoSolver(InternalSolver):
         self._pyomo_solver: SolverFactory = solver_factory
         self._obj_sense: str = "min"
         self._varname_to_var: Dict[bytes, pe.Var] = {}
+        self._varname_to_idx: Dict[str, int] = {}
         self._cname_to_constr: Dict[str, pe.Constraint] = {}
         self._termination_condition: str = ""
         self._has_lp_solution = False
@@ -84,23 +86,24 @@ class BasePyomoSolver(InternalSolver):
         assert cf.lhs is not None
         assert cf.rhs is not None
         assert self.model is not None
-        for (i, name) in enumerate(cf.names):
-            lhs = 0.0
-            for (varname, coeff) in cf.lhs[i]:
-                var = self._varname_to_var[varname]
-                lhs += var * coeff
+        lhs = cf.lhs.tocsr()
+        for i in range(len(cf.names)):
+            row = lhs[i, :]
+            lhsi = 0.0
+            for j in range(row.getnnz()):
+                lhsi += self._all_vars[row.indices[j]] * row.data[j]
             if cf.senses[i] == b"=":
-                expr = lhs == cf.rhs[i]
+                expr = lhsi == cf.rhs[i]
             elif cf.senses[i] == b"<":
-                expr = lhs <= cf.rhs[i]
+                expr = lhsi <= cf.rhs[i]
             elif cf.senses[i] == b">":
-                expr = lhs >= cf.rhs[i]
+                expr = lhsi >= cf.rhs[i]
             else:
                 raise Exception(f"Unknown sense: {cf.senses[i]}")
-            cl = pe.Constraint(expr=expr, name=name)
-            self.model.add_component(name.decode(), cl)
+            cl = pe.Constraint(expr=expr, name=cf.names[i])
+            self.model.add_component(cf.names[i].decode(), cl)
             self._pyomo_solver.add_constraint(cl)
-            self._cname_to_constr[name] = cl
+            self._cname_to_constr[cf.names[i]] = cl
         self._termination_condition = ""
         self._has_lp_solution = False
         self._has_mip_solution = False
@@ -119,18 +122,18 @@ class BasePyomoSolver(InternalSolver):
         assert cf.lhs is not None
         assert cf.rhs is not None
         assert cf.senses is not None
+        x = [v.value for v in self._all_vars]
+        lhs = cf.lhs.tocsr() * x
         result = []
-        for (i, name) in enumerate(cf.names):
-            lhs = 0.0
-            for (varname, coeff) in cf.lhs[i]:
-                var = self._varname_to_var[varname]
-                lhs += var.value * coeff
-            if cf.senses[i] == "<":
-                result.append(lhs <= cf.rhs[i] + tol)
-            elif cf.senses[i] == ">":
-                result.append(lhs >= cf.rhs[i] - tol)
+        for i in range(len(lhs)):
+            if cf.senses[i] == b"<":
+                result.append(lhs[i] <= cf.rhs[i] + tol)
+            elif cf.senses[i] == b">":
+                result.append(lhs[i] >= cf.rhs[i] - tol)
+            elif cf.senses[i] == b"=":
+                result.append(abs(cf.rhs[i] - lhs[i]) < tol)
             else:
-                result.append(abs(cf.rhs[i] - lhs) < tol)
+                raise Exception(f"unknown sense: {cf.senses[i]}")
         return result
 
     @overrides
@@ -163,15 +166,17 @@ class BasePyomoSolver(InternalSolver):
     ) -> Constraints:
         model = self.model
         assert model is not None
-
         names: List[str] = []
         rhs: List[float] = []
-        lhs: List[List[Tuple[bytes, float]]] = []
         senses: List[str] = []
         dual_values: List[float] = []
         slacks: List[float] = []
+        lhs_row: List[int] = []
+        lhs_col: List[int] = []
+        lhs_data: List[float] = []
+        lhs: Optional[coo_matrix] = None
 
-        def _parse_constraint(c: pe.Constraint) -> None:
+        def _parse_constraint(c: pe.Constraint, row: int) -> None:
             assert model is not None
             if with_static:
                 # Extract RHS and sense
@@ -192,30 +197,31 @@ class BasePyomoSolver(InternalSolver):
 
                 if with_lhs:
                     # Extract LHS
-                    lhsc = []
                     expr = c.body
                     if isinstance(expr, SumExpression):
                         for term in expr._args_:
                             if isinstance(term, MonomialTermExpression):
-                                lhsc.append(
-                                    (
-                                        term._args_[1].name.encode(),
-                                        float(term._args_[0]),
-                                    )
+                                lhs_row.append(row)
+                                lhs_col.append(
+                                    self._varname_to_idx[term._args_[1].name]
                                 )
+                                lhs_data.append(float(term._args_[0]))
                             elif isinstance(term, _GeneralVarData):
-                                lhsc.append((term.name.encode(), 1.0))
+                                lhs_row.append(row)
+                                lhs_col.append(self._varname_to_idx[term.name])
+                                lhs_data.append(1.0)
                             else:
                                 raise Exception(
                                     f"Unknown term type: {term.__class__.__name__}"
                                 )
                     elif isinstance(expr, _GeneralVarData):
-                        lhsc.append((expr.name.encode(), 1.0))
+                        lhs_row.append(row)
+                        lhs_col.append(self._varname_to_idx[expr.name])
+                        lhs_data.append(1.0)
                     else:
                         raise Exception(
                             f"Unknown expression type: {expr.__class__.__name__}"
                         )
-                    lhs.append(lhsc)
 
             # Extract dual values
             if self._has_lp_solution:
@@ -225,20 +231,26 @@ class BasePyomoSolver(InternalSolver):
             if self._has_mip_solution or self._has_lp_solution:
                 slacks.append(model.slack[c])
 
-        for constr in model.component_objects(pyomo.core.Constraint):
+        curr_row = 0
+        for (i, constr) in enumerate(model.component_objects(pyomo.core.Constraint)):
             if isinstance(constr, pe.ConstraintList):
                 for idx in constr:
-                    names.append(f"{constr.name}[{idx}]")
-                    _parse_constraint(constr[idx])
+                    names.append(constr[idx].name)
+                    _parse_constraint(constr[idx], curr_row)
+                    curr_row += 1
             else:
                 names.append(constr.name)
-                _parse_constraint(constr)
+                _parse_constraint(constr, curr_row)
+                curr_row += 1
+
+        if len(lhs_data) > 0:
+            lhs = coo_matrix((lhs_data, (lhs_row, lhs_col))).tocoo()
 
         return Constraints(
             names=_none_if_empty(np.array(names, dtype="S")),
             rhs=_none_if_empty(np.array(rhs, dtype=float)),
             senses=_none_if_empty(np.array(senses, dtype="S")),
-            lhs=_none_if_empty(lhs),
+            lhs=lhs,
             slacks=_none_if_empty(np.array(slacks, dtype=float)),
             dual_values=_none_if_empty(np.array(dual_values, dtype=float)),
         )
@@ -264,7 +276,7 @@ class BasePyomoSolver(InternalSolver):
             for index in var:
                 if var[index].fixed:
                     continue
-                solution[f"{var}[{index}]".encode()] = var[index].value
+                solution[var[index].name.encode()] = var[index].value
         return solution
 
     @overrides
@@ -289,9 +301,9 @@ class BasePyomoSolver(InternalSolver):
 
                 # Variable name
                 if idx is None:
-                    names.append(str(var))
+                    names.append(var.name)
                 else:
-                    names.append(f"{var}[{idx}]")
+                    names.append(var[idx].name)
 
                 if with_static:
                     # Variable type
@@ -556,12 +568,14 @@ class BasePyomoSolver(InternalSolver):
         self._all_vars = []
         self._bin_vars = []
         self._varname_to_var = {}
+        self._varname_to_idx = {}
         for var in self.model.component_objects(Var):
             for idx in var:
-                varname = f"{var.name}[{idx}]".encode()
-                if idx is None:
-                    varname = var.name.encode()
-                self._varname_to_var[varname] = var[idx]
+                varname = var.name
+                if idx is not None:
+                    varname = var[idx].name
+                self._varname_to_var[varname.encode()] = var[idx]
+                self._varname_to_idx[varname] = len(self._all_vars)
                 self._all_vars += [var[idx]]
                 if var[idx].domain == pyomo.core.base.set_types.Binary:
                     self._bin_vars += [var[idx]]
@@ -575,7 +589,7 @@ class BasePyomoSolver(InternalSolver):
         for constr in self.model.component_objects(pyomo.core.Constraint):
             if isinstance(constr, pe.ConstraintList):
                 for idx in constr:
-                    self._cname_to_constr[f"{constr.name}[{idx}]"] = constr[idx]
+                    self._cname_to_constr[constr[idx].name] = constr[idx]
             else:
                 self._cname_to_constr[constr.name] = constr
 

miplearn/solvers/tests/__init__.py

@@ -5,6 +5,7 @@
 from typing import Any, List
 
 import numpy as np
+from scipy.sparse import coo_matrix
 
 from miplearn.solvers.internal import InternalSolver, Variables, Constraints
 
@@ -55,15 +56,7 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
         Constraints(
             names=np.array(["eq_capacity"], dtype="S"),
             rhs=np.array([0.0]),
-            lhs=[
-                [
-                    (b"x[0]", 23.0),
-                    (b"x[1]", 26.0),
-                    (b"x[2]", 20.0),
-                    (b"x[3]", 18.0),
-                    (b"z", -1.0),
-                ],
-            ],
+            lhs=coo_matrix([[23.0, 26.0, 20.0, 18.0, -1.0]]),
             senses=np.array(["="], dtype="S"),
         ),
     )
@@ -162,7 +155,7 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
     # Build new constraint and verify that it is violated
     cf = Constraints(
         names=np.array(["cut"], dtype="S"),
-        lhs=[[(b"x[0]", 1.0)]],
+        lhs=coo_matrix([[1.0, 0.0, 0.0, 0.0, 0.0]]),
         rhs=np.array([0.0]),
         senses=np.array(["<"], dtype="S"),
     )
@@ -177,18 +170,12 @@ def run_basic_usage_tests(solver: InternalSolver) -> None:
             Constraints(
                 names=np.array(["eq_capacity", "cut"], dtype="S"),
                 rhs=np.array([0.0, 0.0]),
-                lhs=[
+                lhs=coo_matrix(
                     [
-                        (b"x[0]", 23.0),
-                        (b"x[1]", 26.0),
-                        (b"x[2]", 20.0),
-                        (b"x[3]", 18.0),
-                        (b"z", -1.0),
-                    ],
-                    [
-                        (b"x[0]", 1.0),
-                    ],
-                ],
+                        [23.0, 26.0, 20.0, 18.0, -1.0],
+                        [1.0, 0.0, 0.0, 0.0, 0.0],
+                    ]
+                ),
                 senses=np.array(["=", "<"], dtype="S"),
             ),
         ),
@@ -275,6 +262,8 @@ def _equals_preprocess(obj: Any) -> Any:
             return np.round(obj, decimals=6).tolist()
         else:
             return obj.tolist()
+    elif isinstance(obj, coo_matrix):
+        return obj.todense().tolist()
     elif isinstance(obj, (int, str, bool, np.bool_, np.bytes_, bytes, bytearray)):
         return obj
     elif isinstance(obj, float):

tests/features/test_extractor.py

@@ -9,6 +9,7 @@ from typing import Any
 
 import numpy as np
 import gurobipy as gp
+from scipy.sparse import coo_matrix
 
 from miplearn.features.extractor import FeaturesExtractor
 from miplearn.features.sample import Hdf5Sample, MemorySample
@@ -76,18 +77,10 @@ def test_knapsack() -> None:
         sample.get_array("static_constr_names"),
         np.array(["eq_capacity"], dtype="S"),
     )
-    # assert_equals(
-    #     sample.get_array("static_constr_lhs"),
-    #     [
-    #         [
-    #             ("x[0]", 23.0),
-    #             ("x[1]", 26.0),
-    #             ("x[2]", 20.0),
-    #             ("x[3]", 18.0),
-    #             ("z", -1.0),
-    #         ],
-    #     ],
-    # )
+    assert_equals(
+        sample.get_sparse("static_constr_lhs"),
+        [[23.0, 26.0, 20.0, 18.0, -1.0]],
+    )
     assert_equals(
         sample.get_array("static_constr_rhs"),
         np.array([0.0]),
@@ -321,20 +314,13 @@ def test_constraint_getindex() -> None:
         names=np.array(["c1", "c2", "c3"], dtype="S"),
         rhs=np.array([1.0, 2.0, 3.0]),
         senses=np.array(["=", "<", ">"], dtype="S"),
-        lhs=[
-            [
-                (b"x1", 1.0),
-                (b"x2", 1.0),
-            ],
+        lhs=coo_matrix(
             [
-                (b"x2", 2.0),
-                (b"x3", 2.0),
-            ],
-            [
-                (b"x3", 3.0),
-                (b"x4", 3.0),
-            ],
-        ],
+                [1, 2, 3],
+                [4, 5, 6],
+                [7, 8, 9],
+            ]
+        ),
     )
     assert_equals(
         cf[[True, False, True]],
@@ -342,16 +328,12 @@ def test_constraint_getindex() -> None:
             names=np.array(["c1", "c3"], dtype="S"),
             rhs=np.array([1.0, 3.0]),
             senses=np.array(["=", ">"], dtype="S"),
-            lhs=[
-                [
-                    (b"x1", 1.0),
-                    (b"x2", 1.0),
-                ],
+            lhs=coo_matrix(
                 [
-                    (b"x3", 3.0),
-                    (b"x4", 3.0),
-                ],
-            ],
+                    [1, 2, 3],
+                    [7, 8, 9],
+                ]
+            ),
         ),
     )