Add IIS-based MUS algorithm

[hbierlee]

Use computeIIS of Gurobi to compute a MUS. We still need to make an assump model to ensure that each soft constraint, represented as a group of linear constraints, is individually turned off or on in the IIS. I don't think there is a way around it.

It would be very interesting to benchmark it (I did on pigeonhole v pysat, but there is no contest there of course).

Unlike the other MUS functions, it also has time/mem limit. Let me know if either of these should be omitted to be consistent with the other approaches which do not have time/mem limits, or if I should add a separate PR to add such limits to the other algorithms. Obviously, at least a time limit would be a usefull feature, I think. [edit: remove mem_limit for simplicity]

[OrestisLomis]

I agree a time limit is useful, I have implemented a version with a time limit as well for the generalised MUS paper where I could strip the necessary changes from

[OrestisLomis]

I seem to have run into an edgecase.

import cpmpy as cp
from cpmpy.tools.explain.mus import mus, mus_iis
from cpmpy.transformations.get_variables import get_variables_model
import time

model = cp.Model.from_file("/from/local/path")

print(model)

for var in get_variables_model(model):
    print(f"{var}: {var.lb=}, {var.ub=}")

res_iis = mus_iis(model.constraints, solver="gurobi")
print(res_iis)
print(len(res_iis))

res_true = mus(res_iis, solver="exact")
print(res_true)
assert len(res_true) == len(res_iis)

This give the following:

Constraints:
    (x[0] <= 0) or (x[1] <= 0)
    (x[0] <= 1) or (x[2] <= -1)
    (x[0] <= 2) or (x[3] <= -1)
    (x[1] <= 0) or (x[2] <= 0)
    (x[0] <= -1) or (x[1] <= 1)
    (x[1] <= 2) or (x[3] <= 0)
    (x[0] <= 0) or (x[2] <= 0)
    (x[1] <= -1) or (x[2] <= 1)
    (x[2] <= 2) or (x[3] <= 1)
    (x[0] <= 1) or (x[3] <= 0)
    (x[1] <= 1) or (x[3] <= 1)
    (x[2] <= 1) or (x[3] <= 2)
    (x[0] != 0) or (x[1] == 1)
    (x[0] != 1) or (x[2] == 0)
    (x[0] != 2) or (x[3] == 0)
    (x[1] != 0) or (x[2] == 1)
    (x[1] != 1) or (x[0] == 0)
    (x[1] != 2) or (x[3] == 1)
    (x[2] != 0) or (x[0] == 1)
    (x[2] != 1) or (x[1] == 0)
    (x[2] != 2) or (x[3] == 2)
    (x[3] != 0) or (x[0] == 2)
    (x[3] != 1) or (x[1] == 2)
    (x[3] != 2) or (x[2] == 2)
Objective: None
x[0]: var.lb=0, var.ub=2
x[1]: var.lb=0, var.ub=2
x[2]: var.lb=0, var.ub=2
x[3]: var.lb=0, var.ub=2
[(x[0] <= 0) or (x[1] <= 0), (x[0] <= 0) or (x[2] <= 0), (x[1] <= -1) or (x[2] <= 1), (x[1] != 0) or (x[2] == 1), (x[2] != 0) or (x[0] == 1), (x[2] != 1) or (x[1] == 0), (x[3] != 0) or (x[0] == 2), (x[3] != 1) or (x[1] == 2), (x[3] != 2) or (x[2] == 2)]
9

[(x[3] != 2) or (x[2] == 2), (x[0] <= 0) or (x[2] <= 0), (x[2] != 0) or (x[0] == 1), (x[3] != 0) or (x[0] == 2), (x[2] != 1) or (x[1] == 0), (x[3] != 1) or (x[1] == 2), (x[0] <= 0) or (x[1] <= 0), (x[1] <= -1) or (x[2] <= 1)]
Traceback (most recent call last):
  File "/home/orestis_ubuntu/work/muser_PB/scripts/bug_iis.py", line 26, in <module>
    assert len(res_true) == len(res_iis)
AssertionError

I've been able to reproduce this bug on other instances that also have disjunctive constraints.

[hbierlee]

Thanks, there was indeed a bug! Fixed now.

Although your test does not exactly show this, it is good to have a real MUS check. My test suite is currently using:

# Helper function to verify MUS properties
def is_valid_mus(mus_result, hard=[]):
    """
    Verify that a MUS result is:
    1. Hard constraints alone are SAT (precondition for valid MUS)
    2. UNSAT (the constraints together are unsatisfiable)
    3. Minimal (removing any single constraint makes it SAT)
    """
    # Check hard constraints alone are SAT
    if hard and not cp.Model(hard).solve():
        return False, "Hard constraints alone are UNSAT (invalid problem)"

    # Check UNSAT
    if cp.Model(hard + mus_result).solve():
        return False, "MUS is not UNSAT"

    # Check minimality
    for i in range(len(mus_result)):
        subset = mus_result[:i] + mus_result[i+1:]
        if not cp.Model(hard + subset).solve():
            return False, f"MUS is not minimal: removing constraint {i} still UNSAT"

    return True, "Valid MUS"

It may be good (in a new PR) to adopt my test suite for test_tools_mus, since it includes parametrized tests for all MUS algorithms, more diverse instances, and this proper MUS check.

[hbierlee]

I removed the time limit from this PR for simplicity. I will make a separate issue to discuss MUS time limits, which may be closed as no-fix

[hbierlee]

A little description of how it now works, suppose we have:

mus_iis(hard=[c0], soft=[c1, c2 & c3])

Then we compute an IIS by introducing one assumption variable per soft constraint. Since Gurobi does not work with assumptions, we essentially compute IIS/MUS on this:

hard=[c0, s1 -> c1, s2 -> (c2 & c3)], soft=[s1, s2]

Note only the assumption variables are soft constraints, while hard constraints now include the reified soft constraints. One potential performance improvement we could do is that no assumption variable is needed for any soft constraint represented by a single linear constraint, e.g. we could compute IIS on:

hard=[c0, s2 -> (c2 & c3)], soft=[c1, s2]

[OrestisLomis]

As discussed I found some bugs while benchmarking this.

There are a couple pb instances which gurobi reports to be satisfiable contradicting other solvers and also given an error when trying to compute the IIS. Interestingly the actual error is not always the same. (This also shows gurobi does an initial check, although it doesn't seem to worsen performance on the instances where it hurts other solvers..)

There are also some XCSP instances where I get an error that a variable has not yet been added to the model. This does seem like a bug in the IIS algorithm perhaps. One such instance is added as an attachment.
DiamondFree-005.xml.lzma.zip

[hbierlee]

The "added variables" bug has been fixed, the incorrect solution from gurobi is an unrelated issue for which there is now issue #887

Indeed, gurobi will give an exception if a feasible model is given. Not sure if it's done upfront as an initial check, or if it's something that can be determined as the algorithm goes (if so, then that might explain why performance is not worsened as much).

I've tried to put in the above mentioned optimization (#880 (comment)) but it's turning out to be non-trivial. Perhaps ok to leave for now, the results were decent on the PB competition instances, namely slightly worse than deletion-based/base Exact, however much worse on XCSP. I think from that point of view, this PR should be mergeable.

From here, then we re-benchmark once we have resolved the tolerance bug, perhaps using expression trees of gurobi which may improve performance further along with other linearization improvements. If there's further interest, we could improve the IIS-based method by generalizing SAT-based MUS extensions such as e.g. model rotation, either by implementing IIS ourselves, or seeing if gurobi's callbacks are flexible enough for this purpose. The available callback evens are here, interestingly, they also have estimates on the min/max size of the IIS, meaning this approach could have different anytime qualities compared to our other MUS methods (discussed a bit in in #881)

[hbierlee]

From offline discussion

• still need to make this part of native_mus introduced in add native mus computation for exact #909
• the newly proposed non-group/group MUS approach is more complex
• however, the current approach is not appropriate for ILP without the optimization for singleton groups. This needs to be added (and benchmarked to see if it's effective)

[hbierlee]

In a weird twist, our "safe" assumption method is also not clearing that particular unit test. This may be an upstream bug. I've made a forum post. In the meantime, it should be ok to continue with the optimization.

[hbierlee]

Now added the singleton group optimization, requiring an added private function to get to the Gurobi constraints. The alternative there would be to use getConstr in-between adding soft constraints, which would then require grb_model.update, which is actually quite slow.

I think it can be reviewed, but if approved, we should probably wait with merging until we resolved the (potentially upstream) issue that I linked in the previous comment.

[tias]

good; I do have some requests wrt code readability/maintainability

[tias]

I'm confused here... you are already adding constraints on line 549?

or no wait, you only add the reified, not the soft one?

OK, this is too far away from the real logic...

you declare grb_soft on line 541 before the loop, so please build it during that loop... I mean, add s._add_transfo(c) to grb_soft instead of first adding a c and then overwriting it here 10 lines below

[hbierlee]

I've now done some effort to avoid this to make the code more readable, see latest version. There is a good reason we first just collect soft constraints without adding them to Gurobi, namely because I first add all hard constraints so I can set their IISConstrForce attribute by iterating over all Gurobi constraints, and then add the soft constraints.

Please see for an updated version which does it the old way below. I'm fine with either, some parts are simpler, others more complex. Note we have to should call model.update() as it can be expensive (but is required before changing constraint attributes). Let me know which version you prefer, the current one or the old one (in which case I can revert).

cpmpy/cpmpy/solvers/gurobi.py

Lines 535 to 593 in 8b72ba3

	hard = toplevel_list(hard)
	soft = toplevel_list(soft, merge_and=False)
	# post hard constraints to Gurobi
	s = cls(cp.Model(hard))
	# The Gurobi IIS algorithm minimizes constraints directly, unlike assumption-based solvers. However, a user-level constraint may be transformed to a group of multiple Gurobi constraints. In this case, we have to represent this group by a *single* soft constraint, otherwise the Gurobi IIS may not map to the user-level constraint MUS. We collect `tf_soft` so that `tf_soft[i]` is a single soft constraint representing `soft[i]`.
	tf_soft = []
	# note: we collect `tf_soft` here but do not add them to the solver yet, so we can first enable the `IISConstrForce` for all hard constraints
	# note: we use `s.transform`, then `_add_transformed`, to add some of the constraints to the solver. This bypasses normal creation of `user_vars`. This is safe to do since `user_vars` are not used in this algorithm, and the solver object does not leave this function.
	for con in soft:
	# manually transform the constraint so we can see whether `con` is represented by more than one constraint
	tf_cons = s.transform(con)
	if len(tf_cons) == 0:
	# this uncommon case ensures `tf_soft` maps to `soft`
	tf_soft.append(cp.BoolVal(True))
	elif len(tf_cons) == 1:
	# if `con` represented by a single transformed constraint, it can be added as-is
	tf_soft.append(tf_cons[0])
	else:
	# we represent the group of multiple Gurobi constraints `tf_cons` with a new auxiliary assumption variable `a` and adding *hard* constraint `a -> /\ C`
	assumption = cp.boolvar()
	# adding `a -> /\ C` may require re-transform due to the added implication
	s.add(assumption.implies(cp.all(tf_cons)))
	# then `a>=1` will be the single soft constraint to group `C` is in the MUS
	tf_soft.append(assumption >= 1)
	# update required to avoid `gurobipy._exception.GurobiError: GenConstr has not yet been added to the model`
	s.native_model.update()
	# force all hard constraints (including the reified soft constraints) into the IIS by enabling the `IISConstrForce`. Different constraint types have different names for this attritube.
	for con in s.native_model.getConstrs():
	con.setAttr("IISConstrForce", 1)
	for con in s.native_model.getGenConstrs():
	con.setAttr("IISGenConstrForce", 1)
	# we don't use these latter two at the time of writing, but future-proof anyway
	for con in s.native_model.getQConstrs():
	con.setAttr("IISQConstrForce", 1)
	for con in s.native_model.getSOSs():
	con.setAttr("IISSOSForce", 1)
	# now add each soft constraint or the assumption representing it
	grb_soft = [s._add_transformed(c) for c in tf_soft]
	# compute IIS (conveniently fails if original model was SAT since it will solve the model)
	try:
	s.native_model.computeIIS()
	except gp.GurobiError as e:
	if e.errno == gp.GRB.Error.IIS_NOT_INFEASIBLE:
	raise AssertionError("MUS: model must be UNSAT")
	raise
	mus = []
	for soft_i, grb_soft_i in zip(soft, grb_soft):
	# if grb_soft_i has the `IISConstr` attribute enabled, then `soft_i` is in the MUS. Again, the exact attribute depends on the constraint type.
	if any(getattr(grb_soft_i, attr, False) for attr in ("IISConstr", "IISGenConstr", "IISQConstr", "IISSOS")):
	mus.append(soft_i)
	return mus

[hbierlee]

@tias I've processed the review, see one outstanding comment. I think you can just say A or B, then approve, then I will merge.

[hbierlee]

There is now a ticket (https://support.gurobi.com/hc/en-us/requests/116323) for the confirmed IIS bug in Gurobi (wrt forcing constraints into the IIS). May be worth just waiting until it's resolved before we merge/release IIS-based MUS. Even if fixed in gurobi 13.0.2, if we this, it may be worth putting a warning in case someone has gurobi<13.0.2 (or more extremely, upping the bound for gurobi>=13.0.2, but it'd prefer the warning).