Creation date: 2019-09-16

Series: Constraint Solver in Julia

This is part 5 of the series about: How to build a Constraint Programming solver?

Let's do a small recap of what we have done so far. We've created a julia package in the first post and implemented the most basic backtracking for solving a sudoku puzzle. We understood that this is pretty slow as we could use the information after we fixed a cell to a value to find out the values of cells which are connected by a constraint to that cell. In part 3 I described that these new fixes itself can be used further on which gave us a performance improvement again.

In the last post we tried not to make backtracking faster but instead tried to avoid it as much as possible by solving the sudoku in a more human like fashion.

We could now further try to reduce the number of steps we call the now slower `alldifferent`

constraint or make it faster or change the ordering of the calls or fix values in a different ordering, as at the moment we only decide which cell we want to tackle next in backtracking but not on the ordering of how we fix the number in a cell i.e if it has the possible values 5,6,7 it might be reasonable to start with the 7 if we think that it's more likely based on some heuristics.

What I want to do instead in this post is changing the data structure a bit. We are currently working on a grid which seems to be a good choice for a sudoku and the 8-queen problem as well as magic squares but it's not always the case for example it's stupid for graph coloring so maybe we should just have linear indices internally instead of cartesian indices. Another point is that we are working with a dictionary to represent the values in our search space for each cell but adding and removing to a dictionary sometimes seems quite stupid as we set a cell in backtracking and then reset it and have to fill all the values again it had before. Maybe there is a better data structure. Besides that it seems to be reasonable to have a `struct`

for a variable which can hold some information of that variable maybe something like a name for a country in graph coloring?

Additionally I think it's quite weird to have a grid and a search space separately if I think about it. The variable could hold the information about being fixed or not.

I found a new idea at MiniCP Part 2 (unfortunately all sites of the website have the same url :/). MiniCP is a constraint solver used in teaching constraint solving so you might in general be interested at that. Pascal Van Hentenryck who also is the professor of the online about Discrete optimization at coursera which I did three years ago works on that. Definitely check out his course!

Back to the idea: Instead of having a dictionary we hold the possibilities in a vector.

`[1,2,3,4,5,6,7,8,9]`

now removing from a vector is much slower than from a dictionary so we have to think about that. In the reverse pruning part we would also add values back to the vector which I think is a bit stupid in the dictionary version as we kind of removing some values all the time and then want to add them again if we made a mistake.

Instead we could use a pointer which tells us how many values are still possible i.e at the beginning we say we have 9 possible values and if we remove the 9 we would have

`[1,2,3,4,5,6,7,8,9], 8`

where the 8 indicates that we only have 8 possible values.

Okay sure super easy if we always remove the 9 first and then 8 and then 7 and so on.

Obviously it's not so easy if we want to remove the 3 now.

Let's start with our 9 values again and we have two vectors now. The values (`values`

) shown in the first row and the indices (`idxs`

) shown in the second row.

If we want to know at which position the value 7 is, we would call `idxs[7]`

and we would get `7`

back as the value `7`

is at position `7`

in `values`

.

What happens now if we want to remove the 3?

We have the red barrier which tells us that only the left part is still possible and we swapped the value which was behind that barrier but is still possible `9`

with the value we want to remove `3`

.

Now we need to figure out the changes in the indices such that we still know where the value `3`

is as we later want to have an easy check whether it's still possible or not.

So if we want to check whether the 3 is a possible value we would do `idxs[3]`

get 9 back and would check whether 9 is behind that barrier and see: `3`

is not a possible value.

Pretty simple, right? Let's remove value `5`

:

Straightforward I would assume. It gets a bit more tricky when we want to remove a value which isn't as it's initial position anymore i.e 9.

If we remove 9:

Before this the `values`

looked exactly the same as the `idxs`

but here that wouldn't work anymore but besides that we did what we did before as well.

We moved the red barrier then set the value which we want to remove behind that barrier and by move I mean swap its position such that we still have the values 1-9.

Therefore we had to look at which position the `9`

is which was at position `3`

and change the `idxs`

vector in such a way that at the position where we now have the `9`

(position `7`

) points to the position where the 9 was as we swapped the values.

Setting the variable to a value can be done in a similar fashion by moving the barrier to directly behind the first value and then swapping the desired value to the front.

This is what they do with MiniCP and I have a little problem with that as in backtracking we set a value and later set it to something else which is not easy using this method as we move the barrier more than one step.

My approach instead is to have two barriers one from left and one from the right. This might be a stupid idea and I see myself hating it in the future but at the moment it seems interesting so I'm sorry I do it my way...

We create a `struct`

for our `Variable`

in `src/Variable.jl`

```
mutable struct Variable
idx :: Int
first_ptr :: Int
last_ptr :: Int
values :: Vector{Int}
indices :: Vector{Int}
end
```

Now we need some functions like:

`has`

, `rm!`

, `fix!`

We first assume that our values are an interval starting with 1 so we need later something which works with different intervals but **later**.

```
function has(v::Variable, x::Int)
if x > length(v.values) || x < 1
return false
end
ind = v.indices[x]
return v.first_ptr <= ind <= v.last_ptr
end
```

Then our remove function which I explained at the beginning. Here we don't check whether the value is actually inside the range for speed so whenever we remove a value we should be sure that it isn't already removed.

```
function rm!(v::Variable, x::Int)
ind = v.indices[x]
v.indices[x], v.indices[v.values[v.last_ptr]] = v.indices[v.values[v.last_ptr]], v.indices[x]
v.values[ind], v.values[v.last_ptr] = v.values[v.last_ptr], v.values[ind]
v.last_ptr -= 1
end
```

here we also just set the barriers without checking whether it makes sense because later in backtracking we fix a value several times which normally isn't done in constraint programming as if a value is set it should be set.

```
function fix!(v::Variable, x::Int)
ind = v.indices[x]
v.last_ptr = ind
v.first_ptr = ind
end
```

we should also have a function which checks whether a variable is fixed:

```
function isfixed(v::Variable)
return v.last_ptr == v.first_ptr
end
```

Okay now to the actual changes:

We want to have only `search_space`

and not `grid`

and `search_space`

and we build the `search_space`

by simply adding more variables to it.

The internal structure is a simple vector then and the user can have any structure for example a grid as before for the sudoku.

I know these are substantial changes but as mentioned I'm writing this along coding and therefore some bugs and bigger changes will happen that's somehow also the goal of the series to see the changes while they happen instead of following a step by step guide on how to build a constraint solver where everything works :D

Okay we need a `addVar!`

function which adds it to our model and returns the `Variable`

object back to the user who can use it in it's own data structure.

I'm thinking about:

```
com_grid = Array{CS.Variable, 2}(undef, 9, 9)
for (ind,val) in enumerate(grid)
if val == 0
com_grid[ind] = CS.addVar!(com, 1, 9)
else
com_grid[ind] = CS.addVar!(com, 1, 9; fix=val)
end
end
```

we provide a start and end value and can fix it while creating the variable.

```
function addVar!(com::CS.CoM, from::Int, to::Int; fix=nothing)
ind = length(com.search_space)+1
var = Variable(ind, 1, to-from+1, 1:to-from+1, 1:to-from+1)
if fix !== nothing
fix!(var, fix)
end
push!(com.search_space, var)
push!(com.subscription, Int[])
return var
end
```

I also changed the `init`

function as we build the search space on the go now:

```
function init()
com = CoM()
com.constraints = Vector{Constraint}()
com.subscription = Vector{Vector{Int}}()
com.search_space = Vector{Variable}()
com.bt_infeasible = Vector{Int}()
com.info = CSInfo(0, false, 0, 0)
return com
end
```

For the constraint functions we now not just give the cartesian indices we give a list of variables over which the constraint should hold i.e `CS.add_constraint!(com, CS.all_different, com_grid[CartesianIndices((rc:rc,1:9))])`

instead of `CS.add_constraint!(com, CS.all_different, CartesianIndices((rc:rc,1:9)))`

which means that this is our new `add_constraint!`

function:

```
"""
add_constraint!(com::CS.CoM, fct, variables)
Add a constraint using a function name and the variables over which the constraint should hold
"""
function add_constraint!(com::CS.CoM, fct, variables)
current_constraint_number = length(com.constraints)+1
indices = vec([v.idx for v in variables])
constraint = Constraint(current_constraint_number, fct, indices)
push!(com.constraints, constraint)
for i in indices
push!(com.subscription[i], current_constraint_number)
end
end
```

**Attention:** This will be changed later based on other problems ;)

Now to a more interesting part as this is unfortunately quite boring to read I assume :/ Let's have a look at the `all_different`

function.

Normally we saved which indices changed, how we pruned (which values got removed) and which cells got fixed. Now as we don't really remove any values we simply move a pointer and do some swapping in the reverse pruning step we can just move the `last_ptr`

to the right (don't need to swap as the ordering is not important) and therefore we only need to save how many values we have pruned but not which ones. We also don't have to save which values we fixed which we only needed in reverse pruning.

For now I will still use a dictionary for `changed`

as I think it's easier to handle.

This at the beginning of `all_different`

```
changed = Dict{CartesianIndex, Bool}()
pruned = Dict{CartesianIndex,Vector{Int}}()
fixed = Dict{CartesianIndex, Bool}()
```

gets changed to:

```
changed = Dict{Int, Bool}()
pruned = zeros(Int, length(indices))
```

Then I changed the `fixed_vs_unfixed`

function which had as input the `grid`

which doesn't exist anymore and `not_val`

which is now only part of the user interface but not part of our solver anymore.

There was

```
if grid[i] != not_val
push!(fixed_vals, grid[i])
```

which is now

```
if isfixed(search_space[i])
push!(fixed_vals, value(search_space[i]))
```

which is much easier to read well we have to write the `value`

function :D

```
function value(v::Variable)
@assert v.last_ptr == v.first_ptr
return v.values[v.last_ptr]
end
```

but that was simple. I check here that it's actually fixed again which we might want to remove when we are sure that everything works.

I encountered some "problems" during the changing process that I want to share and not show you every change I made as that might become boring.

We don't have `pval`

anymore so I decided that each constraint gets a range of values which it operates on instead and `all_different`

now gets the `Constraint`

object and not just the indices. That makes the `add_constraint!`

function a little bit bigger see GitHub

`reverse_pruning`

gets quite short now as we don't have to check anymore which values we removed or fixed in our pruning step. We just have to increase the `last_ptr`

.

So this:

```
for local_ind in constraint.indices
if !haskey(constraint_output.pruned, local_ind)
continue
end
if haskey(constraint_output.fixed, local_ind)
com.grid[local_ind] = com.not_val
end
if !haskey(com.search_space, local_ind) && length(constraint_output.pruned[local_ind]) > 0
com.search_space[local_ind] = Dict{Int,Bool}()
end
for pval in constraint_output.pruned[local_ind]
com.search_space[local_ind][pval] = true
end
end
```

is now:

```
for (i,local_ind) in enumerate(constraint.indices)
com.search_space[local_ind].last_ptr += constraint_output.pruned[i]
end
```

I've added an `issue`

to GitHub that we need to be careful if a constraint directly fixes a value and not removes all other values like we do now as then this would not work. Then we also need to move the `first_ptr`

of the variable accordingly.

In our backtracking function we now save the position of `last_ptr`

before we set a value and in the end before `return :Infeasible`

we reset it to what it was before.

We have the `fulfills_constraints`

function which we call before setting a value in `all_different`

but as we don't have the `search_space & grid`

structure anymore it isn't fulfilled as we have set it already. Therefore the function now checks every index which is not the current one.

Checking whether our problem is solved is done by:

```
if all(v->isfixed(v), com.search_space)
return :Solved
end
```

now which takes a bit time compared to `length(com.search_space) == 0`

but it's okay. We don't call it often.

You remember that I said we care about variables with a range different to `1-x`

**later**? That's now:

So we simply store an offset in our `Variable`

struct. If we have `0-8`

the offset is 1. Which means if we want to remove 5. We look at index `5+offset = 6`

of the indices vector to get the index of the value where 5 is stored. That's it.

You might think:

That was hard work

We should have done that before

Well yes that might be correct but I think we kind of went from a somehow specialized solver to a more abstract solver and we improved the performance ;)

We moved from light blue to dark blue which isn't too bad.

As mentioned in the last post our `all_different`

constraint only works if we have the same number of indices as number of variables but how about if we need to choose 9 values from 15?

We still check the maximum matching against the number of indices (`9`

) but we need to find

an even alternating path starting at a free vertex

according to Berge's lemma. A free vertex is now a value which wasn't chosen in the maximum matching. Now we could just try to find all even alternating path starting at a free vertex but that doesn't sound practical so instead we add another node to and connect it in such a way that we can still find all strongly connected components as we do anyway.

As an example:

The upper blue lines would be the even alternating path starting at the free vertex `10`

and the two blue lines at the bottom make it to a cycle which will be detected by our strongly connected components.

We don't care about these extra edges when we remove values so we need to check there whether the edge has an end point at the new vertex and if that's the case we ignore it.

I've removed the test case with the 42,44,... sudoku as I think it's not very helpful and we could still solve it by having all a range `42-x`

and remove the odd values `43,45`

at the start.

I've added a test which has a normal sudoku but some cells don't have a range from 1-9 but `9-11`

or something. This adds more possibilities to solving the sudoku but adds functionality we didn't have before with our `build_search_space`

structure.

I'm thinking about solving Killer sudokus next which works with a `sum`

constraint as well as the `all_different`

constraint.

Maybe backtracking shouldn't be implemented as recursion to never run in the problem of getting a recursion limit error.

Should the variable check whether it's feasible to set it? So whenever a variable is fixed it might should call `fulfills_constraints`

itself and `return`

whether it's successful as I think the person who implements constraints shouldn't handle that this directly.

Stay tuned ;)

**Here it is:**

And as always:

Thanks for reading!

If you enjoy the blog in general please consider a donation via Patreon to keep this project running.

If you have any questions, comments or improvements on the way I code or write please comment here and or write me a mail o.kroeger@opensourc.es and for code feel free to open a PR on GitHub ConstraintSolver.jl

Additionally if you want to practice code reviews or something like that I would be more than happy to send you what I have before I publish. ;)

I'll keep you updated if there are any news on this 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.

© Ole Kröger. Last modified: February 27, 2022. Website built with Franklin.jl.