adding QuickXplain (#416)

* add quickxplain

* doc updates

* typo

* remove order argument, rely on ordering of list

* fix import

* fix import

* add quickxplain to tools and refactor algo

* whoops

* naive version of quickxplain

* update tests

* update example

* update tests

* update docs

cpmpy/tools/mus.py

@@ -6,11 +6,12 @@
 
 """
 import numpy as np
-from cpmpy import *
+import cpmpy as cp
 from cpmpy.expressions.variables import NDVarArray
 from cpmpy.transformations.get_variables import get_variables
 from cpmpy.transformations.normalize import toplevel_list
 
+
 def mus(soft, hard=[], solver="ortools"):
     """
         A CP deletion-based MUS algorithm using assumption variables
@@ -36,28 +37,28 @@ def mus(soft, hard=[], solver="ortools"):
     # order so that constraints with many variables are tried and removed first
     candidates = sorted(soft, key=lambda c: -len(get_variables(c)))
 
-    assump = boolvar(shape=len(soft), name="assump")
+    assump = cp.boolvar(shape=len(soft), name="assump")
     if len(soft) == 1:
         assump = NDVarArray(shape=1, dtype=object, buffer=np.array([assump]))
 
-    m = Model(hard+[assump.implies(candidates)]) # each assumption variable implies a candidate
-    s = SolverLookup.get(solver, m)
+    m = cp.Model(hard + [assump.implies(candidates)])  # each assumption variable implies a candidate
+    s = cp.SolverLookup.get(solver, m)
     assert not s.solve(assumptions=assump), "MUS: model must be UNSAT"
 
     mus = []
-    core = sorted(s.get_core()) # start from solver's UNSAT core
+    core = sorted(s.get_core())  # start from solver's UNSAT core
     for i in range(len(core)):
-        subassump = mus + core[i+1:]  # check if all but 'i' makes constraints SAT
-        
+        subassump = mus + core[i + 1:]  # check if all but 'i' makes constraints SAT
+
         if s.solve(assumptions=subassump):
             # removing it makes it SAT, must keep for UNSAT
             mus.append(core[i])
         # else: still UNSAT so don't need this candidate
-    
+
     # create dictionary from assump to candidate
-    dmap = dict(zip(assump,candidates))
+    dmap = dict(zip(assump, candidates))
     return [dmap[assump] for assump in mus]
-    
+
 
 def mus_naive(soft, hard=[], solver="ortools"):
     """
@@ -76,18 +77,117 @@ def mus_naive(soft, hard=[], solver="ortools"):
     # ensure toplevel list
     soft = toplevel_list(soft, merge_and=False)
 
-    m = Model(hard+soft)
+    m = cp.Model(hard + soft)
     assert not m.solve(solver=solver), "MUS: model must be UNSAT"
 
     mus = []
     # order so that constraints with many variables are tried and removed first
     core = sorted(soft, key=lambda c: -len(get_variables(c)))
     for i in range(len(core)):
-        subcore = mus + core[i+1:]  # check if all but 'i' makes core SAT
-        
-        if Model(hard+subcore).solve(solver=solver):
+        subcore = mus + core[i + 1:]  # check if all but 'i' makes core SAT
+
+        if cp.Model(hard + subcore).solve(solver=solver):
             # removing it makes it SAT, must keep for UNSAT
             mus.append(core[i])
         # else: still UNSAT so don't need this candidate
-    
+
     return mus
+
+
+def quickxplain(soft, hard=[], solver="ortools"):
+    """
+        Find a preferred minimal unsatisfiable subset of constraints, based on the ordering of constraints.
+
+        A total order is imposed on the constraints using the ordering of `soft`.
+        Constraints with lower index are preferred over ones with higher index
+
+        Assumption-based implementation for solvers that support s.solve(assumptions=...) and s.get_core()
+        More naive version available as `quickxplain_naive` to use with other solvers.
+
+        CPMpy implementation of the QuickXplain algorithm by Junker:
+            Junker, Ulrich. "Preferred explanations and relaxations for over-constrained problems." AAAI-2004. 2004.
+            https://cdn.aaai.org/AAAI/2004/AAAI04-027.pdf
+    """
+
+    soft = toplevel_list(soft, merge_and=False)
+
+    assump = cp.boolvar(shape=len(soft), name="assump")
+    if len(soft) == 1:
+        assump = NDVarArray(shape=1, dtype=object, buffer=np.array([assump]))
+
+    m = cp.SolverLookup.get(solver, cp.Model(hard))
+    m += assump.implies(soft)
+
+    assert m.solve(assumptions=assump) is False, "The model should be UNSAT!"
+    dmap = dict(zip(assump, soft))
+
+    # the recursive call
+    def do_recursion(soft, hard, delta):
+
+        if len(delta) != 0 and m.solve(assumptions=hard) is False:
+            # conflict is in hard constraints, no need to recurse
+            return []
+
+        if len(soft) == 1:
+            # conflict is not in hard constraints, but only 1 soft constraint
+            return list(soft)  # base case of recursion
+
+        split = len(soft) // 2  # determine split point
+        more_preferred, less_preferred = soft[:split], soft[split:]  # split constraints into two sets
+
+        # treat more preferred part as hard and find extra constants from less preferred
+        delta2 = do_recursion(less_preferred, hard + more_preferred, more_preferred)
+        # find which preferred constraints exactly
+        delta1 = do_recursion(more_preferred, hard + delta2, delta2)
+        return delta1 + delta2
+
+    # optimization: find max index of solver core
+    solver_core = set(m.get_core())
+    max_idx = max(i for i,a in enumerate(assump) if a in solver_core)
+
+    core = do_recursion(list(assump)[:max_idx+1], [], [])
+    return [dmap[a] for a in core]
+
+
+def quickxplain_naive(soft, hard=[], solver="ortools"):
+    """
+        Find a preferred minimal unsatisfiable subset of constraints, based on the ordering of constraints.
+
+        A total order is imposed on the constraints using the ordering of `soft`.
+        Constraints with lower index are preferred over ones with higher index
+
+        Naive implementation, re-solving the model from scratch.
+        Can be slower depending on the number of global constraints used and solver support for reified constraints.
+
+        CPMpy implementation of the QuickXplain algorithm by Junker:
+            Junker, Ulrich. "Preferred explanations and relaxations for over-constrained problems." AAAI-2004. 2004.
+            https://cdn.aaai.org/AAAI/2004/AAAI04-027.pdf
+    """
+
+
+    soft = toplevel_list(soft, merge_and=False)
+    assert cp.Model(hard + soft).solve(solver) is False, "The model should be UNSAT!"
+
+    # the recursive call
+    def do_recursion(soft, hard, delta):
+
+        m = cp.Model(hard)
+        if len(delta) != 0 and m.solve(solver) is False:
+            # conflict is in hard constraints, no need to recurse
+            return []
+
+        if len(soft) == 1:
+            # conflict is not in hard constraints, but only 1 soft constraint
+            return list(soft)  # base case of recursion
+
+        split = len(soft) // 2  # determine split point
+        more_preferred, less_preferred = soft[:split], soft[split:]  # split constraints into two sets
+
+        # treat more preferred part as hard and find extra constants from less preferred
+        delta2 = do_recursion(less_preferred, hard + more_preferred, more_preferred)
+        # find which preferred constraints exactly
+        delta1 = do_recursion(more_preferred, hard + delta2, delta2)
+        return delta1 + delta2
+
+    core = do_recursion(soft, hard, [])
+    return core

examples/advanced/quickxplain.py

@@ -0,0 +1,71 @@
+"""
+    CPMpy implementation of the QuickXplain algorithm by Junker.
+        Junker, Ulrich. "Preferred explanations and relaxations for over-constrained problems." AAAI-2004. 2004.
+        https://cdn.aaai.org/AAAI/2004/AAAI04-027.pdf
+"""
+
+import cpmpy as cp
+from cpmpy.transformations.get_variables import get_variables
+from cpmpy.transformations.normalize import toplevel_list
+
+
+def quickXplain(soft, hard=[], solver="ortools"):
+    """
+        Find a preferred minimal unsatisfiable subset of constraints
+        A partial order is imposed on the constraints using the order of `soft`.
+        Constraints with lower index are preferred over ones with higher index
+
+        :param: soft: soft constraints, list of expressions
+        :param: hard: hard constraints which cannot be relaxed, optional, list of expressions
+        :param: solver: name of a solver, see SolverLookup.solvernames()
+            "z3", "pysat" and "gurobi" are incremental, "ortools" restarts the solver
+    """
+
+    soft = toplevel_list(soft, merge_and=False)
+    assump = cp.boolvar(shape=len(soft))
+    solver = cp.SolverLookup.get(solver, cp.Model(hard))
+    solver += assump.implies(soft)
+
+    assert solver.solve(assumptions=assump) is False, "The model should be UNSAT!"
+
+    dmap = dict(zip(assump, soft))
+    core = recurse_explain(list(assump), [], [], solver=solver)
+    return [dmap[a] for a in core]
+
+
+def recurse_explain(soft, hard, delta, solver):
+    if len(delta) != 0 and solver.solve(assumptions=hard) is False:
+        # conflict is in hard constraints, no need to recurse
+        return []
+
+    if len(soft) == 1:
+        # conflict is not in hard constraints, but only 1 soft constraint
+        return list(soft)  # base case of recursion
+
+    split = len(soft) // 2  # determine split point
+    more_preferred, less_preferred = soft[:split], soft[split:]  # split constraints into two sets
+
+    # treat more preferred part as hard and find extra constants from less preferred
+    delta2 = recurse_explain(less_preferred, hard + more_preferred, more_preferred, solver=solver)
+    # find which preferred constraints exactly
+    delta1 = recurse_explain(more_preferred, hard + delta2, delta2, solver=solver)
+    return delta1 + delta2
+
+
+if __name__ == "__main__":
+    # example of quickXplain paper
+    options = {"roof racks": 500,
+               "CD-player": 500,
+               "one additional seat": 800,
+               "metal color": 500,
+               "special luxury version": 2600}
+
+    preferences = [cp.boolvar(name=name) for name in options.keys()]
+    p1, p2, p3, p4, p5 = preferences
+
+    hard = cp.sum(var * options[var.name] for var in preferences) <= 3000
+    core1 = quickXplain([p1, p2, p3, p4, p5], hard)
+    print("One cannot combine these options:", core1)
+
+    core2 = quickXplain([p3, p1, p2, p5, p4], hard)
+    print("One cannot combine these options:", core2)

tests/test_tools_mus.py

@@ -1,13 +1,18 @@
 import unittest
 from unittest import TestCase
 
-from cpmpy import *
-from cpmpy.tools.mus import mus, mus_naive
+import cpmpy as cp
+from cpmpy.tools.mus import mus, mus_naive, quickxplain, quickxplain_naive
 
 class MusTests(TestCase):
 
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.mus_func = mus
+        self.naive_func = mus_naive
+
     def test_circular(self):
-        x = intvar(0, 3, shape=4, name="x")
+        x = cp.intvar(0, 3, shape=4, name="x")
         # circular "bigger then", UNSAT
         cons = [
             x[0] > x[1], 
@@ -18,64 +23,94 @@ def test_circular(self):
             (x[3] > x[1]).implies((x[3] > x[2]) & ((x[3] == 3) | (x[1] == x[2])))
         ]
 
-        self.assertEqual(mus(cons), cons[:3])
-        self.assertEqual(mus_naive(cons), cons[:3])
+        self.assertEqual(self.mus_func(cons), cons[:3])
+        self.assertEqual(self.naive_func(cons), cons[:3])
 
     def test_bug_191(self):
         """
         Original Bug request: https://github.com/CPMpy/cpmpy/issues/191
         When assum is a single boolvar and candidates is a list (of length 1), it fails.
         """
-        bv = boolvar(name="x")
+        bv = cp.boolvar(name="x")
         hard = [~bv]
         soft = [bv]
 
-        mus_cons = mus(soft=soft, hard=hard) # crashes
+        mus_cons = self.mus_func(soft=soft, hard=hard) # crashes
         self.assertEqual(set(mus_cons), set(soft))
-        mus_naive_cons = mus_naive(soft=soft, hard=hard) # crashes
+        mus_naive_cons = self.naive_func(soft=soft, hard=hard) # crashes
         self.assertEqual(set(mus_naive_cons), set(soft))
 
     def test_bug_191_many_soft(self):
         """
         Checking whether bugfix 191  doesn't break anything in the MUS tool chain,
         when the number of soft constraints > 1.
         """
-        x = intvar(-9, 9, name="x")
-        y = intvar(-9, 9, name="y")
+        x = cp.intvar(-9, 9, name="x")
+        y = cp.intvar(-9, 9, name="y")
         hard = [x > 2]
         soft = [
             x + y < 6,
             y == 4
         ]
 
-        mus_cons = mus(soft=soft, hard=hard) # crashes
+        mus_cons = self.mus_func(soft=soft, hard=hard) # crashes
         self.assertEqual(set(mus_cons), set(soft))
-        mus_naive_cons = mus_naive(soft=soft, hard=hard) # crashes
+        mus_naive_cons = self.naive_func(soft=soft, hard=hard) # crashes
         self.assertEqual(set(mus_naive_cons), set(soft))
 
     def test_wglobal(self):
-        x = intvar(-9, 9, name="x")
-        y = intvar(-9, 9, name="y")
+        x = cp.intvar(-9, 9, name="x")
+        y = cp.intvar(-9, 9, name="y")
 
         cons = [
-            x < 0, 
-            x < 1,
+            x < 0,
             x > 2,
+            x < 1,
             y > 0,
             y == 4, 
             (x + y > 0) | (y < 0),
             (y >= 0) | (x >= 0),
             (y < 0) | (x < 0),
             (y > 0) | (x < 0),
-            AllDifferent(x,y) # invalid for musx_assum
+            cp.AllDifferent(x,y)
         ]
 
         # non-determinstic
         #self.assertEqual(set(mus(cons)), set(cons[1:3]))
-        ms = mus(cons)
+        ms = self.mus_func(cons)
+        self.assertLess(len(ms), len(cons))
+        self.assertFalse(cp.Model(ms).solve())
+        ms = self.naive_func(cons)
         self.assertLess(len(ms), len(cons))
-        self.assertFalse(Model(ms).solve())
-        self.assertEqual(set(mus_naive(cons)), set(cons[1:3]))
+        self.assertFalse(cp.Model(ms).solve())
+        # self.assertEqual(set(self.naive_func(cons)), set(cons[:2]))
+
+
+class QuickXplainTests(MusTests):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.mus_func = quickxplain
+        self.naive_func = quickxplain_naive
+
+    def test_prefered(self):
+
+        a,b,c,d = [cp.boolvar(name=n) for n in "abcd"]
+
+        mus1 = [b,d]
+        mus2 = [a,b,c]
+
+        hard = [~cp.all(mus1), ~cp.all(mus2)]
+        subset = self.mus_func([a,b,c,d],hard)
+        self.assertSetEqual(set(subset), {a,b,c})
+        subset2 = self.mus_func([d,c,b,a], hard)
+        self.assertSetEqual(set(subset2), {b,d})
+
+        subset = self.naive_func([a, b, c, d], hard)
+        self.assertSetEqual(set(subset), {a, b, c})
+        subset2 = self.naive_func([d, c, b, a], hard)
+        self.assertSetEqual(set(subset2), {b, d})
+
 
 if __name__ == '__main__':
     unittest.main()