diff --git a/miplearn/features/extractor.py b/miplearn/features/extractor.py index 872767d..7cf0c2f 100644 --- a/miplearn/features/extractor.py +++ b/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) diff --git a/miplearn/features/sample.py b/miplearn/features/sample.py index 71d76e3..50d9f72 100644 --- a/miplearn/features/sample.py +++ b/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) diff --git a/miplearn/solvers/gurobi.py b/miplearn/solvers/gurobi.py index 536750e..d49906a 100644 --- a/miplearn/solvers/gurobi.py +++ b/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) diff --git a/miplearn/solvers/internal.py b/miplearn/solvers/internal.py index 1865159..9688cb0 100644 --- a/miplearn/solvers/internal.py +++ b/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): """ diff --git a/miplearn/solvers/pyomo/base.py b/miplearn/solvers/pyomo/base.py index 3e306a5..5292eb0 100644 --- a/miplearn/solvers/pyomo/base.py +++ b/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 diff --git a/miplearn/solvers/tests/__init__.py b/miplearn/solvers/tests/__init__.py index f4abe83..3bc74d3 100644 --- a/miplearn/solvers/tests/__init__.py +++ b/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): diff --git a/tests/features/test_extractor.py b/tests/features/test_extractor.py index 858d976..26cd7a0 100644 --- a/tests/features/test_extractor.py +++ b/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], + ] + ), ), )