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)

[OrestisLomis]

Some small doc comments, looks good otherwise when keeping in mind that gurobi (and IIS) may give wrong results

[OrestisLomis]

I find this wording slightly weird, which I guess is for future proofing reasons. It reads as several solvers may be used, but actually this implementation is gurobi specific.

[hbierlee]

I think I prefer the future-proof version, it is meant as an IIS based alg, not a Gurobi alg.

[OrestisLomis]

I'm adding another bugged instance here. Specifically this instance has the behaviour that the precision problems (or some other problem) is only present in the IIS computation itself. As checking the mus using gurobi does also confirm that IIS finds wrong MUSes on this instance.
I tested this without the 'dummy' assumption variables as well since it is a PB problem and hence there is no need for them. However, the problem remains.

normalized-16array_alg_ineq5.zip

[hbierlee]

Thanks, it's good to have the instance, since we suspect it is the tolerance issue, I think it's ok to just check back on it when we have fixed the tolerance issues we already know about. So I'll ping Tias and to give a final look before merging.

Note I've made the MIP->ILP changes a bit more project-wide, if this is no good, I can just revert 7c7505e

[tias]

Please put the doc changes in a separate PR which will be trivial to accept (making reviewing this one cleaner too).

An ILP modeller would never introduce BigM's on their constraints to then use a native IIS method on only auxiliary 01 vars... thats not the spirit of ILP modellling.

Let's think it more through... so in CPMpy we can have constraint groups (either because the current linearisation does not keep it as one Generalized, or because it really is a conjunction of other constraints).

So how can we support native IIS on linear/generalized constraints, while still allowing groups of constraints? Can we post-process the result, instead of pre-processing the input?

[hbierlee]

It'd be great to compute the MUS/IIS on soft constraints directly on root-level, but other than the partial optimisation I proposed, I'm not sure there is a way to make that correct? This is actually what I was hoping to do at first, I had tried out some variation of this:

# Instantiate solver (will check if solver is installed and licensed)
s = cp.SolverLookup.get(solver, cp.Model(hard))
grb_model = s.grb_model

# Force all hard constraints into the IIS (1 = force in, 0 = force out, -1 (default) = soft), as well as all `a -> c` for each assumption `a` and constraint `c` of each group/soft constraint in `soft`
grb_model.update()  # update required ; otherwise `getConstrs` can return empty
for hard_constraint in grb_model.getConstrs():
hard_constraint.IISConstrForce = 1
for hard_constraint in grb_model.getGenConstrs():  # CPMpy also posts general constraints
hard_constraint.IISGenConstrForce = 1

grb_softs = [(s_, s.add(s_)) for s_ in soft]  # hacked gurobi interface to return the solver constraints, so we have for each soft constraint the list of gurobi constraints it is represented by

import gurobipy  # Safe to import gurobipy since instantiating the Gurobi solver succeeded
try:
    grb_model.computeIIS()
except gurobipy.GurobiError as e:
    if e.errno == gurobipy.GRB.Error.IIS_NOT_INFEASIBLE:
        raise AssertionError("MUS: model must be UNSAT")
raise

# Find which constraints are in the IIS/MUS
return [soft for soft, grb_cons in grb_softs if any(getattr(grb_con, "IISConstr", False) or getattr(grb_con, "IISGenConstr", False) for grb_con in grb_cons)]

This was passing many of my tests, but then Orestis found a failing case in a benchmark instance. I've now come up with a minimal example which I think gets at the same problem. The above approach can piece together a MUS with transformed constraints from a set of user-level constraints which are not a MUS. Given two constraints c1 and c2 which (somehow) transform a & (~a) and a, leading to the transformed model a & (~a), a. At least, the above approach doesn't work here, because the MUS on the transformed constraints could return transformed constraints 2 and 3, thus a MUS of c1 and c2 (but the only real MUS is c1). It'd be great if there's a fix but I'm not sure there is?

I've been able to repro it in a test:

        valid, msg = is_valid_mus(result, hard)
>       assert valid, msg
E       AssertionError: MUS is not minimal: removing constraint 1 still UNSAT
E       assert False

tests/test_mus.py:261: AssertionError
--------------------------------------------------------------------------------------------------------------------- Captured stdout call ---------------------------------------------------------------------------------------------------------------------
HARD
[]
SOFT
[p <= 0, (p >= 1) and (p <= 0)]
VARIABLES
['p in 0..1']
MUS
[p <= 0, (p >= 1) and (p <= 0)]
============================================================================================================== 1 failed, 129 deselected in 0.06s ===============================================================================================================

I've been working on the expression tree stuff: while gurobi has a native AND General constraint, weirdly this cannot be placed in the expression tree (otherwise that'd be the way). Side-note, most of these become (or should become) indicator constraints, not big-M (perhaps internally this makes no difference).