Pysat wsum

cpmpy/solvers/pysat.py

@@ -33,7 +33,7 @@
 from ..expressions.core import Expression, Comparison, Operator
 from ..expressions.variables import _BoolVarImpl, NegBoolView, boolvar
 from ..expressions.utils import is_any_list, is_int
-from ..transformations.get_variables import get_variables
+from ..transformations.get_variables import get_variables, get_variables_model
 from ..transformations.to_cnf import to_cnf
 
 class CPM_pysat(SolverInterface):
@@ -64,6 +64,15 @@ def supported():
         except ImportError as e:
             return False
 
+    @staticmethod
+    def pb_supported():
+        try:
+            from pypblib import pblib
+            from pysat.pb import PBEnc
+            import pysat
+            return True
+        except ImportError as e:
+            return False
 
     @staticmethod
     def solvernames():
@@ -114,6 +123,7 @@ def __init__(self, cpm_model=None, subsolver=None):
         self.pysat_solver = Solver(use_timer=True, name=subsolver)
 
         # initialise everything else and post the constraints/objective
+        assert all(v.is_bool() for v in get_variables_model(cpm_model)), "Only support boolean variables in model"
         super().__init__(name="pysat:"+subsolver, cpm_model=cpm_model)
 
 
@@ -243,40 +253,41 @@ def _post_constraint(self, cpm_expr):
                 raise NotImplementedError(
                     f"Automatic conversion of Operator {cpm_expr} to CNF not yet supported, please report on github.")
         elif isinstance(cpm_expr, Comparison):
+            left = cpm_expr.args[0] # left-hand side, sum/wsum/mul
+            bound = cpm_expr.args[1] # right-hand side, constant
+
             # only handle cardinality encodings (for now)
-            if isinstance(cpm_expr.args[0], Operator) and cpm_expr.args[0].name == "sum" and all(
-                    isinstance(v, _BoolVarImpl) for v in cpm_expr.args[0].args):
-                lits = self.solver_vars(cpm_expr.args[0].args)
-                bound = cpm_expr.args[1]
+            if isinstance(left, Operator) and left.name == "sum" and all(isinstance(v, _BoolVarImpl) for v in left.args):
+                lits = self.solver_vars(left.args)
                 if not is_int(bound):
                     raise NotImplementedError(f"PySAT sum: rhs `{bound}` type {type(bound)} not supported")
 
                 clauses = []
                 if cpm_expr.name == "<":
-                    clauses += CardEnc.atmost(lits=lits, bound=bound - 1, vpool=self.pysat_vpool).clauses
+                    clauses += CardEnc.atmost(lits=lits, bound=bound-1, vpool=self.pysat_vpool).clauses
                 elif cpm_expr.name == "<=":
                     clauses += CardEnc.atmost(lits=lits, bound=bound, vpool=self.pysat_vpool).clauses
                 elif cpm_expr.name == ">=":
                     clauses += CardEnc.atleast(lits=lits, bound=bound, vpool=self.pysat_vpool).clauses
                 elif cpm_expr.name == ">":
-                    clauses += CardEnc.atleast(lits=lits, bound=bound + 1, vpool=self.pysat_vpool).clauses
+                    clauses += CardEnc.atleast(lits=lits, bound=bound+1, vpool=self.pysat_vpool).clauses
                 elif cpm_expr.name == "==":
                     clauses += CardEnc.equals(lits=lits, bound=bound, vpool=self.pysat_vpool).clauses
                 elif cpm_expr.name == "!=":
                     # special cases with bounding 'hardcoded' for clarity
                     if bound <= 0:
-                        clauses += CardEnc.atleast(lits=lits, bound=bound + 1, vpool=self.pysat_vpool).clauses
+                        clauses += CardEnc.atleast(lits=lits, bound=bound+1, vpool=self.pysat_vpool).clauses
                     elif bound >= len(lits):
-                        clauses += CardEnc.atmost(lits=lits, bound=bound - 1, vpool=self.pysat_vpool).clauses
+                        clauses += CardEnc.atmost(lits=lits, bound=bound-1, vpool=self.pysat_vpool).clauses
                     else:
                         ## add implication literal to atleast/atmost
                         is_atleast = self.solver_var(boolvar())
                         clauses += [atl + [-is_atleast] for atl in
-                                    CardEnc.atleast(lits=lits, bound=bound + 1, vpool=self.pysat_vpool).clauses]
+                                    CardEnc.atleast(lits=lits, bound=bound+1, vpool=self.pysat_vpool).clauses]
 
                         is_atmost = self.solver_var(boolvar())
                         clauses += [atm + [-is_atmost] for atm in
-                                    CardEnc.atmost(lits=lits, bound=bound - 1, vpool=self.pysat_vpool).clauses]
+                                    CardEnc.atmost(lits=lits, bound=bound-1, vpool=self.pysat_vpool).clauses]
 
                         ## add is_atleast or is_atmost
                         clauses.append([is_atleast, is_atmost])
@@ -285,9 +296,51 @@ def _post_constraint(self, cpm_expr):
 
                 # post the clauses
                 self.pysat_solver.append_formula(clauses)
-            else:
-                raise NotImplementedError(f"Non-operator constraint {cpm_expr} not supported by CPM_pysat")
 
+            # WEIGHTED !
+            elif isinstance(left, Operator) and (left.name in ["wsum", "mul"]) and is_int(bound):
+
+                if not CPM_pysat.pb_supported():
+                    raise ImportError("Please install PyPBLib: pip install pypblib")
+                from pysat.pb import PBEnc
+
+                if left.name == "mul":
+                    # single weight,value pair, in list
+                    weights = [left.args[0]]
+                    lits = [self.solver_var(left.args[1])]
+                else: # wsum
+                    weights = left.args[0]
+                    lits = self.solver_vars(left.args[1])
+
+                clauses = []
+                if cpm_expr.name == "<" or (cpm_expr.name == "!=" and bound >= len(lits)):
+                    clauses += PBEnc.leq(lits=lits, weights=weights, bound=bound-1, vpool=self.pysat_vpool).clauses
+                elif cpm_expr.name == "<=":
+                    clauses += PBEnc.leq(lits=lits, weights=weights, bound=bound,vpool=self.pysat_vpool).clauses
+                elif cpm_expr.name == ">" or (cpm_expr.name == "!=" and bound <= 0):
+                    clauses += PBEnc.geq(lits=lits, weights=weights, bound=bound+1, vpool=self.pysat_vpool).clauses
+                elif cpm_expr.name == ">=":
+                    clauses += PBEnc.geq(lits=lits, weights=weights, bound=bound ,vpool=self.pysat_vpool).clauses
+                elif cpm_expr.name == "==":
+                    clauses +=PBEnc.equals(lits=lits, weights=weights, bound=bound, vpool=self.pysat_vpool)
+
+                elif cpm_expr.name == "!=":
+                    # BUG with pblib solved in Pysat dev 0.1.7.dev12
+                    is_atleast = self.solver_var(boolvar())
+                    clauses += [atl + [-is_atleast] for atl in
+                                PBEnc.geq(lits=lits, weights=weights, bound=bound+1, vpool=self.pysat_vpool).clauses]
+                    is_atmost = self.solver_var(boolvar())
+                    clauses += [atm + [-is_atmost] for atm in
+                                PBEnc.leq(lits=lits, weights=weights, bound=bound-1, vpool=self.pysat_vpool).clauses]
+
+                    ## add is_atleast or is_atmost
+                    clauses.append([is_atleast, is_atmost])
+                else:
+                    raise NotImplementedError(f"Comparison: {cpm_expr} not supported by CPM_pysat")
+
+                self.pysat_solver.append_formula(clauses)
+            else:
+                raise NotImplementedError(f"Comparison: {cpm_expr} not supported by CPM_pysat")
         else:
             raise NotImplementedError(f"Non-operator constraint {cpm_expr} not supported by CPM_pysat")
 

cpmpy/transformations/flatten_model.py

@@ -81,6 +81,8 @@
 import math
 import numpy as np
 from ..expressions.core import *
+from ..expressions.core import _wsum_should
+from ..expressions.core import _wsum_make
 from ..expressions.variables import _NumVarImpl, _IntVarImpl, _BoolVarImpl, NegBoolView
 from ..expressions.utils import is_num, is_any_list
 
@@ -185,6 +187,17 @@ def flatten_constraint(expr):
     - Numeric inequality (>=,>,<,<=,): Numexpr >=< Var     (CPMpy class 'Comparison')
     - Reification (double implication): Boolexpr == Var    (CPMpy class 'Comparison')
         """
+
+        left, right = expr.args[0], expr.args[1]
+
+        # Flatten a complex weighted sum where any of the sub expressions is a sum, mul, neg, ...
+        # e.g. bv0 - 3 * (bv2 + 2 * bv1)
+        if isinstance(left, Operator) and \
+            (left.name == "wsum" or
+            any(_wsum_should_flatten(sub_expr) for sub_expr in left.args) ):
+            w_new, x_new = _wsum_make_flatten(left)
+            return [Comparison(expr.name, Operator("wsum",[w_new, x_new] ),right)]
+
         if all(__is_flat_var(arg) for arg in expr.args):
             return [expr]
 
@@ -393,7 +406,6 @@ def normalized_boolexpr(expr):
     """
     assert(not __is_flat_var(expr))
     assert(expr.is_bool()) 
-
     if isinstance(expr, Operator):
         # and, or, ->
 
@@ -621,3 +633,46 @@ def negated_normal(expr):
             return negated_normal(Operator('and', expr.decompose()))
         else:
             raise NotImplementedError("negate_normal {}".format(expr))
+
+
+def _wsum_should_flatten(arg):
+    # in flatten we also turn negations (in sums) into weighted sums,
+    # for stronger expression simplification
+    if isinstance(arg, Operator) and arg.name == "-":
+        # TODO: HANDLE the ABS operator 
+        if any(isinstance(sub_arg, Operator) and sub_arg.name == "abs" for sub_arg in arg.args):
+            return False
+        return True
+    return _wsum_should(arg)
+
+def _wsum_make_flatten(sub_expr):
+    # Mixed operator weighted sums that can be flattened into one,
+    # where any of the sub expressions is a sum, mul, neg, ...
+    # e.g. bv0 - 3 * (bv2 + 2 * bv1)
+    if sub_expr.name == 'wsum':
+        w_new, x_new = [], []
+        for wi, xi in zip(sub_expr.args[0], sub_expr.args[1]):
+            wni, xni = _wsum_make_flatten(xi)
+            wni = [wnij * wi for wnij in wni]
+            w_new += wni
+            x_new += xni
+        return w_new, x_new
+    elif sub_expr.name == 'mul':
+        w = sub_expr.args[0]
+        wi, x = _wsum_make_flatten(sub_expr.args[1])
+        wi_new = [wij*w for wij in wi]
+        return wi_new, x
+    elif sub_expr.name == "sum":
+        w_new, x_new = [], []
+        for xi in sub_expr.args:
+            wni, xni = _wsum_make_flatten(xi)
+            w_new += wni
+            x_new += xni
+        return w_new, x_new
+    elif sub_expr.name == "-" and isinstance(sub_expr.args[0], Operator):
+        # - (3 * y)
+        w, x = _wsum_make_flatten(sub_expr.args[0])
+        return [-i for i in w], x
+
+    ## base case
+    return _wsum_make(sub_expr)

tests/test_constraints.py

@@ -15,7 +15,7 @@
 # Exclude certain operators for solvers.
 # Not all solvers support all operators in CPMpy
 EXCLUDE_OPERATORS = {"gurobi": {"mod"},
-                     "pysat": {"sum", "wsum", "sub", "mod", "div", "pow", "abs", "mul","-"}}
+                     "pysat": {"sub", "mod", "div", "pow", "abs", "mul","-"}}
 
 # Some solvers only support a subset of operators in imply-constraints
 # This subset can differ between left and right hand side of the implication

tests/test_pysat_wsum.py

@@ -0,0 +1,78 @@
+import unittest
+import cpmpy as cp 
+from cpmpy import *
+from cpmpy.solvers.pysat import CPM_pysat
+from cpmpy.transformations.to_cnf import to_cnf
+
+class TestEncodeLinearConstraint(unittest.TestCase):
+    def setUp(self):
+        self.bv = boolvar(shape=3)
+
+    def test_pysat_simple_atmost(self):
+
+        atmost = cp.Model(
+            ## < This does not work
+            - 2 * self.bv[0] < 3,
+            ## <=
+            - 3 * self.bv[1] <= 3,
+            ## >
+            2 * self.bv[2] > 1,
+            ## >=
+            4 * self.bv[2] >= 3,
+        )
+        ps = CPM_pysat(atmost)
+        ps.solve()
+
+    def test_pysat_boolean_linear_sum(self):
+        ls = cp.Model(
+            2 * self.bv[0] + 3 * self.bv[1] <= 3,
+            self.bv[1] == 1,
+        )
+
+        ps = CPM_pysat(ls)
+        ps.solve()
+
+    def test_pysat_complex_expressions(self):
+        ls = cp.Model(
+            2 * self.bv[0] + 3 * self.bv[1] <= 3,
+            self.bv[1] == 1,
+            self.bv[2] == 0
+        )
+
+        ps = CPM_pysat(ls)
+        ps.solve()
+
+
+    def test_pysat_unsat(self):
+        ls = cp.Model(
+            2 * self.bv[0] + 3 * self.bv[1] <= 3,
+            self.bv[0] == 1,
+            self.bv[1] == 1
+        )
+
+        ps = CPM_pysat(ls)
+        solved = ps.solve()
+        self.assertFalse(solved)
+
+    def test_encode_linear_expressions(self):
+        expressions = [
+            - self.bv[2] == -1,
+            - 2 * self.bv[2] == -2,
+            self.bv[0] - self.bv[2] > 0,
+            -self.bv[0] + self.bv[2] > 0,
+            2 * self.bv[0] + 3 * self.bv[2] > 0,
+            2 * self.bv[0] - 3 * self.bv[2] + 2 * self.bv[1] > 0,
+            self.bv[0] - 3 * self.bv[2] > 0,
+            self.bv[0] - 3 * (self.bv[2] + 2 * self.bv[1])> 0,
+        ]
+
+        ## check all types of linear constraints are handled
+        for expression in expressions:
+            ps = CPM_pysat(Model(
+                expression
+            ))
+            ps.solve()
+
+if __name__ == '__main__':
+    unittest.main()
+