Table constraint

This is part 20 of an ongoing series about: How to build a constraint solver? You might want to check out earlier posts if you're new. See first post.

I didn't only code a Sudoku solver in Python two years ago. Basically I tried a lot of what I do now again but this time I'm more structured and Julia feels like a much better language for it. We need the table constraint to solve variants of Sudoku puzzles this time and I'll show you in this post how the constraint works in it's basics and how you can add your own constraint to the solver.

Normally I have a puzzle in mind which I want to solve and build the constraint around it. This time I looked which constraints might be useful and wanted to do something about scheduling problems as constraint solvers are the most used tool for those problems.

Unfortunately I wasn't able to find the best strategy to have interval variables inside JuMP. I will think about it next again. The next thing I found was the table constraint which is a really powerful constraint.

Currently we can limit the search space of a single variable easily by just setting bounds like:

@variable(m, 1 <= x[1:5] <= 7, Int)
# or
@variable(m, x[1:5], CS.Integers([1,2,3,4,5,6,7]))

with the table constraint you can define a set of possibilities for several variables combined. This can be also used to create other constraints by just enumerating all possibilities like all possible combinations which fulfill the all different constraint.

A small example constraint:

table = [
    1 2 3 1 1;
    1 3 3 2 1;
    1 1 3 2 1;
    1 1 1 2 4;
    1 8 3 1 2;
    4 5 5 3 4;
]

@constraint(model, x in CS.TableSet(table))

This concrete set of options can't be transformed into ==, <= or all different constraints in an easy way.

You'll see in the next posts how this constraint can be used for solving other puzzles.

Naive way of implementing the constraint

We always need a way to check feasibility and prune values if we got new knowledge. A very naive way would be to run over all possibilities and check whether the current search space is feasible.

Let's have a look at the table again:

table = [
    1 2 3 1 1;
    1 3 3 2 1;
    1 1 3 2 1;
    1 1 1 2 4;
    1 8 3 1 2;
    4 5 5 3 4;
]

For pruning: if setting x[1] we can remove the last row and need to check which possible values are only used in this row i.e x[2] = 5 is not used by any other row therefore we can remove it from the search space of x[2].

The goal is always to prune as much as possible and if we have millions of rows it takes a long time to do these checks and keep track of which rows are possible and which aren't. We can implement a way to not iterate over rows which aren't possible anymore and so on but in general it feels like there is probably a better way to do it.

After thinking about sparse matrices and some other stuff for not too long I decided to just read the paper I've found beforehand.

Reasonable way of how the constraint can be implemented

I've found out about the table constraint on minicp.org which I mentioned here before. They've linked to the following paper: Compact-Table: Efficiently Filtering Table Constraints with Reversible Sparse Bit-Sets.

I recommend you to read that paper if you want to know all the details but I will also explain the rough idea:

We'll do quite a few things when the constraint is added to the model to make it faster in the feasibility and pruning steps.

• Remove rows which can not be fulfilled because of the initial variable constraints

• Compute a so called support for each variable and value combination

In the above example:

@variable(m, 1 <= x[1:5] <= 7, Int)
table = [
    1 2 3 1 1;
    1 3 3 2 1;
    1 1 3 2 1;
    1 1 1 2 4;
    1 8 3 1 2;
    4 5 5 3 4;
]
@constraint(model, x in CS.TableSet(table))

We will first iterate over all rows and check whether the corresponding variable (column) has the value associated with that row. This will remove 1 8 3 1 2 as x[2] is bounded above by 7.

You might wonder why the user gives the stupid table in the first place but it's quite handy if you want to easily define your model and let the solver do the work.

I've tested my implementation on some small examples like the one above first and then used it to solve sudokus. What else? :D

There I defined the table once with:

table = permutedims(reduce(hcat, collect(permutations(1:9))))

and used it for all 27 constraints and let the solver remove rows which based on the start setting of the sudoku aren't possible anymore.

Support

Next we define the support on the remaining rows:

[
    1 2 3 1 1;
    1 3 3 2 1;
    1 1 3 2 1;
    1 1 1 2 4;
    4 5 5 3 4;
]

in the following way: We list all 5 variables and their values in one column:

var = val	support
x[1] = 1	1 1 1 1 0
x[1] = 2	0 0 0 0 0
x[1] = 3	0 0 0 0 0
x[1] = 4	0 0 0 0 1
x[1] = 5	0 0 0 0 0
x[1] = 6	0 0 0 0 0
x[1] = 7	0 0 0 0 0
x[2] = 1	0 0 1 1 0
x[2] = 2	1 0 0 0 0
x[2] = 3	0 1 0 0 0
x[2] = 4	0 0 0 0 0
x[2] = 5	0 0 0 0 1
x[2] = 6	0 0 0 0 0
x[2] = 7	0 0 0 0 0

... the same for the other three variables.

For each variable value pair the support shows in which rows of the table constraint this pair is set i.e x[2] = 5 is only true in the last row which is represented by 0 0 0 0 1 whereas x[1] = 1 is true in all but the last row => 1 1 1 1 0. Which means that x[1] = 1, x[2] = 5 is never possible.

You can see that there are quite a lot of supports which are 0 0 0 0 0 so no support which is the case because that setting is never possible. These pairs can be removed immediately.

Now those permutations are not actually stored like you might imagine (in a boolean vector) but instead in a vector of UInt64 and we use the bits. This later has some advancements in checking fast whether a whole chunk of rows are invalid.

Feasibility checking

If we want to check whether a setting 1:3, 2 is valid we can use | over the supports of x[1]={1,2,3} and & them with x[2] = 2.

In our case for x[1]:

_ 1 1 1 1 0 
| 0 0 0 0 0
| 0 0 0 0 0
= 1 1 1 1 0

for x[2]

1 0 0 0 0

and then &

_ 1 1 1 1 0
& 1 0 0 0 0
= 1 0 0 0 0
!= 0 0 0 0 0
=> feasible

This assumes that we currently haven't removed any other row of the initial table which here is a column. If we removed the first row already otherwise we would compute:

_ 0 1 1 1 1 
& 1 1 1 1 0
& 1 0 0 0 0
= 0 0 0 0 0
=> infeasible

Which brings us to some things we store inside our constraint at all times and update when something changes:

• current possible rows are represented by a UInt64 vector called words

  • it has some specialties which I explain after the pruning part

• a mask (UInt64 vector of the same length) which stores the in between results

So the above example more precisely:

mask___: 1 1 1 1 1
______ & 1 1 1 1 0
______ & 1 0 0 0 0
words_ & 0 1 1 1 1
words` & 0 0 0 0 0
=> infeasible

Pruning

In pruning we want to remove values from the search space based on new information about the search space from other constraints or from trial and error in the backtracking phase.

First we:

• check which variables have changed since last time by checking whether the size of the possible values have changed

• update the domain size for next time

• check which variables are not fixed

For the changed variables we compute the mask for all possible values to check feasibility as described before.

Next step is to filter the domains based on the new knowledge which uses all unfixed variables. For each variable value pair we check feasibility. If not feasible we remove the value from the search space.

Now this can be slow which brings us to the next point.

How to store extra information to make it fast?

We kind of assume/hope that we remove quite a lot of rows in the first few backtracking steps ideally in the root node. Let's say we have 1million rows in our table. If we can't remove any of them in init the UInt64 vector has about 16000 elements. In the first prune run it might be the case that we can remove 90% of the rows. Then it would be quite stupid to iterate over the 16000 elements every time and do some & operation if words[idx] == 0. Instead we can use the a similar trick as with our general variables where we have values and indices to store which values are still possible. Here we also have words and indices and indices simply shows which indices of words are non zero.

A small example:

We start with: (1 in words simply indicates that it is not zero.)

words__: [1,1,1,1,1,1]
indices: [1,2,3,4,5,6]

and then we remove some rows such that some words (64 rows in the same block) get removed/0.

words__: [1,1,0,1,1,0]
indices: [1,2,4,5,3,6]

Now we iterate over the indices with for idx in indices and then use words[idx]. Additionally we save how many words are non zero to stop early.

Another thing we can improve is speed up our check if a value is feasible/possible is that we check whether there exists at least one still feasible row which has the variable/value pair which we currently check.

For this a naive approach is to use words & support[var, val] and check whether it is the zero vector but if it's not the zero vector because the last possible index is non zero we would basically iterate over a lot of nonsense again every time we check that. But if we forget reverse pruning for now we will never change the words in a way such that it will be non zero before the last index if it is the last index now. That means that we can save the position of the last non zero and every time start from there and update that last position if it changes.

If you want to have all the details I suggest that you read the paper as it's not that hard to understand.

I sure made a lot of bugs in implementing it and fixes over fixes but that's just me :D

How to add it to the ConstraintSolver.jl package

If you want to add your own constraint to the package this is the correct section. I'll also add some docs in the official docs but this is a rough overview on how I think the workflow should be at the moment.

This will also give you a better overview of how the constraint solver works internally if you want to add some code to it.

Currently new constraints are created in form of x in Set(y) style where x is a variable or variables Set is the name of the constraint (more or less) and y can hold some other information needed. One Set is the AllDifferentSet and now the TableSet.

For different constraints you need different information for the table constraint we need words, mask, ...

Therefore we create a new ConstraintType in src/types.jl.

mutable struct TableConstraint <: Constraint
    std::ConstraintInternals
    current::RSparseBitSet
    ...
end

the std::ConstraintInternals has to exist as this will store information which every constraint needs like indices to know which constraints must be called if a variable changed.

Additionally you need to define the set which is the user interface. In this case we need TableSet(table) whereas in the AllDifferentSet case we don't have any parameters.

Nevertheless both need a dimension field which was the reason why you had to write x[1:9] in AllDifferentSet(9) before instead of x[1:9] in AllDifferentSet().

This means we as a developer need to define two types:

struct TableSetInternal <: MOI.AbstractVectorSet
    dimension :: Int
    table     :: Array{Int, 2}
end

struct TableSet <: JuMP.AbstractVectorSet 
    table     :: Array{Int, 2}
end
JuMP.moi_set(ts::TableSet, dim) = TableSetInternal(dim, ts.table)

The internal type is as the name suggests for the internal structure which holds the information of the dimension.

JuMP.moi_set(ts::TableSet, dim) = TableSetInternal(dim, ts.table)

is used to add the dimension.

How do we define the constraint?

All constraints are in src/MOI_wrapper/constraints.jl where you need to implement at least two functions:

First you need to tell that there is a new constraint which this solver supports:

MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorOfVariables},
    ::Type{TableSetInternal},
) = true

as mentioned before you now only need TableSetInternal.

Then you need to define the function add_constraint.

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::TableSetInternal,
)
    com = model.inner

    internals = ConstraintInternals(
        length(com.constraints) + 1,
        vars,
        set,
        Int[v.value for v in vars.variables]
    )

    constraint = TableConstraint(
        internals,
        ... # define the initials of all other fields
    )

    push!(com.constraints, constraint)
    for (i, ind) in enumerate(constraint.std.indices)
        push!(com.subscription[ind], constraint.std.idx)
    end
    # this is to show the number of constraints in each type. 
    # Go ahead and add yours in `types.jl` see: NumberConstraintTypes
    com.info.n_constraint_types.table += 1

    return MOI.ConstraintIndex{MOI.VectorOfVariables,TableSetInternal}(length(com.constraints))
end

The internals should be as mentioned here. Important is that the index is just one more than the previous and that vars, set and indices are defined.

push!(com.constraints, constraint)

pushes the new constraint to the list of constraints. Currently one has to fill in the subscription fields manually... I should change that ;)

In the end you have to return a MOI.ConstraintIndex.

Next up you have to define a still_feasible method in a new file in src/constraints/ with the following arguments:

function still_feasible(
    com::CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
    value::Int,
    index::Int,
)

and return true if setting index (represents the global index of the variable) to value has still a feasible solution and false otherwise.

That check is often quite easy. Harder is to prune the constraint:

function prune_constraint!(
    com::CS.CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal;
    logs = true,
)

Here you also have to return true if this is feasible and false if not but the idea is to prune as much as possible if it doesn't take tooo much time :)

How to remove a value?

You can remove a value from the search space with:

if !rm!(com, com.search_space[vidx], value)
    feasible = false
end

this removes value from variables[vidx]. If this makes the problem infeasible feasible is set to false. You can either return false immediately or do some stuff before you return if it's necessary for your data structure.

This can already be enough for basic constraints as it was for everything besides the table constraint.

There are a number of other methods which might need to be implemented by your constraint.

function init_constraint!(
    com::CS.CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
)

return true/false (feasibility)
end

Will be called before the first pruning.

function finished_pruning_constraint!(com::CS.CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal)

This function will be called after all pruning steps are done and before the next backtracking part starts. Here you can save the state at the end of a node which can be handy for backtracking.

function update_best_bound_constraint!(com::CS.CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
    var_idx::Int,
    lb::Int,
    ub::Int
)

It might be quite handy to provide some extra information to the lp solver if one exists i.e the range of the sum of variables.

This function gets one variable index and the lower und upper bound it will be set to. The goal here is change an lp variable which can be created in init to update the bounds of the variable (which might equal the sum) to give some hints to the lp solver.

You'll find more about this in the documentation (ones I'm done with it :D). At the moment you can check out the source code if you want to add something like this.

Currently there are three more functions which are there if you save some information in the constraint which depends on the node and might need to be reversed if we backtrack.

function single_reverse_pruning_constraint!(
    com::CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
    var_idx::Int,
    changes::Vector{Tuple{Symbol,Int,Int,Int}}
)

Will be called when reversing a node with the var_idx and a list of changes.

Each tuple holds the following information of what your constraint did in the pruning step.

# [1] :fix, :rm, :rm_below, :rm_above
# [2] To which value got it fixed, which value was removed, which value was the upper/lower bound
# [3] Only if fixed it saves the last ptr to revert the changes otherwise 0
# [4] How many values got removed (0 for fix)

i.e (:rm, 2, 0, 1) means that you removed the value 2. In this function you need to reverse this (or you can use the next function).

A more general function which gets called after single_reverse_pruning_constraint!:

function reverse_pruning_constraint!(
    com::CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
    backtrack_id::Int
)

There you can revert the whole backtrack_id in one go. If you are interested in the parent backtrack_id you can use:

parent_id = com.backtrack_vec[backtrack_id].parent_idx

And the last function:

function restore_pruning_constraint!(
    com::CoM,
    constraint::TableConstraint,
    fct::MOI.VectorOfVariables,
    set::TableSetInternal,
    prune_steps::Union{Int, Vector{Int}}
)

This gets called if we need to prune on some nodes which we already pruned before. This is the case if jump from one part of the tree to another probably because that node has a better bound.

prune_steps includes a list of all backtrack_ids on the path down to the one we want to reach. If you saved all the information you need in finished_pruning_constraint! you only need to care about the last element in prune_steps which you can get with last(prune_steps) which works also if prune_steps is an integer.

Thanks for reading and special thanks to my five patrons!

List of patrons paying more than 4$ per month:

• Site Wang

• Gurvesh Sanghera

• Szymon Bęczkowski

Currently I get 18$ 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).

Any questions, comments or improvements? Please comment here and/or write me an e-mail o.kroeger@opensourc.es and for code feel free to open a PR/issue on GitHub ConstraintSolver.jl

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

Want to be updated? Consider subscribing on Patreon for free