ConstraintSolver.jl Reified constraint

Creation date: 2020-07-13

Tags: julia, constraint-programming, jump

I think this is part 25 of this long long series. If you haven't checked out the series before: Welcome :) If you're interested in constraint programming you will have a lot to read. Start here: Solving sudokus.

There and in the sidebar you will find all the other posts in this series. I haven't worked on the solver for quite a bit and therefore the last blog post is a while back. Some changes were a bit more complicated than expected and I needed a break. Additionally, for the last two constraints, indicator and now this one, I would like to have some puzzle that can be solved with it to test the constraints more and be able to have a visual appealing blog post :D

I know that these blog posts might be a bit boring for some of you, especially if you're not fully into constraint programming. I often try to make them such that there is some value for other readers as well. This post is no different! In this post you can learn: How to extend JuMP.jl in a way it supports your own type of constraint. In the past we have seen several set constraints.

Recap indicator constraint

This way we can support constraints of the form:

@constraint(m, [x,y] in Set())

where Set can be CS.AllDifferentSet as an example or CS.TableSet which needs an argument like CS.TableSet(table).

Another thing we implemented was the support for != which is in the same section of "How to use JuMP for your own solver"

We need to support a simple looking != as this is not a standard constraint in (non-)linear solvers, which are normally used with JuMP.

This time I wanted to support a so call reified constraint which is quite similar to the indicator constraint, which I have discussed in the last two posts.

The indicator constraint looks like:

@constraint(m, a => {x + y <= 12})

which reads as if a == 1 then make sure that x+y <= 12 otherwise don't care about x and y. Here a needs to be a binary value and the inner constraint (don't know whether there is a good name for it actually...) must be a constraint and surrounded by {}.

One thing this doesn't support is: Iff a constraint is satisfied a binary variable should be active.

Iff stands for if and only if.

Extend JuMP on a different level

This can be done by the reified constraint which looks similar:

@constraint(m, a := {x + y <= 12})

So this means if a is active then x+y <= 12 and if x+y <= 12 then a is active. Active means it's set to true or to false if you write it:

@constraint(m, !a := {x + y <= 12})

Indicator constraints are partially supported by JuMP. With partially I mean they allow the above version but not:

@constraint(m, a => {[x,y] in CS.AllDifferentSet()})

which I made available for the constraint solver and described in the last post. Reified constraints are not supported by JuMP but since v0.21.3 it is possible define them relatively easy. Here is the PR that made it possible JuMP.jl#2228.

That PR makes it possible to use your own symbol between a left and a rhs, so in our case :=.

The method declaration for this:

function JuMP.parse_constraint_head(_error::Function, ::Val{:(:=)}, lhs, rhs)

For dispatching reasons the symbol is surrounded by Val.

The inside of the function is basically the same as for the indicator constraint. The only difference is that we want to create a ReifiedSet and not an IndicatorSet.

We need two _build_reified_constraint methods for this as we want to support scalar inner constraints as well as the more general set type of constraints.

function _build_reified_constraint(
    _error::Function, variable::JuMP.AbstractVariableRef,
    constraint::JuMP.ScalarConstraint, ::Type{CS.ReifiedSet{A}}) where A
    set = ReifiedSet{A}(JuMP.jump_function(constraint), JuMP.moi_set(constraint), 2)
    return JuMP.VectorConstraint([variable, JuMP.jump_function(constraint)], set)

for vector constraints so something like a => { [x,y] in CS.AllDifferentSet()}.

function _build_reified_constraint(
    _error::Function, variable::JuMP.AbstractVariableRef,
    constraint::JuMP.VectorConstraint, ::Type{CS.ReifiedSet{A}}) where A
    set = CS.ReifiedSet{A}(MOI.VectorOfVariables(constraint.func), constraint.set, 1+length(constraint.func))
    vov = VariableRef[variable]
    append!(vov, constraint.func)
    return JuMP.VectorConstraint(vov, set)

⚠ Important note

In the indicator post I mentioned that I save the activation variable in the IndicatorSet but don't use it. I found out that I can't use it because the index is not the correct index.

It does not work because JuMP has an inner version of the model at the beginning and then copies it over to the solver once optimize! is called. This copy changes the indices which means that the saved index is not compatible. It is therefore important to have the indices only inside the return statement. This makes sure that JuMP handles the indices correctly in a later process and when our MOI.add_constraint methods are called the indices are correct and can be used as we wish.

Another problem I encountered was that everything is done inside the @constraint macro which is basically called when including the file and not when calling the function and it has it's own quirks. If you want to work with macros you should read the documentation and probably don't trust me.

I encountered that I'll need a call like

buildcall = :(_build_indicator_constraint($_error, $(esc(variable)), $rhs_buildcall, $S))

where the buildcall is returned later and will be used to build the constraint once the function is actually called. As this is inside a JuMP function which I extended the call to _build_reified_constraint would try to call a JuMP function with that name which doesn't exist and I can't create a new function of a different package (without using eval). I can only create a new method of an existing function. That meant I was a little stuck and asked on one of the help places I mentioned in my basics posts.

If you encounter this issue as well you might want to just check the discourse thread.

I learned a bit about macro hygiene and the useful macro @macroexpand but I don't feel qualified to explain anything in more detail.

The result was this line:

buildcall = :($(esc(:(CS._build_reified_constraint)))($_error, $(esc(variable)), $rhs_buildcall, $S))

So I head to escape my function name such that there is no error that CS does not exist inside JuMP but don't escape the other variables and calls as they do exist in JuMP.

Actual implementation of the constraint

Feel free to jump to the code using See on GitHub to see the full code. In this section I want to give a brief overview of how this constraint is currently implemented in the solver. If you are only interested in how to extend JuMP this part might not be that interesting for you.

We first of all need the add_constraint method. I've explained that in previous posts and I think it's not that complicated. Feel free to reach out if you have questions ;)

function MOI.add_constraint(
) where {A, T<:Real, RS<:ReifiedSet{A}}

Afterwards I needed to implement all the constraint functions like:

function init_constraint!(
    fct::Union{MOI.VectorOfVariables, VAF{T}},
) where {A, T<:Real, RS<:ReifiedSet{A}}

if you check it out on GitHub, you'll see all the implementations in that file.

What the constraint does:

you might see that this is not the best way of pruning in some cases. There might be ways to anti-prune a constraint in the inactive case but this is not supported by the constraints and therefore not by the reified constraint. This can be done in future versions.

Then there are methods like single_reverse_pruning_constraint which reverses a single variable when jumping through the backtrack tree. This needs to be called for the inner constraint if it is affected.

For some constraint not all of those methods are needed but as the inner constraint can be any constraint the reified constraint needs to implement them and call the corresponding inner constraint method when needed.

The final step is always to create more test cases, which are sometimes also created in between and at least a few in the beginning before anything works.

I hope that I fixed all the bugs as sometimes the table constraint with all its methods for reversing is a bit tricky.

As mentioned earlier this constraint type needs a newer version of JuMP as the other constraints before. Fortunately it's not a problem because I can specify the versions in the Project.toml:

Formatting = "^0.4.1"
JSON = "~0.18, ~0.19, ~0.20, ~0.21"
JuMP = "^0.21.3"
MathOptInterface = "^0.9.11"
MatrixNetworks = "^1"
julia = "1"

When changing that section one just needs to change the minor version number such that this is now available in v0.3.0.

Next up

There are still a lot of things to work on for the constraint solver:

There is also one issue which exists since the early beginning and was suggested by my former supervisor at LANL: Having a FlatZinc reader #27.

Which seems like a very useful issue but also complicated, which is the reason why I haven't touched it yet.

I'm thinking about doing the refactoring in a live stream as this could also be a good way of explaining the general structure of the code. Maybe someone is interested in the project and would like to add something to it. ;)

Additionally, I'm very happy to say that 10 of you now support me and this blog. Thank you so much! It was unimaginable a year ago.

Thanks for reading and special thanks to my 10 patrons!

List of patrons paying more than 4$ per month:

Currently I get 28$ per Month using Patreon and PayPal.

For a donation of a single dollar per month you get early access to the posts. Try it out at the start of a month and if you don't enjoy it just cancel your subscription before pay day (end of month).

I'll keep you updated on Twitter OpenSourcES as well as my more personal one: Twitter Wikunia_de

Want to be updated? Consider subscribing and receiving a mail whenever a new post comes out.

Powered by Buttondown.

Subscribe to RSS