diff --git a/cpmpy/solvers/highs.py b/cpmpy/solvers/highs.py index 52399111b..68f864f49 100644 --- a/cpmpy/solvers/highs.py +++ b/cpmpy/solvers/highs.py @@ -42,6 +42,7 @@ import time import warnings +import cpmpy as cp import numpy as np import numpy.typing as npt @@ -236,6 +237,57 @@ def transform(self, cpm_expr): cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearisation, rewrite ~bv into 1-bv return cpm_cons + def _add_transformed(self, con): + """ + Post a single already-transformed constraint to the HiGHS model. + + Returns the HiGHS row indices created by the constraint. This is also + used by ``mus_native`` to map HiGHS IIS rows back to CPMpy constraints. + """ + if isinstance(con, Comparison): + lhs, rhs = con.args + if not is_num(rhs): + raise AssertionError(f"HiGHS: unexpected non-numeric RHS in comparison {con}, please report on our issue tracker.") + + indices, values, const = self._row_from_linexpr(lhs) + # effective rhs: lhs_vars + const rhs => lhs_vars rhs - const + bound = rhs - const + + if con.name == "<=": + lower = -self._inf + upper = bound + elif con.name == ">=": + lower = bound + upper = self._inf + elif con.name == "==": + lower = bound + upper = bound + else: + raise NotSupportedError(f"HiGHS: unexpected comparison operator after linearization: {con.name}, please report on our issue tracker.") + + row_idx = self.highs.getNumRow() + self.highs.addRow(lower, upper, len(indices), indices, values) + return [row_idx] + + elif isinstance(con, BoolVal): + if con.args[0] is False: + # add an always-infeasible row: 1 <= 0 + indices = np.empty(0, dtype=np.int32) + values = np.empty(0, dtype=np.float64) + row_idx = self.highs.getNumRow() + self.highs.addRow(1, 0, 0, indices, values) + return [row_idx] + # BoolVal(True) is a tautology; nothing to post + return [] + + elif isinstance(con, DirectConstraint): + # delegate to user-provided function with native model + con.callSolver(self, self.highs) + return [] + + else: + raise NotImplementedError(f"HiGHS: unexpected transformed constraint {con}, please report on our issue tracker.") + def add(self, cpm_expr): """ Eagerly add a constraint to the underlying solver. @@ -247,43 +299,7 @@ def add(self, cpm_expr): get_variables(cpm_expr, collect=self.user_vars) for con in self.transform(cpm_expr): - if isinstance(con, Comparison): - lhs, rhs = con.args - if not is_num(rhs): - raise AssertionError(f"HiGHS: unexpected non-numeric RHS in comparison {con}, please report on our issue tracker.") - - indices, values, const = self._row_from_linexpr(lhs) - # effective rhs: lhs_vars + const rhs => lhs_vars rhs - const - bound = rhs - const - - if con.name == "<=": - lower = -self._inf - upper = bound - elif con.name == ">=": - lower = bound - upper = self._inf - elif con.name == "==": - lower = bound - upper = bound - else: - raise NotSupportedError(f"HiGHS: unexpected comparison operator after linearization: {con.name}, please report on our issue tracker.") - - self.highs.addRow(lower, upper, len(indices), indices, values) - - elif isinstance(con, BoolVal): - if con.args[0] is False: - # add an always-infeasible row: 1 <= 0 - indices = np.empty(0, dtype=np.int32) - values = np.empty(0, dtype=np.float64) - self.highs.addRow(1, 0, 0, indices, values) - # BoolVal(True) is a tautology; nothing to post - - elif isinstance(con, DirectConstraint): - # delegate to user-provided function with native model - con.callSolver(self, self.highs) - - else: - raise NotImplementedError(f"HiGHS: unexpected transformed constraint {con}, please report on our issue tracker.") + self._add_transformed(con) return self @@ -438,3 +454,61 @@ def solve(self, time_limit=None, **kwargs): return has_sol + @classmethod + def mus_native(cls, soft, hard=[]): + """ + Compute a MUS using HiGHS' native IIS row extractor. + + A CPMpy soft constraint can transform to multiple rows, but this is not supported. + Also HiGHS' IIS support is currently LP-only and works at the level of + native rows. Even when there is a bijection of CPMpy constraints to HiGHS constraints + the returned IIS may be invalid and hence raise an + """ + import highspy + + soft_cons = toplevel_list(soft, merge_and=False) + hard_cons = toplevel_list(hard, merge_and=False) + + s = cls() + for cpm_con in s.transform(hard_cons): + raise NotSupportedError("HiGHS does not support hard constraints for MUS extraction.") + s._add_transformed(cpm_con) + + soft_rows = [] + for soft_con in soft_cons: + soft_con_tf = s.transform(soft_con) + + if len(soft_con_tf) == 0: + soft_con_rep = cp.BoolVal(True) + elif len(soft_con_tf) == 1: + soft_con_rep = soft_con_tf[0] + else: + raise NotSupportedError("HiGHS only supports MUS extraction for linear constraints.") + assumption = cp.boolvar() + additional_hard_constraint = assumption.implies(cp.all(soft_con_tf)) + for tf_con in s.transform(additional_hard_constraint): + s._add_transformed(tf_con) + soft_con_rep = assumption >= 1 + + soft_rows.append(s._add_transformed(soft_con_rep)) + + if s.solve() is not False: + raise AssertionError("MUS: model must be UNSAT") + + s.highs.setOptionValue("iis_strategy", 15) # whole LP + IIS reduction + + + status, iis = s.highs.getIis() + if status == highspy.HighsStatus.kError: + raise NotSupportedError("HiGHS: native MUS extraction failed.") + if not getattr(iis, "valid_", False): + raise NotSupportedError("HiGHS: native MUS extraction did not return a valid MUS. The infeasibility may rely on integrality.") + + iis_rows = frozenset(iis.row_index_) + native_core = [ + soft_con + for soft_con, rows in zip(soft_cons, soft_rows) + if any(row in iis_rows for row in rows) + ] + + return native_core diff --git a/cpmpy/tools/explain/mus.py b/cpmpy/tools/explain/mus.py index f583807ea..b51d3ba55 100644 --- a/cpmpy/tools/explain/mus.py +++ b/cpmpy/tools/explain/mus.py @@ -7,6 +7,7 @@ - Native MUS for given solvers: - Exact: deletion-based MUS extraction - Gurobi: IIS-based MUS extraction + - HiGHS: IIS-based MUS extraction """ import warnings import numpy as np @@ -69,7 +70,7 @@ def mus_native(soft, hard=[], solver="exact"): :param soft: soft constraints, list of expressions :param hard: hard constraints, optional, list of expressions - :param solver: which solver to use (`exact` or `gurobi`) + :param solver: which solver to use (`exact`, `gurobi` or `highs`) """ # get solver class @@ -346,4 +347,3 @@ def optimal_mus_naive(soft, hard=[], weights=None, solver="ortools", hs_solver=" Naive implementation of `optimal_mus` without assumption variables and incremental solving """ return ocus_naive(soft, hard, weights, meta_constraint=True, solver=solver, hs_solver=hs_solver) -