Series: Constraint Solver in Julia
- Part One: Setup and backtracking
- Part Two: Pruning
- Part Three: Pruning+Benchmarking
- Part Four: Alldifferent constraint
- Part Five: Better data structure
- Part Six: Sum constraint first part
This is part 7 of the series about: How to build a Constraint Programming solver?
Last time we introduced the equal sum constraint (
eq_sum) and figured out that it isn't that fast to solve a hard killer sudoku with it.
Edit 02.10.2019: I just found out that each cage (colored region) has a all different constraint as well. I don't use that here as I didn't read it :/ but it's also reasonable for the solver and this post that we concentrate on the sum constraint.
Last time we ended up with this if we don't backtrack:
- a killer sudoku has the same constraints as a normal sudoku
- plus the constraint that each colored region must add up to the number shown in the top left of that region
I talked about some ideas of how we might be able to further prune the search space:
- Combining sum constraint with all different constraint
- i.e bottom left 7,8,9 + 7,8,9 should add up to 16 but it is inside same all different => no 8 possible
- Use the sum constraint of row, column or block add up to 45
- bottom right block: We can have an additional constraint that the red marked cells add up to 9
First what is our current benchmark result for this problem? My setup:
- 6 core Intel(R) Xeon(R) CPU E5-1650 0 @ 3.20GHz
- Julia 1.2
julia> @benchmark main(;benchmark=true) BenchmarkTools.Trial: memory estimate: 44.83 GiB allocs estimate: 317941771 -------------- minimum time: 42.523 s (18.96% GC) median time: 42.523 s (18.96% GC) mean time: 42.523 s (18.96% GC) maximum time: 42.523 s (18.96% GC) -------------- samples: 1 evals/sample: 1
and this is the info:
Info: pre_backtrack_calls = 77 backtracked = true backtrack_counter = 87829 in_backtrack_calls = 4136981
Edit 02.10.2019: I've realized that I never showed the solution to the hard puzzle but I think it might be interesting. As far as I can see the solution does fulfill the actual constraints of killer sudokus (each colored region should have unique digits) as well.
Okay let's start with a very simple version:
If there are only two unfixed variables in a sum constraint we test all options and check whether it's feasible or not i.e if the possibilities are
[1,2,3,4] and the sum should equal 5 we currently don't prune anything as we normally only prune the outside values but we clearly can reduce the possibilities of the second variable to
[1,4] as well.
We check first how many unfixed variables there are and save the indices (both the normal index of the variable and the local index for
# if there are only two left check all options n_unfixed = 0 unfixed_ind = zeros(Int,2) unfixed_local_ind = zeros(Int,2) unfixed_rhs = constraint.rhs li = 0 for i in indices li += 1 if !isfixed(search_space[i]) n_unfixed += 1 if n_unfixed <= 2 unfixed_ind[n_unfixed] = i unfixed_local_ind[n_unfixed] = li end else unfixed_rhs -= value(search_space[i]) end end
Then we check whether we have
n_unfixed == 2 and if that is the case:
for v in 1:2 if v == 1 this, local_this = unfixed_ind, unfixed_local_ind other, local_other = unfixed_ind, unfixed_local_ind else other, local_other = unfixed_ind, unfixed_local_ind this, local_this = unfixed_ind, unfixed_local_ind end for val in values(search_space[this]) if !has(search_space[other], unfixed_rhs-val) rm!(search_space[this], val) changed[this] = true pruned[local_this] += 1 end end end
This means we iterate through all possibilities for the first variable and check whether the rest is a possibility in the second variable and then we do it the other way around. If the other variable doesn't have the possibility we can remove the possibility from the current checking variable.
This gives us a nice reduction but not enough:
julia> @benchmark main(;benchmark=true) BenchmarkTools.Trial: memory estimate: 33.59 GiB allocs estimate: 244525754 -------------- minimum time: 32.301 s (18.77% GC) median time: 32.301 s (18.77% GC) mean time: 32.301 s (18.77% GC) maximum time: 32.301 s (18.77% GC) -------------- samples: 1 evals/sample: 1
The next obvious step is to reduce it further by checking whether the two variables are also in the same all different constraint.
Then we can remove values from the
[7,8,9] + [7,8,9] == 16 case as the 8 obviously can't be chosen but only because of the combination with the all different constraint.
Therefore inside the
if n_unfixed == 2 we add:
intersect_cons = intersect(com.subscription[unfixed_ind], com.subscription[unfixed_ind]) is_all_different = false for constraint_idx in intersect_cons if nameof(com.constraints[constraint_idx].fct) == :all_different is_all_different = true break end end
This gets all constraints for both variables and checks whether there is an all different constraint which connects those two.
for val in values(search_space[this]) we add
if is_all_different && unfixed_rhs-val == val rm!(search_space[this], val) changed[this] = true pruned[local_this] += 1 if has(search_space[other], unfixed_rhs-val) rm!(search_space[other], val) changed[other] = true pruned[local_other] += 1 end end
as a safety check we check whether both indices still have possible values and if not we return our
if !feasible(search_space[unfixed_ind]) || !feasible(search_space[unfixed_ind]) return ConstraintOutput(false, changed, pruned) end
That gives us another nice speed up:
julia> @benchmark main(;benchmark=true) BenchmarkTools.Trial: memory estimate: 12.82 GiB allocs estimate: 94858954 -------------- minimum time: 12.512 s (17.94% GC) median time: 12.512 s (17.94% GC) mean time: 12.512 s (17.94% GC) maximum time: 12.512 s (17.94% GC) -------------- samples: 1 evals/sample: 1
We can of course do something like this for 3 variables or more but I don't want to bruteforce too much.
Instead I want to think more about the plot I've shown you before:
Here the two red cells should add up to 9 as the bottom right block (as well as all other blocks, rows, cols) adds up to 45 and the red cells plus the block adds up to \(54=13+15+11+15\).
In general I had the following idea:
- If we have all different constraints and sum constraints
- Check if we have #variables == #values for all different
- Compute the sum of all different
- Compute the sum using sum constraint where all variables are inside this all different constraint
- Compute the total sum of all sum constraints which have at least one variable inside all different
- If there are not too many outside variables add a sum constraint
- for the outside variables
- Take care that we don't mess up anything
At the start of solve I call:
# check for better constraints simplify!(com)
First step in that function:
- check if we have all_different and sum constraints
- where all different uses every value
b_all_different = false b_all_different_sum = false b_sum = false for constraint in com.constraints if nameof(constraint.fct) == :all_different b_all_different = true if length(constraint.indices) == length(constraint.pvals) b_all_different_sum = true end elseif nameof(constraint.fct) == :eq_sum b_sum = true end end
We add constraints during this process but to make it not too complicated we will only work on the current set of constraints.
if b_all_different_sum && b_sum n_constraints_before = length(com.constraints) for constraint_idx in 1:length(com.constraints) constraint = com.constraints[constraint_idx] if nameof(constraint.fct) == :all_different if length(constraint.indices) == length(constraint.pvals) ----> ADDING REST here (1) end end end end
During the following process we want to know the sum of the all different constraint (in our case always
sum(1:9) == 45), the sum of sum constraints where all the sums are inside the all different constraint (
in_sum) the total sum we compute which means
all_diff_sum plus the outside variables (
We need to keep track of which indices are outside and which indices we checked already.
all_diff_sum = sum(constraint.pvals) in_sum = 0 total_sum = 0 outside_indices = Int cons_indices_dict = arr2dict(constraint.indices)
If we used a variable already we will remove it from
cons_indices_dict such that we don't use any variable inside the all different constraint more than once.
for variable_idx in keys(cons_indices_dict) for sub_constraint_idx in com.subscription[variable_idx] ----> ADDING REST here (2) end end
We only check constraints for now which were already there before we ran
# don't mess with constraints added later on if sub_constraint_idx > n_constraints_before continue end
Now we want to check whether we are dealing with a
sub_constraint = com.constraints[sub_constraint_idx] if nameof(sub_constraint.fct) == :eq_sum ----> ADDING REST here (3) end
We definitely add the
total_sum and then check whether all of the indices of that constraint are inside the all different constraint. If that is the case we add the
rhs also to
in_sum and delete the variable from
cons_indices_dict such that we don't check the same constraint again later which would result in adding too much.
total_sum += sub_constraint.rhs all_inside = true for sub_variable_idx in sub_constraint.indices if !haskey(cons_indices_dict, sub_variable_idx) all_inside = false push!(outside_indices, sub_variable_idx) else delete!(cons_indices_dict, sub_variable_idx) end end if all_inside in_sum += sub_constraint.rhs end break
The break at the end is because we don't want to check other sub constraints if we found one already.
In Rest (1):
# make sure that there are not too many outside indices if length(outside_indices) < 3 add_constraint!(com, CS.eq_sum, com.search_space[outside_indices]; rhs=total_sum - all_diff_sum) end
We should check first whether it still works by running
This checks whether it gets solved appropriately which is the case. And this is our benchmark run:
julia> @benchmark main(;benchmark=true) BenchmarkTools.Trial: memory estimate: 3.25 GiB allocs estimate: 23988420 -------------- minimum time: 3.183 s (17.82% GC) median time: 3.198 s (18.05% GC) mean time: 3.198 s (18.05% GC) maximum time: 3.214 s (18.29% GC) -------------- samples: 2 evals/sample: 1
This looks much better than our starting point but still worse than I hoped it will look like.
Before we might try something else we should check that this actually works and doesn't mess things up. It's always good to think about problems which are different than the ones you currently try to solve.
What if there is a variable inside an all different constraint which is not part of a sum constraint? I know this is never the case in a killer sudoku but our solver should be a general solver so it should still be able to solve it.
The obvious problem is that we assume that there is a sum constraint for every variable in our all different constraint. What we do is basically saying that the white cells inside the bottom right block add up to 0 instead of we don't know what they add up to.
Therefore we need to check whether each variable is part of a sum constraint and if not we ignore that particular all different constraint.
Two more things for this post:
- Whenever we remove values
- Check if we removed all values or all but one:
- Check whether we are still feasible
- Make sure that we fix values when we can in the sum constraint
More on the second one later. I didn't think about it before and it's kind of a bug in the code...
We do the first part of checking inside
all_different but not inside
eq_sum and we should actually do it inside the
Each of the following functions will now return whether the function made the problem infeasible:
fix!, remove_above!, remove_below!, rm!
We need the
com model for checking this so the first parameter will be
com for all of those.
The problem is that we want to have
::CS.CoM to specify the type but the type doesn't exist as we have to specify the type of
Variable before it. I moved the type declaration into the main file now and include
mutable struct CoM.
If we fix a value we check whether this violates another constraint:
if !fulfills_constraints(com, v.idx, x) return false end
that is the first thing we check in the function and in the end we add a
In general I think it's reasonable to check first before removing anything. This means we have to compute in
remove_above! first how many values will be deleted to check whether only a single or no value is left and then check whether it would be feasible before we actually remove values.
This for example is the start of the
remove_below! function now:
vals = values(var) still_possible = filter(v -> v >= val, vals) if length(still_possible) == 0 return false, 0 elseif length(still_possible) == 1 if !fulfills_constraints(com, var.idx, still_possible) return false, 0 end end
Ah and we can set
bt_infeasible there as well i.e:
vals = values(var) still_possible = filter(v -> v >= val, vals) if length(still_possible) == 0 com.bt_infeasible[var.idx] += 1 return false, 0 elseif length(still_possible) == 1 if !fulfills_constraints(com, var.idx, still_possible) com.bt_infeasible[var.idx] += 1 return false, 0 end end
Later if we remove a value we call:
if !rm!(com, variable, value) # return INFEASIBLE end # otherwise add to prune
Now the next point is the following: if we have variables
a is fixed to
b is either
2 and the sum is 6.
We can easily fix the value of
1 but we currently don't. This is due to the case that \(5+2 \geq 6\) where 2 is the maximum value so it seems to be okay if there are other values which might be able to reduce the sum.
Therefore in the end instead of checking if there are two unfixed variables we also check if there is a single unfixed variable.
if n_unfixed == 1 if !has(search_space[unfixed_ind], unfixed_rhs) com.bt_infeasible[unfixed_ind] += 1 return ConstraintOutput(false, changed, pruned, pruned_below) else changed[unfixed_ind] = true still_feasible, pr_below, pr_above = fix!(com, search_space[unfixed_ind], unfixed_rhs) if !still_feasible return ConstraintOutput(false, changed, pruned, pruned_below) end pruned[unfixed_local_ind] += pr_above pruned_below[unfixed_local_ind] += pr_below end elseif n_unfixed == 2 # ...
This is the first time that we actually fix a variable directly. Hope you can remember that if we fix a variable we move two pointers:
Therefore in reverse pruning we also need to move the
I created a new
pruned_below vector to save when we pruned from the other direction.
I'm not too sure whether there might be better ways to compute the next variable to branch on instead of the
bt_infeasible counter. For this problem it doesn't seem to be the ideal choice but still better than the normal branch on lowest number of possibilities one.
julia> @benchmark main(;benchmark=true, backtrack=true) BenchmarkTools.Trial: memory estimate: 3.57 GiB allocs estimate: 27452510 -------------- minimum time: 3.421 s (16.52% GC) median time: 3.489 s (17.88% GC) mean time: 3.489 s (17.88% GC) maximum time: 3.556 s (19.19% GC) -------------- samples: 2 evals/sample: 1
That is the end result for the hard killer sudoku now but it's quite good that for the easy one from wikipedia we don't need backtracking ;)
I've added some test cases only for the sum constraint to make sure that it does what it should. This way I observed the problem with not fixing a variable.
The last part solves issue #7 A variable can't be fixed in a constraint
What's up next time?
Good question. I'm thinking about:
- Graph coloring
- scheduling problems (check feasibility first)
- getting all solutions
- being able to optimize (will be a harder one based on my current DFS implementation :/)
- Element constraint
Maybe you have some ideas: Please comment!
Stay tuned ;)
A new post is out: Part 8: UI refactoring and bugfix of sum constraint
And as always:
Thanks for reading and special thanks to my first two patrons!
If you want to have early access to the next posts as well and 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 email@example.com 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
|Become a Patron!|