From bb5129ac698ed3601f0059b8720d946314fdb2ab Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Fri, 27 Mar 2026 20:25:19 +0000 Subject: [PATCH 01/23] Add initial tree w/ mostly reification --- cpmpy/solvers/gurobi.py | 249 ++++++++++++++++++++++++++++++++++++-- tests/test_constraints.py | 82 +++++++------ 2 files changed, 285 insertions(+), 46 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 26fc3bda4..a107b63a2 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -42,6 +42,7 @@ ============== """ +import cpmpy as cp from typing import Optional, List from .solver_interface import SolverInterface, SolverStatus, ExitStatus, Callback @@ -191,6 +192,9 @@ def solve(self, time_limit:Optional[float]=None, solution_callback=None, **kwarg for param, val in kwargs.items(): self.grb_model.setParam(param, val) + # write LP file for debugging + self.grb_model.write("/tmp/model.lp") + _ = self.grb_model.optimize(callback=solution_callback) grb_objective = self.grb_model.getObjective() @@ -258,7 +262,7 @@ def solver_var(self, cpm_var): # special case, negative-bool-view. Should be eliminated in linearize if isinstance(cpm_var, NegBoolView): - raise NotSupportedError("Negative literals should not be left as part of any equation. Please report.") + return 1 - self.solver_var(cpm_var._bv) # create if it does not exit if cpm_var not in self._varmap: @@ -362,15 +366,17 @@ def transform(self, cpm_expr): supported=self.supported_global_constraints, supported_reified=self.supported_reified_global_constraints, csemap=self._csemap) - cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap) # flat normal form - cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification - cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != - cpm_cons = linearize_reified_variables(cpm_cons, min_values=2, csemap=self._csemap) - cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) - cpm_cons = only_implies(cpm_cons, csemap=self._csemap) # anything that can create full reif should go above... + # cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap) # flat normal form + # cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification + # cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != + + # cpm_cons = linearize_reified_variables(cpm_cons, min_values=2, csemap=self._csemap) + # cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) + # cpm_cons = only_implies(cpm_cons, csemap=self._csemap) # anything that can create full reif should go above... + # gurobi does not round towards zero, so no 'div' in supported set: https://github.com/CPMpy/cpmpy/pull/593#issuecomment-2786707188 - cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization - cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv + # cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization + # cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv return cpm_cons def add(self, cpm_expr_orig): @@ -391,13 +397,231 @@ def add(self, cpm_expr_orig): :return: self """ - from gurobipy import GRB + from gurobipy import GRB, nlfunc + import gurobipy as gp + + def add(cpm_expr, depth=0): + """Recursively create a Gurobi expression tree from a CPMpy expression.""" + indent = " " * depth + depth+=1 + self.grb_model.update() + print(f"{indent}Con:", cpm_expr, type(cpm_expr)) + # Base case: numbers + if is_num(cpm_expr): + return int(cpm_expr) + + # # Base case: variables + # if isinstance(cpm_expr, NegBoolView): + # assert False + # # return 1 - self.solver_var(cpm_expr) + + def reify(grb_expr): + """Create a binary variable representing the truth value of cpm_expr.""" + + self.grb_model.update() + print(f"{indent}Reif", grb_expr) + + + # # Base case: already a variable + # if isinstance(cpm_expr, _IntVarImpl) and not isinstance(cpm_expr, NegBoolView): + # return self.solver_var(cpm_expr) + + # Base case: already a variable + if isinstance(grb_expr, gp.Var): + return grb_expr + + + r = self.solver_var(cp.boolvar()) + + self.grb_model.update() + print(indent, "R", r, grb_expr) + + # Double-reification: r == TempConstr works for linear constraints + self.grb_model.addConstr(r == grb_expr) + + + # print("R", r, grb_expr) + # if isinstance(grb_expr, gp.TempConstr): + # # self.grb_model.addConstr(r == (grb_expr)) + # self.grb_model.addConstr((r==1) >> (grb_expr)) + # else: + # # For other expressions, add() returns something we can equate + # # grb_expr = add(cpm_expr) + # self.grb_model.addConstr(r == grb_expr) + # return r + + + + return r + + def reify_cmp(sense, lhs_expr, rhs_const): + """Reify a comparison using indicator constraints.""" + r = self.solver_var(cp.boolvar()) + self.grb_model.update() + self.grb_model.addGenConstrIndicator(r, True, lhs_expr, sense, rhs_const) + # Add negation for double-reification + if sense == GRB.LESS_EQUAL: + self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.GREATER_EQUAL, rhs_const + 1) + elif sense == GRB.GREATER_EQUAL: + self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.LESS_EQUAL, rhs_const - 1) + elif sense == GRB.EQUAL: + self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.LESS_EQUAL, rhs_const - 1) + return r + + + # Base case: variables + if isinstance(cpm_expr, NegBoolView): + return reify(1 - self.solver_var(cpm_expr._bv)) + if isinstance(cpm_expr, _NumVarImpl): + return self.solver_var(cpm_expr) + + # def filter_(args): + # return gp.and_(add(arg) for arg in cpm_expr.args) + + if isinstance(cpm_expr, Operator): + match cpm_expr.name: + case "-": + return -add(cpm_expr.args[0]) + case "or": + return reify(gp.or_(add(arg) for arg in cpm_expr.args)) + case "and": + # return reify(gp.and_([0,1])) + args = [add(arg) for arg in cpm_expr.args if not is_num(arg) or arg != 1] + assert len(args) + if not args: + return 1 + elif any(a == 0 for a in args if is_num(a)): + return 0 + else: + return reify(gp.and_(args)) + case "->": + # TODO has to be wrong ; reify(...) + a, b = cpm_expr.args + lhs, rhs = add(a, depth), add(b, depth) + # return add(a, depth + 1) <= add(b, depth + 1) + self.grb_model.update() + # self.grb_model.addConstr((lhs==1) >> (rhs>=1)) + # return lhs + # return reify((lhs==1) >> (rhs>=1)) + return reify_cmp(GRB.LESS_EQUAL, lhs - rhs, 0) + case "not": + # Negation: not(a) is 1 - a + return 1 - add(cpm_expr.args[0], depth) + case "sum": + return gp.quicksum(add(arg, depth) for arg in cpm_expr.args) + case "wsum": + weights, args = cpm_expr.args + return gp.quicksum(w * add(arg, depth) for w, arg in zip(weights, args)) + case "sub": + a, b = cpm_expr.args + return add(a, depth) - add(b, depth) + elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): + match cpm_expr.name: + case "pow": + # y = x^z + x, z = add(cpm_expr.args[0]), add(cpm_expr.args[1]) + # y = self.solver_var(cp.boolvar()) + # self.grb_model.addGenConstrPow(x, y, z) + # self.grb_model.addConstr(y == (x**z)) + # self.grb_model.addGenConstrNL(x, y, z) + # return reify(x**z) + return x**z + case "abs": + return reify(gp.abs_(add(cpm_expr.args[0], depth))) + case "min": + return reify(gp.min_(add(arg, depth) for arg in cpm_expr.args)) + case "max": + return reify(gp.max_(add(arg, depth) for arg in cpm_expr.args)) + case "mul": + a, b = cpm_expr.args + # return reify(add(a)*add(b)) + return add(a) * add(b) + elif isinstance(cpm_expr, cp.expressions.globalconstraints.GlobalConstraint): + match cpm_expr.name: + case "ite": + i, t, e = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth), add(cpm_expr.args[2], depth) + print('ite', i,t,e) + self.grb_model.addConstr((i==1) >> (t>=1)) + self.grb_model.addConstr((i==0) >> (e>=1)) + return i + + if isinstance(cpm_expr, Comparison): + lhs, rhs = cpm_expr.args + grb_lhs, grb_rhs = add(lhs, depth), add(rhs, depth) + self.grb_model.update() # Ensure vars are in model before creating TempConstr + + # Rewrite to form: (grb_lhs - grb_rhs) sense constant + # addGenConstrIndicator requires rhs to be a constant + lhs_expr = grb_lhs - grb_rhs # LinExpr + + match cpm_expr.name: + case "==": + # lhs == rhs: r=1 -> lhs==rhs, r=0 -> lhs!=rhs + r = self.solver_var(cp.boolvar()) + self.grb_model.update() + # r=1 -> lhs == rhs + # self.grb_model.addGenConstrIndicator(r, True, lhs_expr, GRB.EQUAL, 0) + self.grb_model.addConstr((r==1) >> (lhs_expr == 0)) + # r=0 -> lhs != rhs (reuse the != case) + r_neq = add(lhs != rhs) + self.grb_model.addConstr((r==0) >> (r_neq >= 1)) + # self.grb_model.addConstr(r + r_neq >= 1) + return r + case "!=": + # lhs != rhs => (lhs < rhs) | (lhs > rhs) + r_lt = reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) # lhs - rhs <= -1 + r_gt = reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) # lhs - rhs >= 1 + r = self.solver_var(cp.boolvar()) + self.grb_model.update() + self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) + return r + case "<=": + # lhs <= rhs => (lhs - rhs) <= 0 + return lhs_expr <= 0 + + # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) + case ">=": + # lhs >= rhs => (lhs - rhs) >= 0 + return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 0) + case "<": + # lhs < rhs => (lhs - rhs) <= -1 + return reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) + case ">": + # lhs > rhs => (lhs - rhs) >= 1 + return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) + + raise NotImplementedError(f"add() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") + # add new user vars to the set get_variables(cpm_expr_orig, collect=self.user_vars) # transform and post the constraints for cpm_expr in self.transform(cpm_expr_orig): + print("E", cpm_expr) + grb_expr = add(cpm_expr) + self.grb_model.update() + print("out", grb_expr, type(grb_expr)) + # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 + if isinstance(grb_expr, (gp.Var, gp.LinExpr)): + self.grb_model.addConstr(grb_expr >= 1) + elif isinstance(grb_expr, gp.TempConstr): + self.grb_model.addConstr(grb_expr) + elif type(grb_expr).__name__ == 'GenExprAnd': + # GenExprAnd needs addGenConstrAnd: r = AND(args) then r >= 1 + r = self.solver_var(cp.boolvar()) + self.grb_model.addGenConstrAnd(r, list(grb_expr)) + self.grb_model.addConstr(r >= 1) + elif type(grb_expr).__name__ == 'GenExprOr': + # GenExprOr needs addGenConstrOr: r = OR(args) then r >= 1 + r = self.solver_var(cp.boolvar()) + self.grb_model.addGenConstrOr(r, list(grb_expr)) + self.grb_model.addConstr(r >= 1) + else: + r = self.solver_var(cp.boolvar()) + self.grb_model.addConstr(r == grb_expr) + self.grb_model.addConstr(r >= 1) + continue # Comparisons: only numeric ones as 'only_implies()' has removed the '==' reification for Boolean expressions # numexpr `comp` bvar|const @@ -486,8 +710,11 @@ def add(self, cpm_expr_orig): elif isinstance(cpm_expr, DirectConstraint): cpm_expr.callSolver(self, self.grb_model) + elif isinstance(cpm_expr, _BoolVarImpl): + self.grb_model.addConstr(self.solver_var(cpm_expr) >= 1) + else: - raise NotImplementedError(cpm_expr) # if you reach this... please report on github + raise NotImplementedError(f"Please report unsupported constraint in Gurobi interface: {cpm_expr} of type {type(cpm_expr)}") # if you reach this... please report on github return self __add__ = add # avoid redirect in superclass diff --git a/tests/test_constraints.py b/tests/test_constraints.py index c10736014..c922e50a2 100644 --- a/tests/test_constraints.py +++ b/tests/test_constraints.py @@ -17,12 +17,14 @@ # also add exclusions to the 3 EXCLUDE_* below as needed SOLVERNAMES = [name for name, solver in SolverLookup.base_solvers() if solver.supported()] ALL_SOLS = False # test whether all solutions returned by the solver satisfy the constraint +# ALL_SOLS = True # test whether all solutions returned by the solver satisfy the constraint # Exclude some global constraints for solvers NUM_GLOBAL = { "AllEqual", "AllDifferent", "AllDifferentExcept0", "AllDifferentExceptN", "AllEqualExceptN", - "GlobalCardinalityCount", "InDomain", "Inverse","Circuit", + "GlobalCardinalityCount", "InDomain", "Inverse", + # "Circuit", "Table", 'NegativeTable', "ShortTable", "Regular", "Increasing", "IncreasingStrict", "Decreasing", "DecreasingStrict", "Precedence", "Cumulative", "NoOverlap", "CumulativeOptional", "NoOverlapOptional", @@ -37,6 +39,7 @@ EXCLUDE_GLOBAL = { "pysdd": NUM_GLOBAL | {"Xor"}, "minizinc": {"IncreasingStrict"}, # bug #813 reported on libminizinc + "gurobi": {"Circuit", "CumulativeOptional"} } @@ -284,6 +287,14 @@ def global_functions(solver): else: yield cls(NUM_ARGS) +def generate_cases(solver): + yield cp.boolvar(name="x") >= 0 # issue #736 + x, y = cp.intvar(1, 3,shape=2, name=["x", "y"]) + yield x ** 2 - 2*x*y + y**2 <= 3 + + # p, q = cp.intvar(shape=2, name=["p", "q"]) + # yield p + # yield ((cp.boolvar(name="x") >= 0) | (cp.boolvar(name="y") >= 0)) # issue #736 def reify_imply_exprs(solver): """ @@ -304,48 +315,49 @@ def reify_imply_exprs(solver): def verify(cons): - assert argval(cons) - assert cons.value() - -@pytest.mark.generate_constraints.with_args(bool_exprs) + from cpmpy.transformations.get_variables import get_variables + vars_ = get_variables(cons) + assignment = {v.name: v.value() for v in sorted(vars_, key=lambda v: v.name)} + assert argval(cons), f"argval failed for {cons} with assignment {assignment}" + assert cons.value(), f"value() failed for {cons} with assignment {assignment}" + +def all_constraints(solver): + """Combined generator for all constraint types.""" + # yield from bool_exprs(solver) + # yield from comp_constraints(solver) + # yield from reify_imply_exprs(solver) + yield from generate_cases(solver) + +@pytest.mark.generate_constraints.with_args(all_constraints) @skip_on_missing_pblib(skip_on_exception_only=True) -def test_bool_constraints(solver, constraint): +def test_constraints(solver, constraint): """ - Tests boolean constraint by posting it to the solver and checking the value after solve. + Tests constraint by posting it to the solver and checking the value after solve. """ if ALL_SOLS: - n_sols = SolverLookup.get(solver, Model(constraint)).solveAll(display=lambda: verify(constraint)) + n_sols = SolverLookup.get(solver, Model(constraint)).solveAll(display=lambda: verify(constraint), solution_limit=100) assert n_sols >= 1 else: assert SolverLookup.get(solver, Model(constraint)).solve() assert argval(constraint) assert constraint.value() -@pytest.mark.generate_constraints.with_args(comp_constraints) -@skip_on_missing_pblib(skip_on_exception_only=True) -def test_comparison_constraints(solver, constraint): - """ - Tests comparison constraint by posting it to the solver and checking the value after solve. - """ - if ALL_SOLS: - n_sols = SolverLookup.get(solver, Model(constraint)).solveAll(display= lambda: verify(constraint)) - assert n_sols >= 1 - else: - assert SolverLookup.get(solver,Model(constraint)).solve() - assert argval(constraint) - assert constraint.value() - +if __name__ == "__main__": + solver = None # Use None for no solver-specific exclusions + + generators = [ + ("Boolean expressions", bool_exprs), + ("Comparison constraints", comp_constraints), + ("Global constraints", global_constraints), + ("Global functions", global_functions), + ("Reify/imply expressions", reify_imply_exprs), + ] + + for name, gen in generators: + print(f"\n{'='*60}") + print(f"{name}") + print('='*60) + for i, expr in enumerate(gen(solver)): + model = Model(expr) + print(f"{i+1}. {model}") -@pytest.mark.generate_constraints.with_args(reify_imply_exprs) -@skip_on_missing_pblib(skip_on_exception_only=True) -def test_reify_imply_constraints(solver, constraint): - """ - Tests boolean expression by posting it to solver and checking the value after solve. - """ - if ALL_SOLS: - n_sols = SolverLookup.get(solver, Model(constraint)).solveAll(display=lambda: verify(constraint)) - assert n_sols >= 1 - else: - assert SolverLookup.get(solver, Model(constraint)).solve() - assert argval(constraint) - assert constraint.value() From 8e959e81f82e396012fd04cd04e259127a03b7c9 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 31 Mar 2026 08:35:21 +0000 Subject: [PATCH 02/23] Progress --- cpmpy/solvers/gurobi.py | 174 +++++++++++++++++-------- cpmpy/transformations/flatten_model.py | 32 +++-- 2 files changed, 140 insertions(+), 66 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index a107b63a2..651dbc9c8 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -45,6 +45,8 @@ import cpmpy as cp from typing import Optional, List +import math + from .solver_interface import SolverInterface, SolverStatus, ExitStatus, Callback from ..exceptions import NotSupportedError from ..expressions.core import * @@ -52,9 +54,10 @@ from ..expressions.variables import _BoolVarImpl, NegBoolView, _IntVarImpl, _NumVarImpl, intvar from ..expressions.globalconstraints import DirectConstraint from ..transformations.comparison import only_numexpr_equality -from ..transformations.flatten_model import flatten_constraint, flatten_objective +from ..transformations.flatten_model import flatten_constraint, flatten_objective, get_or_make_var_or_list from ..transformations.get_variables import get_variables from ..transformations.linearize import linearize_constraint, linearize_reified_variables, only_positive_bv, only_positive_bv_wsum, decompose_linear, decompose_linear_objective +from ..transformations.decompose_global import decompose_in_tree from ..transformations.normalize import toplevel_list from ..transformations.reification import only_implies, reify_rewrite, only_bv_reifies from ..transformations.safening import no_partial_functions, safen_objective @@ -300,7 +303,7 @@ def objective(self, expr, minimize=True): supported=self.supported_global_constraints, supported_reified=self.supported_reified_global_constraints, csemap=self._csemap) - obj, flat_cons = flatten_objective(obj, csemap=self._csemap) + obj, flat_cons = flatten_objective(obj, csemap=self._csemap, supported={"pow", "mul"}) obj = only_positive_bv_wsum(obj) # remove negboolviews self.add(safe_cons + decomp_cons + flat_cons) @@ -362,11 +365,15 @@ def transform(self, cpm_expr): # expressions have to be linearized to fit in MIP model. See /transformations/linearize cpm_cons = toplevel_list(cpm_expr) cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"}) # linearize and decompose expect safe exprs + cpm_cons = decompose_in_tree(cpm_cons, + supported=self.supported_global_constraints, + supported_reified=self.supported_reified_global_constraints, + csemap=self._csemap) cpm_cons = decompose_linear(cpm_cons, supported=self.supported_global_constraints, supported_reified=self.supported_reified_global_constraints, csemap=self._csemap) - # cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap) # flat normal form + cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap, supported={"mul", "pow"}) # flat normal form # cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification # cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != @@ -429,12 +436,10 @@ def reify(grb_expr): # Base case: already a variable if isinstance(grb_expr, gp.Var): return grb_expr - - r = self.solver_var(cp.boolvar()) + self.grb_model.update() - print(indent, "R", r, grb_expr) # Double-reification: r == TempConstr works for linear constraints self.grb_model.addConstr(r == grb_expr) @@ -456,16 +461,44 @@ def reify(grb_expr): def reify_cmp(sense, lhs_expr, rhs_const): """Reify a comparison using indicator constraints.""" + + # (2x _3y <= 5) + # (L == 2x + 3y, r == (L <= 5)) + + r_ = self.solver_var(cp.boolvar()) + + # r <-> (x = 3) + # L = x, r -> L = 3, ~r -> (L >= 4 \/ L <= 2) + # + + # TODO if not already LinExpr + L = self.grb_model.addVar(vtype=GRB.INTEGER, lb=-GRB.INFINITY,ub=GRB.INFINITY) + self.grb_model.addConstr(L == lhs_expr) r = self.solver_var(cp.boolvar()) - self.grb_model.update() - self.grb_model.addGenConstrIndicator(r, True, lhs_expr, sense, rhs_const) + + self.grb_model.addGenConstrIndicator(r, True, L, sense, rhs_const) # Add negation for double-reification - if sense == GRB.LESS_EQUAL: - self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.GREATER_EQUAL, rhs_const + 1) - elif sense == GRB.GREATER_EQUAL: - self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.LESS_EQUAL, rhs_const - 1) - elif sense == GRB.EQUAL: - self.grb_model.addGenConstrIndicator(r, False, lhs_expr, GRB.LESS_EQUAL, rhs_const - 1) + match sense: + case GRB.LESS_EQUAL: + self.grb_model.addGenConstrIndicator(r, False, L, GRB.GREATER_EQUAL, rhs_const + 1) + case GRB.GREATER_EQUAL: + self.grb_model.addGenConstrIndicator(r, False, L, GRB.LESS_EQUAL, rhs_const - 1) + case GRB.EQUAL: + # self.grb_model.addGenConstrIndicator(r, False, L, GRB.LESS_EQUAL, rhs_const - 1) + # D = self.grb_model.addVar(vtype=GRB.INTEGER, lb=-GRB.INFINITY,ub=GRB.INFINITY) + # self.grb_model.addConstr(D == gp.abs_(L - lhs_expr)) + # self.grb_model.addConstr((r==0) >> (D >= 1)) + + a = reify_cmp(GRB.GREATER_EQUAL, L, rhs_const + 1) + b = reify_cmp(GRB.LESS_EQUAL, L, rhs_const - 1) + print('neq', a,b) + self.grb_model.addConstr((r==0) >> (a + b >= 1)) + + # # self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) + # r_ = self.solver_var(cp.boolvar()) + # self.grb_model.addConstr((r==0) >> r_) + # self.grb_model.addConstr(r_ == ) + return r @@ -480,20 +513,32 @@ def reify_cmp(sense, lhs_expr, rhs_const): if isinstance(cpm_expr, Operator): match cpm_expr.name: - case "-": - return -add(cpm_expr.args[0]) case "or": + # return reify(gp.or_(add(arg) for arg in cpm_expr.args)) + self.grb_model.update() + + # return sum(add(arg, depth=depth) for arg in cpm_expr.args) + # def normalize(arg): + # if isinstance(arg, gp.QuadExpr): + # self.native_model + + # args = [normalize(arg) for arg in args] return reify(gp.or_(add(arg) for arg in cpm_expr.args)) case "and": + self.grb_model.update() + return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) + # return reify(gp.and_([0,1])) - args = [add(arg) for arg in cpm_expr.args if not is_num(arg) or arg != 1] + args = [add(arg, depth=depth) for arg in cpm_expr.args if not is_num(arg) or arg != 1] assert len(args) if not args: return 1 elif any(a == 0 for a in args if is_num(a)): return 0 else: - return reify(gp.and_(args)) + # return gp.and_(args) + # return sum(args) >= len(args) + return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) case "->": # TODO has to be wrong ; reify(...) a, b = cpm_expr.args @@ -507,11 +552,15 @@ def reify_cmp(sense, lhs_expr, rhs_const): case "not": # Negation: not(a) is 1 - a return 1 - add(cpm_expr.args[0], depth) + case "-": + return -add(cpm_expr.args[0], depth=depth) case "sum": + return sum(add(arg, depth) for arg in cpm_expr.args) return gp.quicksum(add(arg, depth) for arg in cpm_expr.args) case "wsum": weights, args = cpm_expr.args - return gp.quicksum(w * add(arg, depth) for w, arg in zip(weights, args)) + return sum(w * add(arg, depth) for w, arg in zip(weights, args)) + # return gp.quicksum(w * add(arg, depth) for w, arg in zip(weights, args)) case "sub": a, b = cpm_expr.args return add(a, depth) - add(b, depth) @@ -526,24 +575,40 @@ def reify_cmp(sense, lhs_expr, rhs_const): # self.grb_model.addGenConstrNL(x, y, z) # return reify(x**z) return x**z - case "abs": - return reify(gp.abs_(add(cpm_expr.args[0], depth))) - case "min": - return reify(gp.min_(add(arg, depth) for arg in cpm_expr.args)) - case "max": - return reify(gp.max_(add(arg, depth) for arg in cpm_expr.args)) case "mul": a, b = cpm_expr.args # return reify(add(a)*add(b)) return add(a) * add(b) - elif isinstance(cpm_expr, cp.expressions.globalconstraints.GlobalConstraint): - match cpm_expr.name: - case "ite": - i, t, e = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth), add(cpm_expr.args[2], depth) - print('ite', i,t,e) - self.grb_model.addConstr((i==1) >> (t>=1)) - self.grb_model.addConstr((i==0) >> (e>=1)) - return i + case _: + # flatargs, flatcons = zip(*[get_or_make_var_or_list(arg, csemap=self._csemap) for arg in cpm_expr.args]) + flatargs = [add(arg, depth) for arg in cpm_expr.args] + self.grb_model.update() + # self.add(flatcons) + match cpm_expr.name: + case "abs": # y = abs(x) + x, = flatargs + x_, = cpm_expr.args + # y = add(cp.intvar(lb=max(0,cpm_expr.args[0].lb), ub=cpm_expr.args[0].ub), depth) + y = add(cp.intvar(lb=x_.lb, ub=x_.ub), depth) + self.grb_model.update() + # self.native_model.addConstr(self.solver_var(y) == gp.abs_(add(x, depth))) + self.native_model.addConstr(y == gp.abs_(x)) + return y + + # return gp.abs_(add(cpm_expr.args[0], depth)) + # return reify(abs(add(cpm_expr.args[0], depth))) + case "min" | "max": + y = add(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) + self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) + return y + # elif isinstance(cpm_expr, cp.expressions.globalconstraints.GlobalConstraint): + # match cpm_expr.name: + # case "ite": + # i, t, e = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth), add(cpm_expr.args[2], depth) + # print('ite', i,t,e) + # self.grb_model.addConstr((i==1) >> (t>=1)) + # self.grb_model.addConstr((i==0) >> (e>=1)) + # return i if isinstance(cpm_expr, Comparison): lhs, rhs = cpm_expr.args @@ -556,16 +621,14 @@ def reify_cmp(sense, lhs_expr, rhs_const): match cpm_expr.name: case "==": - # lhs == rhs: r=1 -> lhs==rhs, r=0 -> lhs!=rhs - r = self.solver_var(cp.boolvar()) - self.grb_model.update() - # r=1 -> lhs == rhs - # self.grb_model.addGenConstrIndicator(r, True, lhs_expr, GRB.EQUAL, 0) + return reify_cmp(GRB.EQUAL, lhs_expr, 0) + self.grb_model.addConstr((r==1) >> (lhs_expr == 0)) # r=0 -> lhs != rhs (reuse the != case) r_neq = add(lhs != rhs) self.grb_model.addConstr((r==0) >> (r_neq >= 1)) # self.grb_model.addConstr(r + r_neq >= 1) + # self.grb_model.addConstr(r == (grb_lhs == grb_rhs)) return r case "!=": # lhs != rhs => (lhs < rhs) | (lhs > rhs) @@ -577,9 +640,9 @@ def reify_cmp(sense, lhs_expr, rhs_const): return r case "<=": # lhs <= rhs => (lhs - rhs) <= 0 - return lhs_expr <= 0 + # return reify_cmp(lhs_expr <= 0) - # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) + return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) case ">=": # lhs >= rhs => (lhs - rhs) >= 0 return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 0) @@ -598,29 +661,32 @@ def reify_cmp(sense, lhs_expr, rhs_const): # transform and post the constraints for cpm_expr in self.transform(cpm_expr_orig): - print("E", cpm_expr) grb_expr = add(cpm_expr) self.grb_model.update() - print("out", grb_expr, type(grb_expr)) # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 if isinstance(grb_expr, (gp.Var, gp.LinExpr)): self.grb_model.addConstr(grb_expr >= 1) elif isinstance(grb_expr, gp.TempConstr): self.grb_model.addConstr(grb_expr) - elif type(grb_expr).__name__ == 'GenExprAnd': - # GenExprAnd needs addGenConstrAnd: r = AND(args) then r >= 1 - r = self.solver_var(cp.boolvar()) - self.grb_model.addGenConstrAnd(r, list(grb_expr)) - self.grb_model.addConstr(r >= 1) - elif type(grb_expr).__name__ == 'GenExprOr': - # GenExprOr needs addGenConstrOr: r = OR(args) then r >= 1 - r = self.solver_var(cp.boolvar()) - self.grb_model.addGenConstrOr(r, list(grb_expr)) - self.grb_model.addConstr(r >= 1) + # elif type(grb_expr).__name__ == 'GenExprAnd': + # # GenExprAnd needs addGenConstrAnd: r = AND(args) then r >= 1 + # r = self.solver_var(cp.boolvar()) + # self.grb_model.addGenConstrAnd(r, list(grb_expr)) + # self.grb_model.addConstr(r >= 1) + # elif type(grb_expr).__name__ == 'GenExprOr': + # # GenExprOr needs addGenConstrOr: r = OR(args) then r >= 1 + # r = self.solver_var(cp.boolvar()) + # self.grb_model.addGenConstrOr(r, list(grb_expr)) + # self.grb_model.addConstr(r >= 1) else: - r = self.solver_var(cp.boolvar()) + # r = self.solver_var(cp.boolvar()) + # self.grb_model.addConstr(r == grb_expr) + # self.grb_model.addConstr(r >= 1) + + r = self.grb_model.addVar(lb=1, ub=GRB.INFINITY) self.grb_model.addConstr(r == grb_expr) - self.grb_model.addConstr(r >= 1) + # self.grb_model.addGenConstrNL(r, grb_expr) + # self.grb_model.addConstr(1 == grb_expr) continue # Comparisons: only numeric ones as 'only_implies()' has removed the '==' reification for Boolean expressions diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 567f65e93..06c209861 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -120,7 +120,7 @@ def flatten_model(orig_model, csemap=None): return cp.Model(*basecons, maximize=newobj) -def flatten_constraint(expr, csemap=None): +def flatten_constraint(expr, csemap=None, supported={}): """ input is any expression; except is_num(), pure _NumVarImpl, or Operator/GlobalConstraint with not is_bool() @@ -143,6 +143,7 @@ def flatten_constraint(expr, csemap=None): lst_of_expr = simplify_boolean(lst_of_expr) # simplify boolean expressions, and ensure types are correct for expr in lst_of_expr: + print("expr", expr, supported) if not expr.has_subexpr(): newlist.append(expr) # no need to do anything continue @@ -155,7 +156,10 @@ def flatten_constraint(expr, csemap=None): Var -> Boolexpr (CPMpy class 'Operator', is_bool()) """ # does not type-check that arguments are bool... Could do now with expr.is_bool()! - if expr.name == 'or': + if expr.name in supported: + print("expr", expr, supported) + newlist.extend(flatten_constraint(expr.args, csemap=csemap)) + elif expr.name == 'or': # rewrites that avoid auxiliary var creation, should go to normalize? # in case of an implication in a disjunction, merge in if builtins.any(isinstance(a, Operator) and a.name == '->' for a in expr.args): @@ -264,7 +268,7 @@ def flatten_constraint(expr, csemap=None): # integer comparison (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap) else: - (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap) + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported) newlist.append(Comparison(exprname, lhs, rvar)) newlist.extend(lcons) @@ -300,7 +304,7 @@ def flatten_objective(expr, supported=frozenset(["sum", "wsum"]), csemap=None): raise Exception(f"Objective expects a single variable/expression, not a list of expressions: {expr}") expr = simplify_boolean([expr])[0] - (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap) # might rewrite expr into a (w)sum + (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap, supported=supported) # might rewrite expr into a (w)sum if isinstance(flatexpr, Expression) and flatexpr.name in supported: return (flatexpr, flatcons) else: @@ -322,7 +326,7 @@ def __is_flat_var_or_list(arg): is_any_list(arg) and all(__is_flat_var_or_list(el) for el in arg) or \ is_star(arg) -def get_or_make_var(expr, csemap=None): +def get_or_make_var(expr, csemap=None, supported={}): """ Must return a variable, and list of flat normal constraints Determines whether this is a Boolean or Integer variable and returns @@ -355,7 +359,7 @@ def get_or_make_var(expr, csemap=None): else: # normalize expr into a numexpr LHS, # then compute bounds and return (newintvar, LHS == newintvar) - (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap) + (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap, supported=supported) lb, ub = flatexpr.get_bounds() if not is_int(lb) or not is_int(ub): @@ -461,7 +465,7 @@ def normalized_boolexpr(expr, csemap=None): exprname = '==' else: # other cases: LHS is numexpr - (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap) + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported) return (Comparison(exprname, lhs, rvar), lcons+rcons) @@ -482,7 +486,7 @@ def normalized_boolexpr(expr, csemap=None): return (newexpr, [c for con in flatcons for c in con]) -def normalized_numexpr(expr, csemap=None): +def normalized_numexpr(expr, csemap=None, supported={}): """ all 'flat normal form' numeric expressions... @@ -511,12 +515,12 @@ def normalized_numexpr(expr, csemap=None): # rewrite const*a into a weighted sum, so it can be used as objective elif expr.name == "mul" and getattr(expr, "is_lhs_num", False): w, e = expr.args - return normalized_numexpr(Operator("wsum", ([w], [e])), csemap=csemap) + return normalized_numexpr(Operator("wsum", ([w], [e])), csemap=csemap, supported=supported) elif isinstance(expr, Operator): # rewrite -a into a weighted sum, so it can be used as objective if expr.name == '-': - return normalized_numexpr(Operator("wsum", _wsum_make(expr)), csemap=csemap) + return normalized_numexpr(Operator("wsum", _wsum_make(expr)), csemap=csemap, supported=supported) if not expr.has_subexpr(): return (expr, []) @@ -528,7 +532,7 @@ def normalized_numexpr(expr, csemap=None): we = [_wsum_make(a) for a in expr.args] w = [wi for w,_ in we for wi in w] e = [ei for _,e in we for ei in e] - return normalized_numexpr(Operator("wsum", (w,e)), csemap=csemap) + return normalized_numexpr(Operator("wsum", (w,e)), csemap=csemap, supported=supported) # wsum needs special handling because expr.args is a tuple of which only 2nd one has exprs if expr.name == 'wsum': @@ -555,15 +559,19 @@ def normalized_numexpr(expr, csemap=None): else: # generic operator # recursively flatten all children - flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap) for arg in expr.args]) + print("gen", expr.args, [a.name for a in expr.args]) + flatvars, flatcons = zip(*[(arg, []) if arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args]) + print("gen", flatvars) newexpr = Operator(expr.name, flatvars) return (newexpr, [c for con in flatcons for c in con]) else: # Globalfunction (examples: Max,Min,Element) + print("ge", expr) # just recursively flatten args, which can be lists if not expr.has_subexpr(): + print("b") return (expr, []) else: # recursively flatten all children From a6657a4ad565ebbddd90802ab4992f22e8d13826 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 31 Mar 2026 13:13:34 +0000 Subject: [PATCH 03/23] Update flatten --- cpmpy/solvers/gurobi.py | 102 ++++++++++++------------- cpmpy/transformations/flatten_model.py | 27 ++++--- 2 files changed, 61 insertions(+), 68 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 651dbc9c8..006170464 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -365,24 +365,25 @@ def transform(self, cpm_expr): # expressions have to be linearized to fit in MIP model. See /transformations/linearize cpm_cons = toplevel_list(cpm_expr) cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"}) # linearize and decompose expect safe exprs - cpm_cons = decompose_in_tree(cpm_cons, - supported=self.supported_global_constraints, - supported_reified=self.supported_reified_global_constraints, - csemap=self._csemap) + # cpm_cons = decompose_in_tree(cpm_cons, + # supported=self.supported_global_constraints, + # supported_reified=self.supported_reified_global_constraints, + # csemap=self._csemap) cpm_cons = decompose_linear(cpm_cons, supported=self.supported_global_constraints, supported_reified=self.supported_reified_global_constraints, csemap=self._csemap) - cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap, supported={"mul", "pow"}) # flat normal form + cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap, supported={"mul", "pow", "-", "sum"}) # flat normal form # cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification # cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != # cpm_cons = linearize_reified_variables(cpm_cons, min_values=2, csemap=self._csemap) - # cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) - # cpm_cons = only_implies(cpm_cons, csemap=self._csemap) # anything that can create full reif should go above... + cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) + cpm_cons = only_implies(cpm_cons, csemap=self._csemap) # anything that can create full reif should go above... # gurobi does not round towards zero, so no 'div' in supported set: https://github.com/CPMpy/cpmpy/pull/593#issuecomment-2786707188 - # cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization + cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"-","sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization + # cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv return cpm_cons @@ -504,7 +505,7 @@ def reify_cmp(sense, lhs_expr, rhs_const): # Base case: variables if isinstance(cpm_expr, NegBoolView): - return reify(1 - self.solver_var(cpm_expr._bv)) + return 1 - self.solver_var(cpm_expr._bv) if isinstance(cpm_expr, _NumVarImpl): return self.solver_var(cpm_expr) @@ -513,7 +514,7 @@ def reify_cmp(sense, lhs_expr, rhs_const): if isinstance(cpm_expr, Operator): match cpm_expr.name: - case "or": + case "or": # TODO only usefull if we can represent an or as a sum of binary vars, which is not generally correct # return reify(gp.or_(add(arg) for arg in cpm_expr.args)) self.grb_model.update() @@ -536,19 +537,11 @@ def reify_cmp(sense, lhs_expr, rhs_const): elif any(a == 0 for a in args if is_num(a)): return 0 else: - # return gp.and_(args) - # return sum(args) >= len(args) return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) - case "->": - # TODO has to be wrong ; reify(...) + case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) a, b = cpm_expr.args - lhs, rhs = add(a, depth), add(b, depth) - # return add(a, depth + 1) <= add(b, depth + 1) - self.grb_model.update() - # self.grb_model.addConstr((lhs==1) >> (rhs>=1)) - # return lhs - # return reify((lhs==1) >> (rhs>=1)) - return reify_cmp(GRB.LESS_EQUAL, lhs - rhs, 0) + is_pos = not isinstance(a, NegBoolView) + return (add(a if is_pos else a._bv, depth) == int(is_pos)) >> add(b, depth) case "not": # Negation: not(a) is 1 - a return 1 - add(cpm_expr.args[0], depth) @@ -556,25 +549,16 @@ def reify_cmp(sense, lhs_expr, rhs_const): return -add(cpm_expr.args[0], depth=depth) case "sum": return sum(add(arg, depth) for arg in cpm_expr.args) - return gp.quicksum(add(arg, depth) for arg in cpm_expr.args) case "wsum": weights, args = cpm_expr.args return sum(w * add(arg, depth) for w, arg in zip(weights, args)) - # return gp.quicksum(w * add(arg, depth) for w, arg in zip(weights, args)) case "sub": a, b = cpm_expr.args return add(a, depth) - add(b, depth) elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): match cpm_expr.name: case "pow": - # y = x^z - x, z = add(cpm_expr.args[0]), add(cpm_expr.args[1]) - # y = self.solver_var(cp.boolvar()) - # self.grb_model.addGenConstrPow(x, y, z) - # self.grb_model.addConstr(y == (x**z)) - # self.grb_model.addGenConstrNL(x, y, z) - # return reify(x**z) - return x**z + return add(cpm_expr.args[0]) ** add(cpm_expr.args[1]) case "mul": a, b = cpm_expr.args # return reify(add(a)*add(b)) @@ -611,16 +595,25 @@ def reify_cmp(sense, lhs_expr, rhs_const): # return i if isinstance(cpm_expr, Comparison): + # cpm_expr = linearize_constraint([cpm_expr], supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization + # print(cpm_expr) lhs, rhs = cpm_expr.args grb_lhs, grb_rhs = add(lhs, depth), add(rhs, depth) self.grb_model.update() # Ensure vars are in model before creating TempConstr # Rewrite to form: (grb_lhs - grb_rhs) sense constant # addGenConstrIndicator requires rhs to be a constant - lhs_expr = grb_lhs - grb_rhs # LinExpr + # lhs_expr = grb_lhs - grb_rhs # LinExpr match cpm_expr.name: case "==": + return grb_lhs == grb_rhs + case "<=": + return grb_lhs <= grb_rhs + case ">=": + return grb_lhs >= grb_rhs + case _: + assert False return reify_cmp(GRB.EQUAL, lhs_expr, 0) self.grb_model.addConstr((r==1) >> (lhs_expr == 0)) @@ -630,28 +623,28 @@ def reify_cmp(sense, lhs_expr, rhs_const): # self.grb_model.addConstr(r + r_neq >= 1) # self.grb_model.addConstr(r == (grb_lhs == grb_rhs)) return r - case "!=": - # lhs != rhs => (lhs < rhs) | (lhs > rhs) - r_lt = reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) # lhs - rhs <= -1 - r_gt = reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) # lhs - rhs >= 1 - r = self.solver_var(cp.boolvar()) - self.grb_model.update() - self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) - return r - case "<=": - # lhs <= rhs => (lhs - rhs) <= 0 - # return reify_cmp(lhs_expr <= 0) - - return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) - case ">=": - # lhs >= rhs => (lhs - rhs) >= 0 - return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 0) - case "<": - # lhs < rhs => (lhs - rhs) <= -1 - return reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) - case ">": - # lhs > rhs => (lhs - rhs) >= 1 - return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) + # case "!=": + # # lhs != rhs => (lhs < rhs) | (lhs > rhs) + # r_lt = reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) # lhs - rhs <= -1 + # r_gt = reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) # lhs - rhs >= 1 + # r = self.solver_var(cp.boolvar()) + # self.grb_model.update() + # self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) + # return r + # case "<=": + # # lhs <= rhs => (lhs - rhs) <= 0 + # # return reify_cmp(lhs_expr <= 0) + # + # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) + # case ">=": + # # lhs >= rhs => (lhs - rhs) >= 0 + # return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 0) + # case "<": + # # lhs < rhs => (lhs - rhs) <= -1 + # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) + # case ">": + # # lhs > rhs => (lhs - rhs) >= 1 + # return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) raise NotImplementedError(f"add() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") @@ -662,6 +655,7 @@ def reify_cmp(sense, lhs_expr, rhs_const): # transform and post the constraints for cpm_expr in self.transform(cpm_expr_orig): grb_expr = add(cpm_expr) + print("out", grb_expr) self.grb_model.update() # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 if isinstance(grb_expr, (gp.Var, gp.LinExpr)): diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 06c209861..8ef26f0db 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -198,11 +198,11 @@ def flatten_constraint(expr, csemap=None, supported={}): elif isinstance(expr.args[0], _BoolVarImpl): # LHS is var, ensure RHS is normalized 'Boolexpr' lhs,lcons = expr.args[0], () - rhs,rcons = normalized_boolexpr(expr.args[1], csemap=csemap) + rhs,rcons = normalized_boolexpr(expr.args[1], csemap=csemap, supported=supported) else: # make LHS normalized 'Boolexpr', RHS must be a var - lhs,lcons = normalized_boolexpr(expr.args[0], csemap=csemap) - rhs,rcons = get_or_make_var(expr.args[1], csemap=csemap) + lhs,lcons = normalized_boolexpr(expr.args[0], csemap=csemap, supported=supported) + rhs,rcons = get_or_make_var(expr.args[1], csemap=csemap, supported=supported) newlist.append(Operator(expr.name, (lhs,rhs))) newlist.extend(lcons) @@ -213,7 +213,7 @@ def flatten_constraint(expr, csemap=None, supported={}): # if none of the above cases + continue matched: # a normalizable boolexpr - (con, flatcons) = normalized_boolexpr(expr, csemap=csemap) + (con, flatcons) = normalized_boolexpr(expr, csemap=csemap, supported=supported) newlist.append(con) newlist.extend(flatcons) @@ -373,7 +373,7 @@ def get_or_make_var(expr, csemap=None, supported={}): csemap[expr] = ivar return ivar, [flatexpr == ivar] + flatcons -def get_or_make_var_or_list(expr, csemap=None): +def get_or_make_var_or_list(expr, csemap=None, supported={}): """ Like get_or_make_var() but also accepts and recursively transforms lists Used to convert arguments of globals """ @@ -381,13 +381,13 @@ def get_or_make_var_or_list(expr, csemap=None): if __is_flat_var_or_list(expr): return (expr,[]) elif is_any_list(expr): - flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap) for arg in expr]) + flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap, supported=supported) for arg in expr]) return (flatvars, [c for con in flatcons for c in con]) else: - return get_or_make_var(expr, csemap=csemap) + return get_or_make_var(expr, csemap=csemap, supported=supported) -def normalized_boolexpr(expr, csemap=None): +def normalized_boolexpr(expr, csemap=None, supported={}): """ input is any Boolean (is_bool()) expression output are all 'flat normal form' Boolean expressions that can be 'reified', meaning that @@ -503,6 +503,7 @@ def normalized_numexpr(expr, csemap=None, supported={}): - base_expr: same as 'expr', but all arguments are variables - base_cons: list of flat normal constraints """ + print("normalized_numexpr", expr.name, expr) # XXX a boolexpr is also a valid numexpr... e.g. 30*(iv > 5) + ... see mario obj. if __is_flat_var(expr): return (expr, []) @@ -559,23 +560,21 @@ def normalized_numexpr(expr, csemap=None, supported={}): else: # generic operator # recursively flatten all children - print("gen", expr.args, [a.name for a in expr.args]) - flatvars, flatcons = zip(*[(arg, []) if arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args]) - print("gen", flatvars) + flatvars, flatcons = zip(*(normalized_numexpr(arg, csemap=csemap, supported=supported) if arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args)) newexpr = Operator(expr.name, flatvars) return (newexpr, [c for con in flatcons for c in con]) else: # Globalfunction (examples: Max,Min,Element) - print("ge", expr) # just recursively flatten args, which can be lists if not expr.has_subexpr(): - print("b") return (expr, []) else: # recursively flatten all children - flatvars, flatcons = zip(*[get_or_make_var_or_list(arg, csemap=csemap) for arg in expr.args]) + flatvars, flatcons = zip(*[ + normalized_numexpr(arg, supported=supported) if arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap ,supported=supported) + for arg in expr.args]) # take copy, replace args newexpr = copy.copy(expr) # shallow or deep? currently shallow From cb13ee5f62100dd1eb488a6264037034616a9fae Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 31 Mar 2026 13:37:02 +0000 Subject: [PATCH 04/23] Clean-up / improve abs --- cpmpy/solvers/gurobi.py | 210 +++++++--------------------------------- 1 file changed, 35 insertions(+), 175 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 006170464..d3b7e704c 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -408,128 +408,45 @@ def add(self, cpm_expr_orig): from gurobipy import GRB, nlfunc import gurobipy as gp - def add(cpm_expr, depth=0): + def add(cpm_expr, depth): """Recursively create a Gurobi expression tree from a CPMpy expression.""" indent = " " * depth depth+=1 self.grb_model.update() print(f"{indent}Con:", cpm_expr, type(cpm_expr)) - # Base case: numbers + if is_num(cpm_expr): return int(cpm_expr) - # # Base case: variables - # if isinstance(cpm_expr, NegBoolView): - # assert False - # # return 1 - self.solver_var(cpm_expr) - def reify(grb_expr): - """Create a binary variable representing the truth value of cpm_expr.""" + """Create a binary variable representing the truth value of the gurobi logical constraint (OR/AND)""" self.grb_model.update() print(f"{indent}Reif", grb_expr) - - # # Base case: already a variable - # if isinstance(cpm_expr, _IntVarImpl) and not isinstance(cpm_expr, NegBoolView): - # return self.solver_var(cpm_expr) - # Base case: already a variable if isinstance(grb_expr, gp.Var): return grb_expr - r = self.solver_var(cp.boolvar()) - - - self.grb_model.update() - - # Double-reification: r == TempConstr works for linear constraints - self.grb_model.addConstr(r == grb_expr) - - - # print("R", r, grb_expr) - # if isinstance(grb_expr, gp.TempConstr): - # # self.grb_model.addConstr(r == (grb_expr)) - # self.grb_model.addConstr((r==1) >> (grb_expr)) - # else: - # # For other expressions, add() returns something we can equate - # # grb_expr = add(cpm_expr) - # self.grb_model.addConstr(r == grb_expr) - # return r - - - return r - - def reify_cmp(sense, lhs_expr, rhs_const): - """Reify a comparison using indicator constraints.""" - - # (2x _3y <= 5) - # (L == 2x + 3y, r == (L <= 5)) - - r_ = self.solver_var(cp.boolvar()) - - # r <-> (x = 3) - # L = x, r -> L = 3, ~r -> (L >= 4 \/ L <= 2) - # - - # TODO if not already LinExpr - L = self.grb_model.addVar(vtype=GRB.INTEGER, lb=-GRB.INFINITY,ub=GRB.INFINITY) - self.grb_model.addConstr(L == lhs_expr) r = self.solver_var(cp.boolvar()) - - self.grb_model.addGenConstrIndicator(r, True, L, sense, rhs_const) - # Add negation for double-reification - match sense: - case GRB.LESS_EQUAL: - self.grb_model.addGenConstrIndicator(r, False, L, GRB.GREATER_EQUAL, rhs_const + 1) - case GRB.GREATER_EQUAL: - self.grb_model.addGenConstrIndicator(r, False, L, GRB.LESS_EQUAL, rhs_const - 1) - case GRB.EQUAL: - # self.grb_model.addGenConstrIndicator(r, False, L, GRB.LESS_EQUAL, rhs_const - 1) - # D = self.grb_model.addVar(vtype=GRB.INTEGER, lb=-GRB.INFINITY,ub=GRB.INFINITY) - # self.grb_model.addConstr(D == gp.abs_(L - lhs_expr)) - # self.grb_model.addConstr((r==0) >> (D >= 1)) - - a = reify_cmp(GRB.GREATER_EQUAL, L, rhs_const + 1) - b = reify_cmp(GRB.LESS_EQUAL, L, rhs_const - 1) - print('neq', a,b) - self.grb_model.addConstr((r==0) >> (a + b >= 1)) - - # # self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) - # r_ = self.solver_var(cp.boolvar()) - # self.grb_model.addConstr((r==0) >> r_) - # self.grb_model.addConstr(r_ == ) + self.grb_model.addConstr(r == grb_expr) return r - # Base case: variables if isinstance(cpm_expr, NegBoolView): return 1 - self.solver_var(cpm_expr._bv) if isinstance(cpm_expr, _NumVarImpl): return self.solver_var(cpm_expr) - # def filter_(args): - # return gp.and_(add(arg) for arg in cpm_expr.args) - if isinstance(cpm_expr, Operator): match cpm_expr.name: case "or": # TODO only usefull if we can represent an or as a sum of binary vars, which is not generally correct - # return reify(gp.or_(add(arg) for arg in cpm_expr.args)) - self.grb_model.update() - - # return sum(add(arg, depth=depth) for arg in cpm_expr.args) - # def normalize(arg): - # if isinstance(arg, gp.QuadExpr): - # self.native_model - - # args = [normalize(arg) for arg in args] - return reify(gp.or_(add(arg) for arg in cpm_expr.args)) + return reify(gp.or_(add(arg, depth) for arg in cpm_expr.args)) case "and": - self.grb_model.update() - return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) + # self.grb_model.update() + # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) - # return reify(gp.and_([0,1])) args = [add(arg, depth=depth) for arg in cpm_expr.args if not is_num(arg) or arg != 1] assert len(args) if not args: @@ -537,9 +454,12 @@ def reify_cmp(sense, lhs_expr, rhs_const): elif any(a == 0 for a in args if is_num(a)): return 0 else: - return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) + return reify(gp.and_(add(arg) for arg in args)) + # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) a, b = cpm_expr.args + assert isinstance(a, _BoolVarImpl), f"Implication constraint {cpm_expr} must have BoolVar as lhs, but had {a}" + assert isinstance(b, Comparison), f"Implication constraint {cpm_expr} must have Comparison as rhs, but had {b}" is_pos = not isinstance(a, NegBoolView) return (add(a if is_pos else a._bv, depth) == int(is_pos)) >> add(b, depth) case "not": @@ -550,111 +470,51 @@ def reify_cmp(sense, lhs_expr, rhs_const): case "sum": return sum(add(arg, depth) for arg in cpm_expr.args) case "wsum": - weights, args = cpm_expr.args - return sum(w * add(arg, depth) for w, arg in zip(weights, args)) + return sum(weight * add(arg, depth) for weight, arg in zip(cpm_expr.args[0], cpm_expr.args[1])) case "sub": a, b = cpm_expr.args return add(a, depth) - add(b, depth) - elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): - match cpm_expr.name: - case "pow": - return add(cpm_expr.args[0]) ** add(cpm_expr.args[1]) - case "mul": - a, b = cpm_expr.args - # return reify(add(a)*add(b)) - return add(a) * add(b) - case _: - # flatargs, flatcons = zip(*[get_or_make_var_or_list(arg, csemap=self._csemap) for arg in cpm_expr.args]) - flatargs = [add(arg, depth) for arg in cpm_expr.args] - self.grb_model.update() - # self.add(flatcons) - match cpm_expr.name: - case "abs": # y = abs(x) - x, = flatargs - x_, = cpm_expr.args - # y = add(cp.intvar(lb=max(0,cpm_expr.args[0].lb), ub=cpm_expr.args[0].ub), depth) - y = add(cp.intvar(lb=x_.lb, ub=x_.ub), depth) - self.grb_model.update() - # self.native_model.addConstr(self.solver_var(y) == gp.abs_(add(x, depth))) - self.native_model.addConstr(y == gp.abs_(x)) - return y - - # return gp.abs_(add(cpm_expr.args[0], depth)) - # return reify(abs(add(cpm_expr.args[0], depth))) - case "min" | "max": - y = add(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) - self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) - return y - # elif isinstance(cpm_expr, cp.expressions.globalconstraints.GlobalConstraint): - # match cpm_expr.name: - # case "ite": - # i, t, e = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth), add(cpm_expr.args[2], depth) - # print('ite', i,t,e) - # self.grb_model.addConstr((i==1) >> (t>=1)) - # self.grb_model.addConstr((i==0) >> (e>=1)) - # return i - - if isinstance(cpm_expr, Comparison): + elif isinstance(cpm_expr, Comparison): # cpm_expr = linearize_constraint([cpm_expr], supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization # print(cpm_expr) lhs, rhs = cpm_expr.args grb_lhs, grb_rhs = add(lhs, depth), add(rhs, depth) - self.grb_model.update() # Ensure vars are in model before creating TempConstr - - # Rewrite to form: (grb_lhs - grb_rhs) sense constant - # addGenConstrIndicator requires rhs to be a constant - # lhs_expr = grb_lhs - grb_rhs # LinExpr match cpm_expr.name: case "==": - return grb_lhs == grb_rhs + # Note: rhs/lhs are reversed since Gurobi functions are called by `y == f(x)`, but CPMpy normalizes to `f(x) == y` (e.g. `abs(x) == IV0`) + return grb_rhs == grb_lhs case "<=": return grb_lhs <= grb_rhs case ">=": return grb_lhs >= grb_rhs case _: - assert False - return reify_cmp(GRB.EQUAL, lhs_expr, 0) - - self.grb_model.addConstr((r==1) >> (lhs_expr == 0)) - # r=0 -> lhs != rhs (reuse the != case) - r_neq = add(lhs != rhs) - self.grb_model.addConstr((r==0) >> (r_neq >= 1)) - # self.grb_model.addConstr(r + r_neq >= 1) - # self.grb_model.addConstr(r == (grb_lhs == grb_rhs)) - return r - # case "!=": - # # lhs != rhs => (lhs < rhs) | (lhs > rhs) - # r_lt = reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) # lhs - rhs <= -1 - # r_gt = reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) # lhs - rhs >= 1 - # r = self.solver_var(cp.boolvar()) - # self.grb_model.update() - # self.grb_model.addGenConstrOr(r, [r_lt, r_gt]) - # return r - # case "<=": - # # lhs <= rhs => (lhs - rhs) <= 0 - # # return reify_cmp(lhs_expr <= 0) - # - # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, 0) - # case ">=": - # # lhs >= rhs => (lhs - rhs) >= 0 - # return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 0) - # case "<": - # # lhs < rhs => (lhs - rhs) <= -1 - # return reify_cmp(GRB.LESS_EQUAL, lhs_expr, -1) - # case ">": - # # lhs > rhs => (lhs - rhs) >= 1 - # return reify_cmp(GRB.GREATER_EQUAL, lhs_expr, 1) - - raise NotImplementedError(f"add() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") - + raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") + elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): + match cpm_expr.name: + case "mul": + return add(cpm_expr.args[0], depth) * add(cpm_expr.args[1], depth) + case "pow": + return add(cpm_expr.args[0], depth) ** add(cpm_expr.args[1], depth) + case _: # other global functions cannot be + # assert isinstance(cpm_expr, Comparison) and cpm_expr.name == "==", f"Expected global function to be transformed to f(x) == y, but was {cpm_expr}" + flatargs = [add(arg, depth) for arg in cpm_expr.args] + match cpm_expr.name: + case "abs": # y = abs(x) + return gp.abs_(add(cpm_expr.args[0], depth)) + case "min" | "max": + y = add(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) + self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) + return y + else: + raise NotImplementedError(f"add() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") # add new user vars to the set get_variables(cpm_expr_orig, collect=self.user_vars) # transform and post the constraints for cpm_expr in self.transform(cpm_expr_orig): - grb_expr = add(cpm_expr) + grb_expr = add(cpm_expr, 0) print("out", grb_expr) self.grb_model.update() # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 From cdf1f880cc5ff5aeedeb7571de5410a17bc10f35 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 31 Mar 2026 15:37:07 +0000 Subject: [PATCH 05/23] Fix conjunction/disjunction --- cpmpy/solvers/gurobi.py | 26 ++++++++++++++------------ cpmpy/transformations/flatten_model.py | 2 +- cpmpy/transformations/linearize.py | 2 +- cpmpy/transformations/reification.py | 7 +++++-- 4 files changed, 21 insertions(+), 16 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index d3b7e704c..8b1d02f49 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -379,10 +379,14 @@ def transform(self, cpm_expr): # cpm_cons = linearize_reified_variables(cpm_cons, min_values=2, csemap=self._csemap) cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) - cpm_cons = only_implies(cpm_cons, csemap=self._csemap) # anything that can create full reif should go above... + cpm_cons = only_implies( + cpm_cons, + csemap=self._csemap, + is_supported=lambda cpm_expr: cpm_expr.name == "==" and (isinstance(cpm_expr.args[1], Operator) and cpm_expr.args[1].name in {"or", "and"}) + ) # anything that can create full reif should go above... # gurobi does not round towards zero, so no 'div' in supported set: https://github.com/CPMpy/cpmpy/pull/593#issuecomment-2786707188 - cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"-","sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization + cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"-","sum", "wsum","->","sub","min","max","mul","abs","pow", "or", "and"}), csemap=self._csemap) # the core of the MIP-linearization # cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv return cpm_cons @@ -442,11 +446,12 @@ def reify(grb_expr): if isinstance(cpm_expr, Operator): match cpm_expr.name: case "or": # TODO only usefull if we can represent an or as a sum of binary vars, which is not generally correct - return reify(gp.or_(add(arg, depth) for arg in cpm_expr.args)) + return gp.or_(add(arg, depth) for arg in cpm_expr.args) case "and": # self.grb_model.update() # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) + return gp.and_(add(arg, depth) for arg in cpm_expr.args) args = [add(arg, depth=depth) for arg in cpm_expr.args if not is_num(arg) or arg != 1] assert len(args) if not args: @@ -454,7 +459,7 @@ def reify(grb_expr): elif any(a == 0 for a in args if is_num(a)): return 0 else: - return reify(gp.and_(add(arg) for arg in args)) + return gp.and_(add(arg) for arg in args) # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) a, b = cpm_expr.args @@ -475,19 +480,16 @@ def reify(grb_expr): a, b = cpm_expr.args return add(a, depth) - add(b, depth) elif isinstance(cpm_expr, Comparison): - # cpm_expr = linearize_constraint([cpm_expr], supported=frozenset({"sum", "wsum","->","sub","min","max","mul","abs","pow"}), csemap=self._csemap) # the core of the MIP-linearization - # print(cpm_expr) - lhs, rhs = cpm_expr.args - grb_lhs, grb_rhs = add(lhs, depth), add(rhs, depth) + a, b = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth) match cpm_expr.name: case "==": - # Note: rhs/lhs are reversed since Gurobi functions are called by `y == f(x)`, but CPMpy normalizes to `f(x) == y` (e.g. `abs(x) == IV0`) - return grb_rhs == grb_lhs + # Note: Gurobi functions are called by `y == f(x)`, like CPMpy boolexprs, but numexprs are transformed to `f(x) == y` (e.g. `abs(x) == IV0`), so need to be flipped + return a == b if cpm_expr.args[0].is_bool() else b == a case "<=": - return grb_lhs <= grb_rhs + return a <= b case ">=": - return grb_lhs >= grb_rhs + return a >= b case _: raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 8ef26f0db..25f7ddf12 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -573,7 +573,7 @@ def normalized_numexpr(expr, csemap=None, supported={}): else: # recursively flatten all children flatvars, flatcons = zip(*[ - normalized_numexpr(arg, supported=supported) if arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap ,supported=supported) + normalized_numexpr(arg, supported=supported) if arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap, supported=supported) for arg in expr.args]) # take copy, replace args diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py index 005e28bc9..c40f0a7b8 100644 --- a/cpmpy/transformations/linearize.py +++ b/cpmpy/transformations/linearize.py @@ -200,7 +200,7 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal assert len(new_expr) == 1 if isinstance(new_expr[0], BoolVal) and new_expr[0].value() is True: continue # skip or([BoolVal(True)]) - newlist.append(Operator("or", new_expr)) + newlist.extend(new_expr) continue diff --git a/cpmpy/transformations/reification.py b/cpmpy/transformations/reification.py index b3d682fad..f83337675 100644 --- a/cpmpy/transformations/reification.py +++ b/cpmpy/transformations/reification.py @@ -30,6 +30,7 @@ from .negation import recurse_negation def only_bv_reifies(constraints, csemap=None): + """Transforms all reifications to ``BV -> BE`` or ``BV == BE``""" newcons = [] for cpm_expr in constraints: @@ -52,7 +53,7 @@ def only_bv_reifies(constraints, csemap=None): newcons.append(cpm_expr) return newcons -def only_implies(constraints, csemap=None): +def only_implies(constraints, csemap=None, is_supported=None): """ Transforms all reifications to ``BV -> BE`` form @@ -73,7 +74,9 @@ def only_implies(constraints, csemap=None): for cpm_expr in constraints: # Operators: check BE -> BV - if cpm_expr.name == '->' and cpm_expr.args[1].name == '==': + if is_supported and is_supported(cpm_expr): + newcons.append(cpm_expr) + elif cpm_expr.name == '->' and cpm_expr.args[1].name == '==': a0,a1 = cpm_expr.args if a1.args[0].is_bool() and a1.args[1].is_bool(): # BV0 -> BV2 == BV3 :: BV0 -> (BV2->BV3 & BV3->BV2) From 8b552189209f9c476c2346d8c9f546acd7f97d92 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 12:29:13 +0000 Subject: [PATCH 06/23] Fix test_constraints except for global_constraints --- cpmpy/solvers/gurobi.py | 28 +++++---- cpmpy/transformations/flatten_model.py | 14 +++-- cpmpy/transformations/linearize.py | 84 +++++++++++++++++--------- cpmpy/transformations/reification.py | 14 +++-- 4 files changed, 93 insertions(+), 47 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 8b1d02f49..5aa4874f0 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -364,31 +364,37 @@ def transform(self, cpm_expr): # apply transformations, then post internally # expressions have to be linearized to fit in MIP model. See /transformations/linearize cpm_cons = toplevel_list(cpm_expr) - cpm_cons = no_partial_functions(cpm_cons, safen_toplevel={"mod", "div", "element"}) # linearize and decompose expect safe exprs - # cpm_cons = decompose_in_tree(cpm_cons, - # supported=self.supported_global_constraints, - # supported_reified=self.supported_reified_global_constraints, - # csemap=self._csemap) + cpm_cons = no_partial_functions(cpm_cons, safen_toplevel=frozenset(["mod", "div", "element"])) # linearize and decompose expect safe exprs cpm_cons = decompose_linear(cpm_cons, supported=self.supported_global_constraints, supported_reified=self.supported_reified_global_constraints, csemap=self._csemap) - cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap, supported={"mul", "pow", "-", "sum"}) # flat normal form - # cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification - # cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != + cpm_cons = flatten_constraint(cpm_cons, csemap=self._csemap, supported=frozenset(["mul", "pow", "-", "sum"])) # flat normal form + cpm_cons = reify_rewrite(cpm_cons, supported=frozenset(['sum', 'wsum']), csemap=self._csemap) # constraints that support reification + cpm_cons = only_numexpr_equality(cpm_cons, supported=frozenset(["sum", "wsum", "sub"]), csemap=self._csemap) # supports >, <, != # cpm_cons = linearize_reified_variables(cpm_cons, min_values=2, csemap=self._csemap) cpm_cons = only_bv_reifies(cpm_cons, csemap=self._csemap) cpm_cons = only_implies( cpm_cons, csemap=self._csemap, - is_supported=lambda cpm_expr: cpm_expr.name == "==" and (isinstance(cpm_expr.args[1], Operator) and cpm_expr.args[1].name in {"or", "and"}) + is_supported=lambda cpm_expr: + ( + cpm_expr.name == "==" and + ( + (isinstance(cpm_expr.args[1], Operator) and cpm_expr.args[1].name in {"or", "and"}) or + isinstance(cpm_expr.args[1], _BoolVarImpl) + ) + ) or ( + cpm_expr.name == "->" and isinstance(cpm_expr.args[1], Comparison) + ) ) # anything that can create full reif should go above... # gurobi does not round towards zero, so no 'div' in supported set: https://github.com/CPMpy/cpmpy/pull/593#issuecomment-2786707188 - cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"-","sum", "wsum","->","sub","min","max","mul","abs","pow", "or", "and"}), csemap=self._csemap) # the core of the MIP-linearization + cpm_cons = linearize_constraint(cpm_cons, supported=frozenset({"-","sum", "wsum","->","sub","min","max","mul","abs","pow","and","or"}), csemap=self._csemap) # the core of the MIP-linearization - # cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv + cpm_cons = only_positive_bv(cpm_cons, csemap=self._csemap) # after linearization, rewrite ~bv into 1-bv + # TODO don't rewrite ~p -> LinCons return cpm_cons def add(self, cpm_expr_orig): diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 25f7ddf12..122a84679 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -263,10 +263,12 @@ def flatten_constraint(expr, csemap=None, supported={}): if exprname == '==' and lexpr.is_bool(): if rvar.is_bool(): # this is a reification - (lhs, lcons) = normalized_boolexpr(lexpr, csemap=csemap) + (lhs, lcons) = normalized_boolexpr(lexpr, csemap=csemap, supported=supported) else: # integer comparison - (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap) + (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap, supported=supported) + elif expr.name == "!=": # TODO ; this risks making a reified LinCons with a QuadExpr + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap) else: (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported) @@ -503,7 +505,6 @@ def normalized_numexpr(expr, csemap=None, supported={}): - base_expr: same as 'expr', but all arguments are variables - base_cons: list of flat normal constraints """ - print("normalized_numexpr", expr.name, expr) # XXX a boolexpr is also a valid numexpr... e.g. 30*(iv > 5) + ... see mario obj. if __is_flat_var(expr): return (expr, []) @@ -560,7 +561,8 @@ def normalized_numexpr(expr, csemap=None, supported={}): else: # generic operator # recursively flatten all children - flatvars, flatcons = zip(*(normalized_numexpr(arg, csemap=csemap, supported=supported) if arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args)) + # TODO should be some isinstance _IntVarImpl? + flatvars, flatcons = zip(*(normalized_numexpr(arg, csemap=csemap, supported=supported) if hasattr(arg, "name") and arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args)) newexpr = Operator(expr.name, flatvars) return (newexpr, [c for con in flatcons for c in con]) @@ -573,7 +575,9 @@ def normalized_numexpr(expr, csemap=None, supported={}): else: # recursively flatten all children flatvars, flatcons = zip(*[ - normalized_numexpr(arg, supported=supported) if arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap, supported=supported) + normalized_numexpr(arg, supported=supported) + if arg.name in supported and expr.name + else get_or_make_var_or_list(arg, csemap=csemap, supported=supported) for arg in expr.args]) # take copy, replace args diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py index c40f0a7b8..1554b82b5 100644 --- a/cpmpy/transformations/linearize.py +++ b/cpmpy/transformations/linearize.py @@ -64,7 +64,8 @@ from ..expressions.core import Comparison, Expression, Operator, BoolVal from ..expressions.globalconstraints import GlobalConstraint, DirectConstraint, AllDifferent, Table, NegativeTable from ..expressions.globalfunctions import GlobalFunction, Element, NValue, Count -from ..expressions.utils import is_bool, is_num, is_int, eval_comparison, get_bounds, is_true_cst, is_false_cst +from ..expressions.utils import is_bool, is_num, is_int, eval_comparison, get_bounds, is_true_cst, is_false_cst, is_boolexpr + from ..expressions.variables import _BoolVarImpl, boolvar, NegBoolView, _NumVarImpl from .int2bool import _encode_int_var @@ -89,7 +90,7 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal if isinstance(cpm_expr, _BoolVarImpl): if "or" in supported: # post clause explicitly (don't use cp.any, which will just return the BoolVar) - newlist.append(Operator("or", [cpm_expr])) + newlist.append(cpm_expr) elif isinstance(cpm_expr, NegBoolView): # might as well remove the negation newlist.append(sum([~cpm_expr]) <= 0) @@ -115,8 +116,11 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal f"calling `linearize_constraint`" if isinstance(cond, _BoolVarImpl) and isinstance(sub_expr, _BoolVarImpl): - # shortcut for BV -> BV, convert to disjunction and apply linearize on it - newlist.append(1 * cond + -1 * sub_expr <= 0) + # shortcut for BV -> BV, convert to disjunction and linearize it (if unsupported) + if "->" in supported: + newlist.append(cond.implies(sub_expr >= 1)) + else: + newlist.append(1 * cond + -1 * sub_expr <= 0) # BV -> LinExpr elif isinstance(cond, _BoolVarImpl): @@ -270,6 +274,24 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal return newlist +def only_positive_bv_args(args, csemap=None): + # other operators in comparison such as "min", "max" + nbv_sel = [isinstance(a, NegBoolView) for a in args] + new_cons = [] + if any(nbv_sel): + new_args = [] + for i, nbv in enumerate(nbv_sel): + if nbv: + aux = cp.boolvar() # TODO not added to CSE?? + new_args.append(aux) + new_cons += [aux + args[i]._bv == 1] # aux == 1 - arg._bv + else: + new_args.append(args[i]) + + return new_args, new_cons + else: + return args, [] + def only_positive_bv(lst_of_expr, csemap=None): """ Replaces :class:`~cpmpy.expressions.comparison.Comparison` containing :class:`~cpmpy.expressions.variables.NegBoolView` with equivalent expression using only :class:`~cpmpy.expressions.variables.BoolVar`. @@ -280,32 +302,42 @@ def only_positive_bv(lst_of_expr, csemap=None): newlist = [] for cpm_expr in lst_of_expr: - if isinstance(cpm_expr, Comparison): + if isinstance(cpm_expr, Operator) and cpm_expr.name in {"or", "and"}: + new_args, new_cons = only_positive_bv_args(cpm_expr.args) + cpm_expr_ = copy.copy(cpm_expr) + cpm_expr_.update_args(new_args) + newlist.extend([cpm_expr_] + new_cons) + elif isinstance(cpm_expr, Comparison): lhs, rhs = cpm_expr.args new_lhs = lhs - new_cons = [] - if isinstance(lhs, _NumVarImpl) or lhs.name in {"sum","wsum"}: + new_cons = [] + if (isinstance(lhs, _NumVarImpl) or lhs.name in {"sum","wsum"}) and not is_boolexpr(rhs): new_lhs, const = only_positive_bv_wsum_const(lhs) rhs -= const + elif isinstance(lhs, _BoolVarImpl): + (new_lhs,), new_cons_ = only_positive_bv_args([lhs]) + new_cons += new_cons_ + else: + new_args, new_cons_ = only_positive_bv_args(lhs.args) + new_lhs = copy.copy(lhs) + new_lhs.update_args(new_args) + new_cons += new_cons_ + + if isinstance(rhs, _BoolVarImpl): + (new_rhs,), new_cons_ = only_positive_bv_args([rhs]) + new_cons += new_cons_ + elif is_boolexpr(rhs): + new_args, new_cons_ = only_positive_bv_args(rhs.args) + new_rhs = copy.copy(rhs) + new_rhs.update_args(new_args) + new_cons += new_cons_ + print("NEW", new_rhs) else: - # other operators in comparison such as "min", "max" - nbv_sel = [isinstance(a, NegBoolView) for a in lhs.args] - if any(nbv_sel): - new_args = [] - for i, nbv in enumerate(nbv_sel): - if nbv: - aux = cp.boolvar() - new_args.append(aux) - new_cons += [aux + lhs.args[i]._bv == 1] # aux == 1 - arg._bv - else: - new_args.append(lhs.args[i]) - - new_lhs = copy.copy(lhs) - new_lhs.update_args(new_args) - - if new_lhs is not lhs: - newlist.append(eval_comparison(cpm_expr.name, new_lhs, rhs)) + new_rhs = rhs + + if new_lhs is not lhs or new_rhs is not rhs: + newlist.append(eval_comparison(cpm_expr.name, new_lhs, new_rhs)) newlist += new_cons # already linear else: newlist.append(cpm_expr) @@ -319,9 +351,7 @@ def only_positive_bv(lst_of_expr, csemap=None): subexpr = only_positive_bv([subexpr], csemap=csemap) newlist += [cond.implies(expr) for expr in subexpr] - elif isinstance(cpm_expr, _BoolVarImpl): - raise ValueError(f"Unreachable: unexpected Boolean literal (`_BoolVarImpl`) in expression {cpm_expr}, perhaps `linearize_constraint` was not called before this `only_positive_bv `call") - elif isinstance(cpm_expr, (GlobalConstraint, BoolVal, DirectConstraint)): + elif isinstance(cpm_expr, (GlobalConstraint, BoolVal, DirectConstraint, _BoolVarImpl)): newlist.append(cpm_expr) else: raise Exception(f"{cpm_expr} is not linear or is not supported. Please report on github") diff --git a/cpmpy/transformations/reification.py b/cpmpy/transformations/reification.py index f83337675..1179f4c77 100644 --- a/cpmpy/transformations/reification.py +++ b/cpmpy/transformations/reification.py @@ -175,10 +175,16 @@ def reify_rewrite(constraints, supported=frozenset(), csemap=None): else: # reification, check for rewrite boolexpr = cpm_expr.args[boolexpr_index] if isinstance(boolexpr, Operator): - # Case 1, BE is Operator (and, or, ->) - # assume supported, return as is - newcons.append(cpm_expr) - # could actually rewrite into list of clauses like to_cnf() does... not for here + if boolexpr.name in supported or cpm_expr.name == "==": + # Case 1, BE is Operator (and, or, ->) + newcons.append(cpm_expr) + # could actually rewrite into list of clauses like to_cnf() does... not for here + elif cpm_expr.name == "->": + (auxvar, cons) = get_or_make_var(boolexpr, csemap=csemap) + newcons += cons + reifexpr = copy.copy(cpm_expr) + reifexpr.args[boolexpr_index] = auxvar + newcons.append(reifexpr) elif isinstance(boolexpr, GlobalConstraint): # Case 2, BE is a GlobalConstraint # replace BE by its decomposition, then flatten From fed45a47335bd290bcb661ed175a8f3ac7830db5 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 13:02:13 +0000 Subject: [PATCH 07/23] Fix XOR issue --- cpmpy/transformations/linearize.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py index 1554b82b5..80ff8ca24 100644 --- a/cpmpy/transformations/linearize.py +++ b/cpmpy/transformations/linearize.py @@ -88,7 +88,8 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal for cpm_expr in lst_of_expr: # Boolean literals are handled as trivial linears or unit clauses depending on `supported` if isinstance(cpm_expr, _BoolVarImpl): - if "or" in supported: + # TODO gurobi specifically cannot do reified or's + if "or" in supported and not reified: # post clause explicitly (don't use cp.any, which will just return the BoolVar) newlist.append(cpm_expr) elif isinstance(cpm_expr, NegBoolView): @@ -176,6 +177,11 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal elif isinstance(cpm_expr, Comparison): lhs, rhs = cpm_expr.args + # BV == and([a,b,c]) + if cpm_expr.name == "==" and hasattr(rhs, "name") and rhs.name in supported: + newlist.append(cpm_expr) + continue + if lhs.name == "sub": # convert to wsum lhs = Operator("wsum", [[1, -1], [lhs.args[0], lhs.args[1]]]) @@ -201,6 +207,7 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal if lhs.name == "sum" and len(lhs.args) == 1 and isinstance(lhs.args[0], _BoolVarImpl) and "or" in supported: # very special case, avoid writing as sum of 1 argument new_expr = simplify_boolean([eval_comparison(cpm_expr.name,lhs.args[0], rhs)]) + new_expr = linearize_constraint(new_expr, supported=supported, reified=reified, csemap=csemap) assert len(new_expr) == 1 if isinstance(new_expr[0], BoolVal) and new_expr[0].value() is True: continue # skip or([BoolVal(True)]) From dfac782c202e490160cb7d8608827a76ff6cabbb Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 13:34:55 +0000 Subject: [PATCH 08/23] Fix Circuit issue --- cpmpy/transformations/flatten_model.py | 84 +++++++++++++------------- 1 file changed, 43 insertions(+), 41 deletions(-) diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 122a84679..336f18467 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -120,7 +120,7 @@ def flatten_model(orig_model, csemap=None): return cp.Model(*basecons, maximize=newobj) -def flatten_constraint(expr, csemap=None, supported={}): +def flatten_constraint(expr, csemap=None, supported={}, reified=False): """ input is any expression; except is_num(), pure _NumVarImpl, or Operator/GlobalConstraint with not is_bool() @@ -143,7 +143,6 @@ def flatten_constraint(expr, csemap=None, supported={}): lst_of_expr = simplify_boolean(lst_of_expr) # simplify boolean expressions, and ensure types are correct for expr in lst_of_expr: - print("expr", expr, supported) if not expr.has_subexpr(): newlist.append(expr) # no need to do anything continue @@ -157,8 +156,7 @@ def flatten_constraint(expr, csemap=None, supported={}): """ # does not type-check that arguments are bool... Could do now with expr.is_bool()! if expr.name in supported: - print("expr", expr, supported) - newlist.extend(flatten_constraint(expr.args, csemap=csemap)) + newlist.extend(flatten_constraint(expr.args, csemap=csemap, reified=reified)) elif expr.name == 'or': # rewrites that avoid auxiliary var creation, should go to normalize? # in case of an implication in a disjunction, merge in @@ -168,7 +166,7 @@ def flatten_constraint(expr, csemap=None, supported={}): if isinstance(a, Operator) and a.name == '->': newargs[i:i+1] = [~a.args[0],a.args[1]] # there could be nested implications - newlist.extend(flatten_constraint(Operator('or', newargs), csemap=csemap)) + newlist.extend(flatten_constraint(Operator('or', newargs), csemap=csemap, reified=reified)) continue # conjunctions in disjunctions could be split out by applying distributivity, # but this would explode the number of constraints in favour of having less auxiliary variables. @@ -179,30 +177,30 @@ def flatten_constraint(expr, csemap=None, supported={}): if expr.args[1].name == 'and': a1s = expr.args[1].args a0 = expr.args[0] - newlist.extend(flatten_constraint([a0.implies(a1) for a1 in a1s], csemap=csemap)) + newlist.extend(flatten_constraint([a0.implies(a1) for a1 in a1s], csemap=csemap, reified=True)) continue # 2) if lhs is 'or' then or([a01..a0n])->a1 :: ~a1->and([~a01..~a0n] and split elif expr.args[0].name == 'or': a0s = expr.args[0].args a1 = expr.args[1] - newlist.extend(flatten_constraint([(~a1).implies(~a0) for a0 in a0s], csemap=csemap)) + newlist.extend(flatten_constraint([(~a1).implies(~a0) for a0 in a0s], csemap=csemap, reified=True)) continue # 2b) if lhs is ->, like 'or': a01->a02->a1 :: (~a01|a02)->a1 :: ~a1->a01,~a1->~a02 elif expr.args[0].name == '->': a01,a02 = expr.args[0].args a1 = expr.args[1] - newlist.extend(flatten_constraint([(~a1).implies(a01), (~a1).implies(~a02)], csemap=csemap)) + newlist.extend(flatten_constraint([(~a1).implies(a01), (~a1).implies(~a02)], csemap=csemap, reified=True)) continue # ->, allows a boolexpr on one side elif isinstance(expr.args[0], _BoolVarImpl): # LHS is var, ensure RHS is normalized 'Boolexpr' lhs,lcons = expr.args[0], () - rhs,rcons = normalized_boolexpr(expr.args[1], csemap=csemap, supported=supported) + rhs,rcons = normalized_boolexpr(expr.args[1], csemap=csemap, supported=supported, reified=True) else: # make LHS normalized 'Boolexpr', RHS must be a var - lhs,lcons = normalized_boolexpr(expr.args[0], csemap=csemap, supported=supported) - rhs,rcons = get_or_make_var(expr.args[1], csemap=csemap, supported=supported) + lhs,lcons = normalized_boolexpr(expr.args[0], csemap=csemap, supported=supported, reified=True) + rhs,rcons = get_or_make_var(expr.args[1], csemap=csemap, supported=supported, reified=True) newlist.append(Operator(expr.name, (lhs,rhs))) newlist.extend(lcons) @@ -213,7 +211,7 @@ def flatten_constraint(expr, csemap=None, supported={}): # if none of the above cases + continue matched: # a normalizable boolexpr - (con, flatcons) = normalized_boolexpr(expr, csemap=csemap, supported=supported) + (con, flatcons) = normalized_boolexpr(expr, csemap=csemap, supported=supported, reified=reified) newlist.append(con) newlist.extend(flatcons) @@ -257,20 +255,20 @@ def flatten_constraint(expr, csemap=None, supported={}): continue # ensure rhs is var - (rvar, rcons) = get_or_make_var(rexpr, csemap=csemap) + (rvar, rcons) = get_or_make_var(rexpr, csemap=csemap, reified=reified) # Reification (double implication): Boolexpr == Var # normalize the lhs (does not have to be a var, hence we call normalize instead of get_or_make_var if exprname == '==' and lexpr.is_bool(): if rvar.is_bool(): # this is a reification - (lhs, lcons) = normalized_boolexpr(lexpr, csemap=csemap, supported=supported) + (lhs, lcons) = normalized_boolexpr(lexpr, csemap=csemap, supported=supported, reified=reified) else: # integer comparison - (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap, supported=supported) + (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap, supported=supported, reified=reified) elif expr.name == "!=": # TODO ; this risks making a reified LinCons with a QuadExpr - (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap) + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, reified=reified) else: - (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported) + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported, reified=reified) newlist.append(Comparison(exprname, lhs, rvar)) newlist.extend(lcons) @@ -311,7 +309,7 @@ def flatten_objective(expr, supported=frozenset(["sum", "wsum"]), csemap=None): return (flatexpr, flatcons) else: # any other numeric expression, - var, cons = get_or_make_var(flatexpr, csemap=csemap) + var, cons = get_or_make_var(flatexpr, csemap=csemap, reified=reified) return (var, cons+flatcons) @@ -328,7 +326,7 @@ def __is_flat_var_or_list(arg): is_any_list(arg) and all(__is_flat_var_or_list(el) for el in arg) or \ is_star(arg) -def get_or_make_var(expr, csemap=None, supported={}): +def get_or_make_var(expr, csemap=None, supported={}, reified=False): """ Must return a variable, and list of flat normal constraints Determines whether this is a Boolean or Integer variable and returns @@ -346,7 +344,7 @@ def get_or_make_var(expr, csemap=None, supported={}): if expr.is_bool(): # normalize expr into a boolexpr LHS, reify LHS == bvar - (flatexpr, flatcons) = normalized_boolexpr(expr, csemap=csemap) + (flatexpr, flatcons) = normalized_boolexpr(expr, csemap=csemap, reified=reified) if isinstance(flatexpr,_BoolVarImpl): # avoids unnecessary bv == bv or bv == ~bv assignments @@ -361,7 +359,7 @@ def get_or_make_var(expr, csemap=None, supported={}): else: # normalize expr into a numexpr LHS, # then compute bounds and return (newintvar, LHS == newintvar) - (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap, supported=supported) + (flatexpr, flatcons) = normalized_numexpr(expr, csemap=csemap, supported=supported, reified=reified) lb, ub = flatexpr.get_bounds() if not is_int(lb) or not is_int(ub): @@ -375,7 +373,7 @@ def get_or_make_var(expr, csemap=None, supported={}): csemap[expr] = ivar return ivar, [flatexpr == ivar] + flatcons -def get_or_make_var_or_list(expr, csemap=None, supported={}): +def get_or_make_var_or_list(expr, csemap=None, supported={}, reified=False): """ Like get_or_make_var() but also accepts and recursively transforms lists Used to convert arguments of globals """ @@ -383,13 +381,13 @@ def get_or_make_var_or_list(expr, csemap=None, supported={}): if __is_flat_var_or_list(expr): return (expr,[]) elif is_any_list(expr): - flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap, supported=supported) for arg in expr]) + flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap, supported=supported, reified=reified) for arg in expr]) return (flatvars, [c for con in flatcons for c in con]) else: - return get_or_make_var(expr, csemap=csemap, supported=supported) + return get_or_make_var(expr, csemap=csemap, supported=supported, reified=reified) -def normalized_boolexpr(expr, csemap=None, supported={}): +def normalized_boolexpr(expr, csemap=None, supported={}, reified=False): """ input is any Boolean (is_bool()) expression output are all 'flat normal form' Boolean expressions that can be 'reified', meaning that @@ -425,18 +423,18 @@ def normalized_boolexpr(expr, csemap=None, supported={}): # apply De Morgan's transform for "implies" if expr.name == '->': # TODO, optimisation if args0 is an 'and'? - (lhs,lcons) = get_or_make_var(expr.args[0], csemap=csemap) + (lhs,lcons) = get_or_make_var(expr.args[0], csemap=csemap, reified=reified) # TODO, optimisation if args1 is an 'or'? - (rhs,rcons) = get_or_make_var(expr.args[1], csemap=csemap) + (rhs,rcons) = get_or_make_var(expr.args[1], csemap=csemap, reified=reified) return ((~lhs | rhs), lcons+rcons) if expr.name == 'not': - flatvar, flatcons = get_or_make_var(expr.args[0], csemap=csemap) + flatvar, flatcons = get_or_make_var(expr.args[0], csemap=csemap, reified=reified) return (~flatvar, flatcons) if not expr.has_subexpr(): return (expr, []) else: # one of the arguments is not flat, flatten all - flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap) for arg in expr.args]) + flatvars, flatcons = zip(*[get_or_make_var(arg, csemap=csemap, reified=reified) for arg in expr.args]) newexpr = Operator(expr.name, flatvars) return (newexpr, [c for con in flatcons for c in con]) @@ -455,19 +453,19 @@ def normalized_boolexpr(expr, csemap=None, supported={}): lexpr, rexpr = rexpr, lexpr # ensure rhs is var - (rvar, rcons) = get_or_make_var(rexpr, csemap=csemap) + (rvar, rcons) = get_or_make_var(rexpr, csemap=csemap, reified=reified) # LHS: check if Boolexpr == smth: if (exprname == '==' or exprname == '!=') and lexpr.is_bool(): # this is a reified constraint, so lhs must be var too to be in normal form - (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap) + (lhs, lcons) = get_or_make_var(lexpr, csemap=csemap, reified=reified) if expr.name == '!=' and rvar.is_bool(): # != not needed, negate RHS variable rvar = ~rvar exprname = '==' else: # other cases: LHS is numexpr - (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported) + (lhs, lcons) = normalized_numexpr(lexpr, csemap=csemap, supported=supported, reified=reified) return (Comparison(exprname, lhs, rvar), lcons+rcons) @@ -488,7 +486,7 @@ def normalized_boolexpr(expr, csemap=None, supported={}): return (newexpr, [c for con in flatcons for c in con]) -def normalized_numexpr(expr, csemap=None, supported={}): +def normalized_numexpr(expr, csemap=None, supported={}, reified=False): """ all 'flat normal form' numeric expressions... @@ -512,17 +510,17 @@ def normalized_numexpr(expr, csemap=None, supported={}): elif expr.is_bool(): # unusual case, but its truth-value is a valid numexpr # so reify and return the boolvar - return get_or_make_var(expr, csemap=csemap) + return get_or_make_var(expr, csemap=csemap, reified=reified) # rewrite const*a into a weighted sum, so it can be used as objective elif expr.name == "mul" and getattr(expr, "is_lhs_num", False): w, e = expr.args - return normalized_numexpr(Operator("wsum", ([w], [e])), csemap=csemap, supported=supported) + return normalized_numexpr(Operator("wsum", ([w], [e])), csemap=csemap, supported=supported, reified=reified) elif isinstance(expr, Operator): # rewrite -a into a weighted sum, so it can be used as objective if expr.name == '-': - return normalized_numexpr(Operator("wsum", _wsum_make(expr)), csemap=csemap, supported=supported) + return normalized_numexpr(Operator("wsum", _wsum_make(expr)), csemap=csemap, supported=supported, reified=reified) if not expr.has_subexpr(): return (expr, []) @@ -534,7 +532,8 @@ def normalized_numexpr(expr, csemap=None, supported={}): we = [_wsum_make(a) for a in expr.args] w = [wi for w,_ in we for wi in w] e = [ei for _,e in we for ei in e] - return normalized_numexpr(Operator("wsum", (w,e)), csemap=csemap, supported=supported) + return normalized_numexpr(Operator("wsum", (w,e)), csemap=csemap, supported=supported, reified=reified) + # wsum needs special handling because expr.args is a tuple of which only 2nd one has exprs if expr.name == 'wsum': @@ -555,14 +554,17 @@ def normalized_numexpr(expr, csemap=None, supported={}): i = i+1 # now flatten the resulting subexprs - flatvars, flatcons = map(list, zip(*[get_or_make_var(arg, csemap=csemap) for arg in sub_exprs])) # also bool, reified... + flatvars, flatcons = map(list, zip(*[get_or_make_var(arg, csemap=csemap, reified=reified) for arg in sub_exprs])) # also bool, reified... newexpr = Operator(expr.name, (weights, flatvars)) return (newexpr, [c for con in flatcons for c in con]) else: # generic operator # recursively flatten all children # TODO should be some isinstance _IntVarImpl? - flatvars, flatcons = zip(*(normalized_numexpr(arg, csemap=csemap, supported=supported) if hasattr(arg, "name") and arg.name in supported else get_or_make_var(arg, csemap=csemap) for arg in expr.args)) + flatvars, flatcons = zip(*( + normalized_numexpr(arg, csemap=csemap, supported=supported, reified=reified) + if not reified and hasattr(arg, "name") and arg.name in supported else + get_or_make_var(arg, csemap=csemap, reified=reified) for arg in expr.args)) newexpr = Operator(expr.name, flatvars) return (newexpr, [c for con in flatcons for c in con]) @@ -575,8 +577,8 @@ def normalized_numexpr(expr, csemap=None, supported={}): else: # recursively flatten all children flatvars, flatcons = zip(*[ - normalized_numexpr(arg, supported=supported) - if arg.name in supported and expr.name + normalized_numexpr(arg, supported=supported, reified=reified) + if not reified and arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap, supported=supported) for arg in expr.args]) From 480bd8061e94556dffd214f79755a7402ad88985 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 13:53:13 +0000 Subject: [PATCH 09/23] Clean up --- cpmpy/solvers/gurobi.py | 198 +++++++++++++++----------------------- tests/test_constraints.py | 2 - 2 files changed, 77 insertions(+), 123 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 5aa4874f0..7f6834bf9 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -418,137 +418,93 @@ def add(self, cpm_expr_orig): from gurobipy import GRB, nlfunc import gurobipy as gp - def add(cpm_expr, depth): - """Recursively create a Gurobi expression tree from a CPMpy expression.""" - indent = " " * depth - depth+=1 - self.grb_model.update() - print(f"{indent}Con:", cpm_expr, type(cpm_expr)) - - if is_num(cpm_expr): - return int(cpm_expr) - - def reify(grb_expr): - """Create a binary variable representing the truth value of the gurobi logical constraint (OR/AND)""" + def add(cpm_expr): + """Recursively create a Gurobi constraint from a CPMpy expression.""" + def add_(cpm_expr, depth): + indent = " " * depth + depth+=1 self.grb_model.update() - print(f"{indent}Reif", grb_expr) - - # Base case: already a variable - if isinstance(grb_expr, gp.Var): - return grb_expr - - r = self.solver_var(cp.boolvar()) - self.grb_model.addConstr(r == grb_expr) - - return r - - # Base case: variables - if isinstance(cpm_expr, NegBoolView): - return 1 - self.solver_var(cpm_expr._bv) - if isinstance(cpm_expr, _NumVarImpl): - return self.solver_var(cpm_expr) - - if isinstance(cpm_expr, Operator): - match cpm_expr.name: - case "or": # TODO only usefull if we can represent an or as a sum of binary vars, which is not generally correct - return gp.or_(add(arg, depth) for arg in cpm_expr.args) - case "and": - # self.grb_model.update() - # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) - - return gp.and_(add(arg, depth) for arg in cpm_expr.args) - args = [add(arg, depth=depth) for arg in cpm_expr.args if not is_num(arg) or arg != 1] - assert len(args) - if not args: - return 1 - elif any(a == 0 for a in args if is_num(a)): - return 0 - else: - return gp.and_(add(arg) for arg in args) - # return math.prod(add(arg, depth=depth) for arg in cpm_expr.args) - case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) - a, b = cpm_expr.args - assert isinstance(a, _BoolVarImpl), f"Implication constraint {cpm_expr} must have BoolVar as lhs, but had {a}" - assert isinstance(b, Comparison), f"Implication constraint {cpm_expr} must have Comparison as rhs, but had {b}" - is_pos = not isinstance(a, NegBoolView) - return (add(a if is_pos else a._bv, depth) == int(is_pos)) >> add(b, depth) - case "not": - # Negation: not(a) is 1 - a - return 1 - add(cpm_expr.args[0], depth) - case "-": - return -add(cpm_expr.args[0], depth=depth) - case "sum": - return sum(add(arg, depth) for arg in cpm_expr.args) - case "wsum": - return sum(weight * add(arg, depth) for weight, arg in zip(cpm_expr.args[0], cpm_expr.args[1])) - case "sub": - a, b = cpm_expr.args - return add(a, depth) - add(b, depth) - elif isinstance(cpm_expr, Comparison): - a, b = add(cpm_expr.args[0], depth), add(cpm_expr.args[1], depth) - - match cpm_expr.name: - case "==": - # Note: Gurobi functions are called by `y == f(x)`, like CPMpy boolexprs, but numexprs are transformed to `f(x) == y` (e.g. `abs(x) == IV0`), so need to be flipped - return a == b if cpm_expr.args[0].is_bool() else b == a - case "<=": - return a <= b - case ">=": - return a >= b - case _: - raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") - elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): - match cpm_expr.name: - case "mul": - return add(cpm_expr.args[0], depth) * add(cpm_expr.args[1], depth) - case "pow": - return add(cpm_expr.args[0], depth) ** add(cpm_expr.args[1], depth) - case _: # other global functions cannot be - # assert isinstance(cpm_expr, Comparison) and cpm_expr.name == "==", f"Expected global function to be transformed to f(x) == y, but was {cpm_expr}" - flatargs = [add(arg, depth) for arg in cpm_expr.args] - match cpm_expr.name: - case "abs": # y = abs(x) - return gp.abs_(add(cpm_expr.args[0], depth)) - case "min" | "max": - y = add(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) - self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) - return y + print(f"{indent}Con:", cpm_expr, type(cpm_expr)) + + if is_num(cpm_expr): + return int(cpm_expr) + elif isinstance(cpm_expr, NegBoolView): + return 1 - self.solver_var(cpm_expr._bv) + elif isinstance(cpm_expr, _NumVarImpl): + return self.solver_var(cpm_expr) + elif isinstance(cpm_expr, Operator): + match cpm_expr.name: + case "or": # TODO only usefull if we can represent an or as a sum of binary vars, which is not generally correct + return gp.or_(add_(arg, depth) for arg in cpm_expr.args) + case "and": + return gp.and_(add_(arg, depth) for arg in cpm_expr.args) + case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) + a, b = cpm_expr.args + assert isinstance(a, _BoolVarImpl), f"Implication constraint {cpm_expr} must have BoolVar as lhs, but had {a}" + assert isinstance(b, Comparison), f"Implication constraint {cpm_expr} must have Comparison as rhs, but had {b}" + is_pos = not isinstance(a, NegBoolView) + return (add_(a if is_pos else a._bv, depth) == int(is_pos)) >> add_(b, depth) + case "not": + return 1 - add_(cpm_expr.args[0], depth) + case "-": + return -add_(cpm_expr.args[0], depth=depth) + case "sum": + return sum(add_(arg, depth) for arg in cpm_expr.args) + case "wsum": + return sum(weight * add_(arg, depth) for weight, arg in zip(cpm_expr.args[0], cpm_expr.args[1])) + case "sub": + return add_(cpm_expr.args[0], depth) - add_(cpm_expr.args[1], depth) + elif isinstance(cpm_expr, Comparison): + a, b = add_(cpm_expr.args[0], depth), add_(cpm_expr.args[1], depth) + match cpm_expr.name: + case "==": + # Note: Gurobi functions are called by `y == f(x)`, like normalized CPMpy boolexprs, but CPMpy numexprs are normalized to `f(x) == y` (e.g. `abs(x) == IV0`), so they need to be flipped + return a == b if cpm_expr.args[0].is_bool() else b == a + case "<=": + return a <= b + case ">=": + return a >= b + case _: + raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") + elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): + match cpm_expr.name: + case "mul": + return add_(cpm_expr.args[0], depth) * add_(cpm_expr.args[1], depth) + case "pow": + return add_(cpm_expr.args[0], depth) ** add_(cpm_expr.args[1], depth) + case _: # other global functions cannot be + match cpm_expr.name: + case "abs": # y = abs(x) + # TODO we could support this inside the expression tree with sqrt(pow(x,2))? + return gp.abs_(add_(cpm_expr.args[0], depth)) + case "min" | "max": + y = add_(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) + flatargs = [add_(arg, depth) for arg in cpm_expr.args] + self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) + return y + else: + raise NotImplementedError(f"add_() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") + + grb_expr = add_(cpm_expr, 0) + if isinstance(grb_expr, (gp.Var, gp.LinExpr)): + # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 + return grb_expr >= 1 + elif isinstance(grb_expr, gp.TempConstr): + return grb_expr else: - raise NotImplementedError(f"add() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") + r = self.grb_model.addVar(lb=1, ub=GRB.INFINITY) + return r == grb_expr # add new user vars to the set get_variables(cpm_expr_orig, collect=self.user_vars) # transform and post the constraints for cpm_expr in self.transform(cpm_expr_orig): - grb_expr = add(cpm_expr, 0) - print("out", grb_expr) + con = add(cpm_expr) self.grb_model.update() - # If add() returned a Gurobi Var (not a constraint), wrap it as >= 1 - if isinstance(grb_expr, (gp.Var, gp.LinExpr)): - self.grb_model.addConstr(grb_expr >= 1) - elif isinstance(grb_expr, gp.TempConstr): - self.grb_model.addConstr(grb_expr) - # elif type(grb_expr).__name__ == 'GenExprAnd': - # # GenExprAnd needs addGenConstrAnd: r = AND(args) then r >= 1 - # r = self.solver_var(cp.boolvar()) - # self.grb_model.addGenConstrAnd(r, list(grb_expr)) - # self.grb_model.addConstr(r >= 1) - # elif type(grb_expr).__name__ == 'GenExprOr': - # # GenExprOr needs addGenConstrOr: r = OR(args) then r >= 1 - # r = self.solver_var(cp.boolvar()) - # self.grb_model.addGenConstrOr(r, list(grb_expr)) - # self.grb_model.addConstr(r >= 1) - else: - # r = self.solver_var(cp.boolvar()) - # self.grb_model.addConstr(r == grb_expr) - # self.grb_model.addConstr(r >= 1) - - r = self.grb_model.addVar(lb=1, ub=GRB.INFINITY) - self.grb_model.addConstr(r == grb_expr) - # self.grb_model.addGenConstrNL(r, grb_expr) - # self.grb_model.addConstr(1 == grb_expr) + print("out", con) + self.grb_model.addConstr(con) continue # Comparisons: only numeric ones as 'only_implies()' has removed the '==' reification for Boolean expressions diff --git a/tests/test_constraints.py b/tests/test_constraints.py index c922e50a2..9050ec4dd 100644 --- a/tests/test_constraints.py +++ b/tests/test_constraints.py @@ -39,8 +39,6 @@ EXCLUDE_GLOBAL = { "pysdd": NUM_GLOBAL | {"Xor"}, "minizinc": {"IncreasingStrict"}, # bug #813 reported on libminizinc - "gurobi": {"Circuit", "CumulativeOptional"} - } # Exclude certain operators for solvers. From f21f3cbf58949f6f0e1d960c545db5cdd5ef5227 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 15:11:39 +0000 Subject: [PATCH 10/23] Fixes --- cpmpy/transformations/flatten_model.py | 2 +- cpmpy/transformations/linearize.py | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 336f18467..9b965d80f 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -578,7 +578,7 @@ def normalized_numexpr(expr, csemap=None, supported={}, reified=False): # recursively flatten all children flatvars, flatcons = zip(*[ normalized_numexpr(arg, supported=supported, reified=reified) - if not reified and arg.name in supported + if not reified and hasattr(arg, "name") and arg.name in supported else get_or_make_var_or_list(arg, csemap=csemap, supported=supported) for arg in expr.args]) diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py index 80ff8ca24..ea7334775 100644 --- a/cpmpy/transformations/linearize.py +++ b/cpmpy/transformations/linearize.py @@ -332,14 +332,15 @@ def only_positive_bv(lst_of_expr, csemap=None): new_cons += new_cons_ if isinstance(rhs, _BoolVarImpl): - (new_rhs,), new_cons_ = only_positive_bv_args([rhs]) + new_args, new_cons_ = only_positive_bv_args([rhs]) + new_rhs = copy.copy(rhs) + new_rhs.update_args(new_args) new_cons += new_cons_ elif is_boolexpr(rhs): new_args, new_cons_ = only_positive_bv_args(rhs.args) new_rhs = copy.copy(rhs) new_rhs.update_args(new_args) new_cons += new_cons_ - print("NEW", new_rhs) else: new_rhs = rhs From 3055e2bd5127ffa5eaa59ad93246460c497350fc Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 1 Apr 2026 15:46:57 +0000 Subject: [PATCH 11/23] Add tests --- tests/test_gurobi.py | 222 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 tests/test_gurobi.py diff --git a/tests/test_gurobi.py b/tests/test_gurobi.py new file mode 100644 index 000000000..2ebb92749 --- /dev/null +++ b/tests/test_gurobi.py @@ -0,0 +1,222 @@ +""" +Tests for Gurobi solver transformations and expression tree support. + +These tests verify that CPMpy correctly transforms constraints to Gurobi's +expression tree format by comparing the generated LP file output. + +https://docs.gurobi.com/projects/optimizer/en/current/features/nonlinear.html +""" + +import pytest +import tempfile +import cpmpy as cp +from cpmpy.solvers.gurobi import CPM_gurobi +from cpmpy.expressions.variables import _IntVarImpl, _BoolVarImpl + + +def get_lp_string(solver): + """Write the Gurobi model to LP format and return as string.""" + solver.grb_model.update() + with tempfile.NamedTemporaryFile(suffix=".lp", delete=False, mode="w") as f: + solver.grb_model.write(f.name) + with open(f.name) as rf: + return rf.read() + + +def extract_constraints(lp_string): + """Extract constraints from 'Subject To' and 'General Constraints' sections as a list.""" + lines = lp_string.split("\n") + in_section = False + constraints = [] + for line in lines: + if line.strip() in ("Subject To", "General Constraints"): + in_section = True + continue + if line.strip() in ("Bounds", "Binaries", "Generals", "End"): + in_section = False + if in_section and line.strip(): + constraints.append(line.strip()) + return constraints + + +def reset_counters(): + _IntVarImpl.counter = 0 + _BoolVarImpl.counter = 0 + + +def expression_tree_cases(): + for c in expression_tree_cases_(): + reset_counters() + yield c + + +def expression_tree_cases_(): + """Generator yielding (name, constraint_func, expected_lp) tuples.""" + + x, y, z = [cp.intvar(-2, 2, name=name) for name in "xyz"] + p, q, r = [cp.boolvar(name=name) for name in "pqr"] + + yield ( + "pow", + x**2 + y == 9, + ["(pow(x,2)) + (y) == 9"], + ["qc0: y + [ x ^2 ] = 9"], + ) + + """Positive implications""" + yield ( + "positive_implication", + p.implies(x + y <= 3), + ["(p) -> ((x) + (y) <= 3)"], + ["GC0: p = 1 -> x + y <= 3"], + ) + + """Negative implications""" + yield ( + "negative_implication", + (~p).implies(x + y <= 3), + ["(~p) -> ((x) + (y) <= 3)"], + ["GC0: p = 0 -> x + y <= 3"], + ) + + """While NL constraint pow can be a expression tree node, it cannot be reified""" + yield ( + "imp_quad", + p.implies(x * y == 3), + ["((x) * (y)) == (IV0)", "(p) -> (sum([IV0]) == 3)"], + ["qc0: IV0 + [ - x * y ] = 0", "GC0: p = 1 -> IV0 = 3"], + ) + + yield ( + "pow_bool", + p**2 + q == 2, + ["(pow(p,2)) + (q) == 2"], + ["qc0: q + [ p ^2 ] = 2"], + ) + + yield ( + "multiplication", + z + x * y == 12, + ["(z) + ((x) * (y)) == 12"], + ["qc0: z + [ x * y ] = 12"], + ) + + yield ( + "nested", + z + x * ((-y) ** 2) == 12, + ["(z) + ((x) * (pow(sum([-1] * [y]),2))) == 12"], + [ + "\\ C3 = z + (sqr(y) * x)", + "GC0: C3 = NL : ( PLUS , -1 , -1 ) ( VARIABLE , z , 0 )", + "( MULTIPLY , -1 , 0 ) ( SQUARE , -1 , 2 ) ( VARIABLE , y , 3 )", + "( VARIABLE , x , 2 )", + ], + ) + + # # TODO needlessly reifying + # yield ( + # "subtract", + # -(x * y) == 12, + # ["(z) + ((x) * (y)) == 12"], + # ["qc0: z + [ x * y ] = 12"], + # ) + + # TODO divide (semantic may be slightly different from gurobi?) + + yield ( + "abs", + cp.Abs(x) + y == 3, + ["(IV0) + (y) == 3", "(abs(x)) == (IV0)"], + ["R0: IV0 + y = 3", "GC0: IV0 = ABS ( x )"], + ) + + """Mul is supported in expression tree, but not abs""" + yield ( + "abs_in_mul", + cp.Abs(x) * y + z == 3, + ["((IV0) * (y)) + (z) == 3", "(abs(x)) == (IV0)"], + ["qc0: z + [ IV0 * y ] = 3", "GC0: IV0 = ABS ( x )"], + ) + + yield ( + "mul_in_abs", + cp.Abs(x * y) + z == 3, + ["(IV1) + (z) == 3", "(abs(IV0)) == (IV1)", "((x) * (y)) == (IV0)"], + ["R0: IV1 + z = 3", "qc0: IV0 + [ - x * y ] = 0", "GC0: IV1 = ABS ( IV0 )"], + ) + + # TODO keep as operator? + yield ( + "minus_in", + z * (x - y) == 1, + ["(z) * (sum([1, -1] * [x, y])) == 1"], + ["qc0: [ z * x - z * y ] = 1"], # TODO not sure how it did this, but happy with it + ) + + yield ( + "minus_out", + z * -(x + y) == 1, + ["(z) * (sum([-1, -1] * [x, y])) == 1"], + ["qc0: [ - z * x - z * y ] = 1"], + ) + + yield ( + "reification", + z * (x == 2) == 1, + [ + "(z) * (BV0) == 1", + "(BV0) -> (sum([x]) == 2)", + "(~BV0) -> (sum([1, -1] * [x, BV1]) <= 1)", + "(~BV0) -> (sum([1, -5] * [x, BV1]) >= -2)", + "(BV0) -> (sum([-1] * [BV1]) >= 0)", # TODO ? + ], + [ + "qc0: [ z * BV0 ] = 1", + "GC0: BV0 = 1 -> x = 2", + "GC1: BV0 = 0 -> x - BV1 <= 1", + "GC2: BV0 = 0 -> x - 5 BV1 >= -2", + "GC3: BV0 = 1 -> - BV1 >= 0", + ], + ) + + yield ( + "disjunction", + p | q, + ["(p) or (q)"], + ["GC0: C2 = OR ( p , q )"], + ) + + # yield ( + # "conjunction", + # p & q, + # ["(p) and (q)"], + # ["R0: BV0 >= 1", "GC0: BV0 = OR ( p , q )"], + # ) + + yield ( + "conjunction_in_disjunction", + (p | (q & r)), + ["(p) or (BV0)", "(BV0) == ((q) and (r))"], + ["GC0: C2 = OR ( p , BV0 )", "GC1: BV0 = AND ( q , r )"], + ) + +@pytest.mark.requires_solver("gurobi") +@pytest.mark.parametrize( + "name,constraint,expected_tf,expected_lp", list(expression_tree_cases()), ids=[c[0] for c in expression_tree_cases()] +) +def test_gurobi_expression_tree(name, constraint, expected_tf, expected_lp): + """Test that Gurobi transformation generates expected LP output.""" + reset_counters() + solver = CPM_gurobi() + transformed = [str(c) for c in CPM_gurobi().transform(constraint)] + print("TF", ", ".join(transformed)) + reset_counters() + + solver += constraint + + lp = get_lp_string(solver) + print(lp) + + constraints = extract_constraints(lp) + assert transformed == expected_tf, f"Generated transformation:\n{transformed}" + assert constraints == expected_lp, f"Generated constraints:\n{constraints}\n\nFull LP:\n{lp}" From a056636b0b2cafb313a244f35892f75323146339 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 7 Apr 2026 13:49:23 +0000 Subject: [PATCH 12/23] Fix all tests for gurobi --- cpmpy/solvers/gurobi.py | 13 +++++++++---- cpmpy/transformations/flatten_model.py | 4 ++-- tests/test_gurobi.py | 10 ++++++---- 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 7f6834bf9..2e6c89606 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -459,8 +459,13 @@ def add_(cpm_expr, depth): a, b = add_(cpm_expr.args[0], depth), add_(cpm_expr.args[1], depth) match cpm_expr.name: case "==": - # Note: Gurobi functions are called by `y == f(x)`, like normalized CPMpy boolexprs, but CPMpy numexprs are normalized to `f(x) == y` (e.g. `abs(x) == IV0`), so they need to be flipped - return a == b if cpm_expr.args[0].is_bool() else b == a + if isinstance(a, gp.NLExpr): + # if flattening led to a non-linear expression, then it has to be constraint `y == f(x)` with `y` a `Var` + return self.grb_model.addVar(lb=b, ub=b) == a + else: + # Else, this is a function constraint + # Note: Gurobi functions are called by `y == f(x)`, like normalized CPMpy boolexprs, but CPMpy numexprs are normalized to `f(x) == y` (e.g. `abs(x) == IV0`), so they need to be flipped + return a == b if cpm_expr.args[0].is_bool() else b == a case "<=": return a <= b case ">=": @@ -479,6 +484,7 @@ def add_(cpm_expr, depth): # TODO we could support this inside the expression tree with sqrt(pow(x,2))? return gp.abs_(add_(cpm_expr.args[0], depth)) case "min" | "max": + # TODO outdated ; should be no need to make our own var y = add_(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) flatargs = [add_(arg, depth) for arg in cpm_expr.args] self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) @@ -493,8 +499,7 @@ def add_(cpm_expr, depth): elif isinstance(grb_expr, gp.TempConstr): return grb_expr else: - r = self.grb_model.addVar(lb=1, ub=GRB.INFINITY) - return r == grb_expr + return self.grb_model.addVar(lb=1, ub=1) == grb_expr # add new user vars to the set get_variables(cpm_expr_orig, collect=self.user_vars) diff --git a/cpmpy/transformations/flatten_model.py b/cpmpy/transformations/flatten_model.py index 9b965d80f..94229155b 100644 --- a/cpmpy/transformations/flatten_model.py +++ b/cpmpy/transformations/flatten_model.py @@ -289,7 +289,7 @@ def flatten_constraint(expr, csemap=None, supported={}, reified=False): return newlist -def flatten_objective(expr, supported=frozenset(["sum", "wsum"]), csemap=None): +def flatten_objective(expr, supported=frozenset(["sum", "wsum"]), csemap=None, reified=False): """ - Decision variable: Var - Linear: @@ -560,7 +560,7 @@ def normalized_numexpr(expr, csemap=None, supported={}, reified=False): else: # generic operator # recursively flatten all children - # TODO should be some isinstance _IntVarImpl? + # TODO make this some isinstance _IntVarImpl i/o hasattr? flatvars, flatcons = zip(*( normalized_numexpr(arg, csemap=csemap, supported=supported, reified=reified) if not reified and hasattr(arg, "name") and arg.name in supported else diff --git a/tests/test_gurobi.py b/tests/test_gurobi.py index 2ebb92749..c3a2a4bce 100644 --- a/tests/test_gurobi.py +++ b/tests/test_gurobi.py @@ -103,13 +103,14 @@ def expression_tree_cases_(): yield ( "nested", - z + x * ((-y) ** 2) == 12, - ["(z) + ((x) * (pow(sum([-1] * [y]),2))) == 12"], + z + (x - 3) * ((-y) ** 2) - 3 == 12, + ["(z) + (((x) + -3) * (pow(sum([-1] * [y]),2))) == 15"], [ - "\\ C3 = z + (sqr(y) * x)", + "\\ C3 = z + (sqr(y) * (-3 + x))", "GC0: C3 = NL : ( PLUS , -1 , -1 ) ( VARIABLE , z , 0 )", + # TODO not totally clean MULTIPLY node? "( MULTIPLY , -1 , 0 ) ( SQUARE , -1 , 2 ) ( VARIABLE , y , 3 )", - "( VARIABLE , x , 2 )", + "( PLUS , -1 , 2 ) ( CONSTANT , -3 , 5 ) ( VARIABLE , x , 5 )", ], ) @@ -200,6 +201,7 @@ def expression_tree_cases_(): ["GC0: C2 = OR ( p , BV0 )", "GC1: BV0 = AND ( q , r )"], ) + @pytest.mark.requires_solver("gurobi") @pytest.mark.parametrize( "name,constraint,expected_tf,expected_lp", list(expression_tree_cases()), ids=[c[0] for c in expression_tree_cases()] From 682c9da0d717ef87d2b87d59f2375756cf867e5b Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Tue, 7 Apr 2026 15:14:45 +0000 Subject: [PATCH 13/23] Update min/max --- cpmpy/solvers/gurobi.py | 19 ++++++++++++------- tests/test_gurobi.py | 7 +++++++ 2 files changed, 19 insertions(+), 7 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 88310da12..f57e23540 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -416,7 +416,6 @@ def add(self, cpm_expr_orig): :return: self """ - from gurobipy import GRB, nlfunc import gurobipy as gp def add(cpm_expr): @@ -456,6 +455,11 @@ def add_(cpm_expr, depth): return sum(weight * add_(arg, depth) for weight, arg in zip(cpm_expr.args[0], cpm_expr.args[1])) case "sub": return add_(cpm_expr.args[0], depth) - add_(cpm_expr.args[1], depth) + case "div": + assert False, "TODO" + # TODO + # if not is_num(lhs.args[1]): + # raise NotSupportedError(f"Gurobi only supports division by constants, but got {lhs.args[1]}") elif isinstance(cpm_expr, Comparison): a, b = add_(cpm_expr.args[0], depth), add_(cpm_expr.args[1], depth) match cpm_expr.name: @@ -484,12 +488,13 @@ def add_(cpm_expr, depth): case "abs": # y = abs(x) # TODO we could support this inside the expression tree with sqrt(pow(x,2))? return gp.abs_(add_(cpm_expr.args[0], depth)) - case "min" | "max": - # TODO outdated ; should be no need to make our own var - y = add_(cp.intvar(lb=min(x.lb for x in cpm_expr.args), ub=max(x.ub for x in cpm_expr.args)), depth) - flatargs = [add_(arg, depth) for arg in cpm_expr.args] - self.native_model.addConstr(y == (gp.min_(flatargs) if cpm_expr.name == "min" else gp.max_(flatargs))) - return y + case _: + args = (add_(arg, depth) for arg in cpm_expr.args) + match cpm_expr.name: + case "min": + return gp.min_(args) + case "max": + return gp.max_(args) else: raise NotImplementedError(f"add_() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") diff --git a/tests/test_gurobi.py b/tests/test_gurobi.py index 67ce17ef5..851f3bbe9 100644 --- a/tests/test_gurobi.py +++ b/tests/test_gurobi.py @@ -101,6 +101,13 @@ def expression_tree_cases_(): ["qc0: z + [ x * y ] = 12"], ) + yield ( + "maximum", + z + cp.Maximum([x, y]) == 12, + ["(z) + (IV0) == 12", "(max(x,y)) == (IV0)"], + ["R0: z + IV0 = 12", "GC0: IV0 = MAX ( x , y )"], + ) + yield ( "nested", z + (x - 3) * ((-y) ** 2) - 3 == 12, From 6f8304757c9f3ecb1b46a263d6aa9809f7080fe6 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:13:53 +0000 Subject: [PATCH 14/23] Avoid compelexity in returning reified BVs, fixing int2bool solvers --- cpmpy/solvers/gurobi.py | 5 +++-- cpmpy/transformations/linearize.py | 4 ++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index f57e23540..3f316fef9 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -442,9 +442,10 @@ def add_(cpm_expr, depth): case "->": # Gurobi indicator constraint: (Var == 0|1) >> (LinExpr sense LinExpr) a, b = cpm_expr.args assert isinstance(a, _BoolVarImpl), f"Implication constraint {cpm_expr} must have BoolVar as lhs, but had {a}" - assert isinstance(b, Comparison), f"Implication constraint {cpm_expr} must have Comparison as rhs, but had {b}" + # To not complicate linearize, we could see a _BoolVarImpl as consequent, which we rewrite here to a unary LinExpr is_pos = not isinstance(a, NegBoolView) - return (add_(a if is_pos else a._bv, depth) == int(is_pos)) >> add_(b, depth) + consequent = add_(b, depth) >= 1 if isinstance(b, _BoolVarImpl) else add_(b, depth) + return (add_(a if is_pos else a._bv, depth) == int(is_pos)) >> consequent case "not": return 1 - add_(cpm_expr.args[0], depth) case "-": diff --git a/cpmpy/transformations/linearize.py b/cpmpy/transformations/linearize.py index cb89e4332..5424c1357 100644 --- a/cpmpy/transformations/linearize.py +++ b/cpmpy/transformations/linearize.py @@ -101,7 +101,7 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal # Boolean literals are handled as trivial linears or unit clauses depending on `supported` if isinstance(cpm_expr, _BoolVarImpl): # TODO gurobi specifically cannot do reified or's - if "or" in supported and not reified: + if "or" in supported: # post clause explicitly (don't use cp.any, which will just return the BoolVar) newlist.append(cpm_expr) elif isinstance(cpm_expr, NegBoolView): @@ -131,7 +131,7 @@ def linearize_constraint(lst_of_expr, supported={"sum","wsum","->"}, reified=Fal if isinstance(cond, _BoolVarImpl) and isinstance(sub_expr, _BoolVarImpl): # shortcut for BV -> BV, convert to disjunction and linearize it (if unsupported) if "->" in supported: - newlist.append(cond.implies(sub_expr >= 1)) + newlist.append(cond.implies(cp.all(linearize_constraint([sub_expr], supported=supported, reified=True, csemap=csemap)))) else: newlist.append(1 * cond + -1 * sub_expr <= 0) From b50668b6d10f19edec15051b4d767b604ac51f14 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:17:13 +0000 Subject: [PATCH 15/23] Fix merge (faulty boolexpr reification) --- cpmpy/transformations/reification.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/cpmpy/transformations/reification.py b/cpmpy/transformations/reification.py index 9b6c67f74..efb309f39 100644 --- a/cpmpy/transformations/reification.py +++ b/cpmpy/transformations/reification.py @@ -176,14 +176,18 @@ def reify_rewrite(constraints, supported=frozenset(), csemap=None): boolexpr = cpm_expr.args[boolexpr_index] if isinstance(boolexpr, Operator): if boolexpr.name in supported or cpm_expr.name == "==": - # Case 1, BE is Operator (and, or, ->) + # Case 1a, BE is Operator (and, or, ->) newcons.append(cpm_expr) # could actually rewrite into list of clauses like to_cnf() does... not for here elif cpm_expr.name == "->": + # Case 1b, BE is an unflattened expression (TODO duplicated from below) + # We have BV -> BE, create BV -> auxvar, auxvar == BE (auxvar, cons) = get_or_make_var(boolexpr, csemap=csemap) newcons += cons reifexpr = copy.copy(cpm_expr) - reifexpr.args[boolexpr_index] = auxvar + args = list(reifexpr.args) + args[boolexpr_index] = auxvar + reifexpr.update_args(tuple(args)) newcons.append(reifexpr) elif isinstance(boolexpr, GlobalConstraint): # Case 2, BE is a GlobalConstraint From f050a9333bc7c1a52ffd6855b7aa7729a11d4cc7 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:52:19 +0000 Subject: [PATCH 16/23] Revert global exclusion --- cpmpy/solvers/gurobi.py | 4 +++- tests/test_constraints.py | 3 +-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 3f316fef9..74775d3ad 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -467,7 +467,9 @@ def add_(cpm_expr, depth): case "==": if isinstance(a, gp.NLExpr): # if flattening led to a non-linear expression, then it has to be constraint `y == f(x)` with `y` a `Var` - return self.grb_model.addVar(lb=b, ub=b) == a + y = b if isinstance(b, gp.Var) else self.grb_model.addVar(lb=b, ub=b) + # TODO unclear why this is no longer normalized to be a constant `b` + return y == a else: # Else, this is a function constraint # Note: Gurobi functions are called by `y == f(x)`, like normalized CPMpy boolexprs, but CPMpy numexprs are normalized to `f(x) == y` (e.g. `abs(x) == IV0`), so they need to be flipped diff --git a/tests/test_constraints.py b/tests/test_constraints.py index 0493e78d3..2089f1e9e 100644 --- a/tests/test_constraints.py +++ b/tests/test_constraints.py @@ -23,8 +23,7 @@ NUM_GLOBAL = { "AllEqual", "AllDifferent", "AllDifferentExcept0", "AllDifferentExceptN", "AllEqualExceptN", - "GlobalCardinalityCount", "InDomain", "Inverse", - # "Circuit", + "GlobalCardinalityCount", "InDomain", "Inverse", "Circuit", "Table", 'NegativeTable', "ShortTable", "Regular", "Increasing", "IncreasingStrict", "Decreasing", "DecreasingStrict", "Precedence", "Cumulative", "NoOverlap", "CumulativeOptional", "NoOverlapOptional", From 2df51ebb17e7993794bd59238203665a9c67835b Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:56:22 +0000 Subject: [PATCH 17/23] Add more tests --- tests/test_gurobi.py | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/tests/test_gurobi.py b/tests/test_gurobi.py index 851f3bbe9..170749819 100644 --- a/tests/test_gurobi.py +++ b/tests/test_gurobi.py @@ -56,6 +56,27 @@ def expression_tree_cases_(): x, y, z = [cp.intvar(-2, 2, name=name) for name in "xyz"] p, q, r = [cp.boolvar(name=name) for name in "pqr"] + yield ( + "BV", + p, + ["p"], + ["R0: p >= 1"], + ) + + yield ( + "True", + cp.BoolVal(True), + ["boolval(True)"], + ["R0: C0 = 1"], + ) + + yield ( + "False", + cp.BoolVal(False), + ["boolval(False)"], + ["R0: C0 = 0"], + ) + yield ( "pow", x**2 + y == 9, From 2236b6c76d131b592aad77d6260f1da0f5e4a9e3 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:56:34 +0000 Subject: [PATCH 18/23] Support DirectConstraint --- cpmpy/solvers/gurobi.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 74775d3ad..d0b73236c 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -498,6 +498,10 @@ def add_(cpm_expr, depth): return gp.min_(args) case "max": return gp.max_(args) + + elif isinstance(cpm_expr, DirectConstraint): + cpm_expr.callSolver(self, self.grb_model) + return True else: raise NotImplementedError(f"add_() not implemented for {cpm_expr}, {type(cpm_expr)}, {getattr(cpm_expr, 'name', None)}") From abe0d084f2425a1591c79a6056e8bdab6504a933 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 10:59:25 +0000 Subject: [PATCH 19/23] Refactor --- cpmpy/solvers/gurobi.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index d0b73236c..681959518 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -481,23 +481,20 @@ def add_(cpm_expr, depth): case _: raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): + args = [add_(a, depht) for a in cpm_expr.args] match cpm_expr.name: case "mul": - return add_(cpm_expr.args[0], depth) * add_(cpm_expr.args[1], depth) + return args[0] * args[1] case "pow": - return add_(cpm_expr.args[0], depth) ** add_(cpm_expr.args[1], depth) - case _: # other global functions cannot be - match cpm_expr.name: - case "abs": # y = abs(x) - # TODO we could support this inside the expression tree with sqrt(pow(x,2))? - return gp.abs_(add_(cpm_expr.args[0], depth)) - case _: - args = (add_(arg, depth) for arg in cpm_expr.args) - match cpm_expr.name: - case "min": - return gp.min_(args) - case "max": - return gp.max_(args) + return args[0] ** args[1] + # remaining global function cannot be part of the expression tree + case "abs": # y = abs(x) + # TODO we could support this inside the expression tree with sqrt(pow(x,2))? + return gp.abs_(args[0]) + case "min": + return gp.min_(args) + case "max": + return gp.max_(args) elif isinstance(cpm_expr, DirectConstraint): cpm_expr.callSolver(self, self.grb_model) From 230cd8b8fce79efacc4081b75ff67b326746e05a Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 11:03:12 +0000 Subject: [PATCH 20/23] Typo --- cpmpy/solvers/gurobi.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 681959518..1aabc770e 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -481,7 +481,7 @@ def add_(cpm_expr, depth): case _: raise Exception(f"Expected comparator to be ==,<=,>= in Comparison expression {cpm_expr}, but was {cpm_expr.name}") elif isinstance(cpm_expr, cp.expressions.globalfunctions.GlobalFunction): - args = [add_(a, depht) for a in cpm_expr.args] + args = [add_(a, depth) for a in cpm_expr.args] match cpm_expr.name: case "mul": return args[0] * args[1] From 01f762819d15c76701a496fc8e818ea79f74b9c0 Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 11:09:14 +0000 Subject: [PATCH 21/23] Fix trans tests --- tests/test_trans_linearize.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tests/test_trans_linearize.py b/tests/test_trans_linearize.py index 358e1f332..0d3144cf7 100644 --- a/tests/test_trans_linearize.py +++ b/tests/test_trans_linearize.py @@ -34,7 +34,7 @@ def test_linearize(self): # implies cons = linearize_constraint([a.implies(b)])[0] - assert "sum([1, -1] * [a, b]) <= 0" == str(cons) + assert "(a) -> (b >= 1)" == str(cons) def test_bug_168(self): from cpmpy.solvers import CPM_gurobi @@ -69,8 +69,8 @@ def test_constraint(self): assert str(linearize_constraint([a | b | c])) == "[sum(a, b, c) >= 1]" assert str(linearize_constraint([a | b | (~c)])) == "[sum(a, b, ~c) >= 1]" # test implies - assert str(linearize_constraint([a.implies(b)])) == "[sum([1, -1] * [a, b]) <= 0]" - assert str(linearize_constraint([a.implies(~b)])) == "[sum([1, -1] * [a, ~b]) <= 0]" + assert str(linearize_constraint([a.implies(b)])) == "[(a) -> (b >= 1)]" + assert str(linearize_constraint([a.implies(~b)])) == "[(a) -> (b <= 0)]" assert str(linearize_constraint([a.implies(x+y+z >= 0)])) == str([]) assert str(linearize_constraint([a.implies(x+y+z >= 2)])) == "[(a) -> (sum(x, y, z) >= 2)]" assert str(linearize_constraint([a.implies(x+y+z > 0)])) == "[(a) -> (sum(x, y, z) >= 1)]" @@ -87,7 +87,7 @@ def test_constraint(self): c1, c2, c3 = linearize_constraint([a.implies(x != y)]) assert str(c1) == "(a) -> (sum([1, -1, -6] * [x, y, BV4]) <= -1)" assert str(c2) == "(a) -> (sum([1, -1, -6] * [x, y, BV4]) >= -5)" - assert str(c3) == "sum([1, -1] * [~a, ~BV4]) <= 0" + assert str(c3) == "(~a) -> (BV4 <= 0)" def test_single_boolvar(self): @@ -95,8 +95,8 @@ def test_single_boolvar(self): p = cp.boolvar(name="p") assert str([p >= 1]) == str(linearize_constraint([p])) assert str([p <= 0]) == str(linearize_constraint([~p])) - assert str([Operator("or", [p])]) == str(linearize_constraint([p], supported={"or"})) - assert str([Operator("or", [~p])]) == str(linearize_constraint([~p], supported={"or"})) + assert str([p]) == str(linearize_constraint([p], supported={"or"})) + assert str([~p]) == str(linearize_constraint([~p], supported={"or"})) def test_neq(self): # not equals is a tricky constraint to linearize, do some extra tests on it here From 86f145d197080018404cdc00396e53d6d55bd84a Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 12:12:20 +0000 Subject: [PATCH 22/23] Fix y/fx in == --- cpmpy/solvers/gurobi.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/cpmpy/solvers/gurobi.py b/cpmpy/solvers/gurobi.py index 1aabc770e..b03c7fb69 100644 --- a/cpmpy/solvers/gurobi.py +++ b/cpmpy/solvers/gurobi.py @@ -468,12 +468,14 @@ def add_(cpm_expr, depth): if isinstance(a, gp.NLExpr): # if flattening led to a non-linear expression, then it has to be constraint `y == f(x)` with `y` a `Var` y = b if isinstance(b, gp.Var) else self.grb_model.addVar(lb=b, ub=b) - # TODO unclear why this is no longer normalized to be a constant `b` + # TODO unclear why this is no longer normalized to be a constant `b` (same below) return y == a else: - # Else, this is a function constraint + # Else, this is a function constraint # Note: Gurobi functions are called by `y == f(x)`, like normalized CPMpy boolexprs, but CPMpy numexprs are normalized to `f(x) == y` (e.g. `abs(x) == IV0`), so they need to be flipped - return a == b if cpm_expr.args[0].is_bool() else b == a + y, fx = (a, b) if cpm_expr.args[0].is_bool() else (b, a) + y = y if isinstance(y, gp.Var) else self.grb_model.addVar(lb=y, ub=y) + return y == fx case "<=": return a <= b case ">=": From 6e6f6bb00952c7437c648db34e19b309e0f882db Mon Sep 17 00:00:00 2001 From: Hendrik 'Henk' Bierlee Date: Wed, 8 Apr 2026 12:18:19 +0000 Subject: [PATCH 23/23] Fix supported in test --- tests/test_transf_reif.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_transf_reif.py b/tests/test_transf_reif.py index 6bc206d67..ebc5a7b1e 100644 --- a/tests/test_transf_reif.py +++ b/tests/test_transf_reif.py @@ -82,7 +82,7 @@ def test_reif_rewrite(self): rv = cp.boolvar(name="rv") arr = cp.cpm_array([0,1,2]) - f = lambda expr : str(reify_rewrite(flatten_constraint(expr))) + f = lambda expr : str(reify_rewrite(flatten_constraint(expr), supported={"or"})) fd = lambda expr : str(reify_rewrite(flatten_constraint(decompose_in_tree(expr))))