CPMpy / cpmpy

Constraint Programming and Modeling library in Python, based on numpy, with direct solver access.
Apache License 2.0
226 stars 23 forks source link

Incomplete functions on rhs of implication #184

Open IgnaceBleukx opened 1 year ago

IgnaceBleukx commented 1 year ago

Some expressions are incomplete, that is, there exist value for their arguments for which there is no function value defined. Two cases I can think of and that raises errors currently:

Element constraint index out of range:

a = boolvar()
arr = intvar(0,3, shape=4)
idx = intvar(0,10)

cons = a.implies(arr[idx] >= 1)
assert Model([cons, idx == 10]).solve() 
cons.value() # index out of range error

This first case can/should be handled when calculating the value of a constraint by lazily evaluating the value of the tail in the constraint.

Division by zero

a = boolvar()
x = intvar(0,10)

from cpmpy.transformations.reification import reify_rewrite
cons = a.implies(10 / x <= 5)
reify_rewrite(cons) # results in division by 0 error

The second case is another annoying bounds calculation bug.

tias commented 1 year ago

Our KR colleague Djordje is doing his entire PhD on partial functions, so you can always catch up with him about the full issue ; )

The top one, the element, the answer of the model is correct nowadays if transformations/reification.py transformations are used. That is because there is no solver supporting this, so I added a 'decompose_comparison' to the Element class and that one loops over all indices in a way that idx can take a value outside of the range of array...

So with that effort done, all that is left is indeed to special-case the value() function of operator for -> to recursively call value as needed...

For the second one... not sure what you would propose. We could also be very strict on not allowing division operators where the numerator has a 0 in its domain... at least one of the solvers does this

tias commented 1 year ago

was still open...

we already have an 'argval' in expr/utils. If we replace arg_vals[0] by argval(self.args[0]) everywhere instead of computing it all upfront, and then we might automatically get it also for and, or, etc if Python does it lazily...

IgnaceBleukx commented 1 year ago

Yes and no... There are still cases where the value computation can go wrong: cond -> subexpr with cond = [1,2,3][l] == 2 and subexpr = p: Valid solution (I think) is l = 10 and p = True. But computing the value of cond will raise an IndexError. What I think would be a good solution to this problem is to introduce a IncompleteFunctionError when calculating the value of cond and catch this when computing the value of the entire expression

JoD commented 1 year ago

I'd like to reopen this issue. Supporting partial functions is a huge pain in the butt. If Microsoft with Z3 cannot do it, we should not either.

One way this manifests is that the way Element is handled right now leads to flattening not being equivalence preserving. E.g., right now, the constraint not (x * element([1,2,3],y) >= 5) is true if y<0 or y>2. Unnesting yields the following constraints: not (x * aux >= 5) and aux == element([1,2,3],y). If y<0 or y>2, these constraints now evaluate to false, making their conjunction unsatisfiable.

I've proposed four ways to make Element total this afternoon. My favorite fix is to give Element a default value (0 or the 0-index expression) if the index is outside of the lists's domain. Another is to treat it the same as division by 0 (i.e., not allow it). A third is to let it take any numerical value, basically creating an unknown numeric variable. A fourth is to only allow Element as part of an equality, turning the constraint into a well-defined Boolean (this is the way the global constraint catalog handles it). Any of these is fine by me.

IgnaceBleukx commented 1 year ago

Aha, I did not think of negative contexts yet... This makes it indeed rather tricky.

It seems Minizinc also goes for the default approach but I'm really not a fan of this. Especially for Element this will be very confusing. For example [0,0,0][y] == 0 & y >= 10 will be satisfiable! Giving it any numerical value has the same problem as the above approach.

Following the GCC and making the element constraint boolean is also not my favorite. As in CPMpy we would like to stick closely to Python/Numpy semantics. Moreover, as we subclass np.ndarray we can do very fancy things with numpy indexing to cleanly write our models.

If we really want Element to be total, we should indeed do it as done with division and not allow a variable with bounds greater than the length of the array or negative lower-bound. Although I'm still not convinced we should make it total.

For every negative or reified context where the function is partial (which we can do at runtime by computing bounds of the index) simply decompose the element constraint into several total functions. Your example would then be:

(y == 0) -> 1 < 5
(y == 1) -> 2 < 5
(y == 2) -> 3 < 5

Which is always true, no matter the value of y. I think @Wout4 is working on streamlining negation during flattening so this fix would be in that pull request (can't find it atm).

We can print a warning to the user the element constraint is decomposed and that he can fix it by setting the bounds of the index properly.

JoD commented 1 year ago

To check whether I understand:

not not ([0,0,0][y] == 0)
y>=10

would decompose to

not not (aux == 0)
y>=10
(y==0) -> aux=0
(y==1) -> aux=0
(y==2) -> aux=0

which is satisfiable.

However, after pushing negation we would get:

[0,0,0][y] == 0
y>=10

which would not be satisfiable, as the decomposition is no longer done.

Letting the semantics of an expression depend on its context is pretty tricky.

I'm ok with the division approach.

IgnaceBleukx commented 1 year ago

OK, new insights from trying to write a recursive decompose_global transformation: its a mess to deal with the partial functions... Lets indeed make Element total by construction. I.e., make sure the index is >0 and <len(arr)

JoD commented 1 year ago

Good! In a nutshell, the argument is: "it is hard to write recursive transformations when individual nodes in the term tree can potentially not have a value/range/etc.".

@tias: is this ok by you?

IgnaceBleukx commented 1 year ago

Ok, still not sure this restrictive version is the way to go... I have two main concerns about the restrictive version:

After reading the Minizinc paper: Stuckey, Peter J., and Guido Tack. "Minizinc with functions." Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems: 10th International Conference, CPAIOR 2013, Yorktown Heights, NY, USA, May 18-22, 2013. Proceedings 10. Springer Berlin Heidelberg, 2013.

We could actually make a fairly neat transformation that fixes all issues with partial functions by introducing a guard for them at top-level.

I am fine with changing to the restrictive version as intermediate solution, so decompose_in_tree can be merged, but we should revisit it afterwards.