Constraint Solver Part 8: UI Refactoring

Creation date: 2019-10-04

Tags: julia, constraint-programming, optimization

Series: Constraint Solver in Julia

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

I published 7 posts on this series which started last month hope to work on it further for a while. Either in the next part or part 10 I think I'll do a recap post which goes through all the parts we've done again. Maybe a recap video which will show the code. Not sure yet. If you have some ideas please comment!

Done: Recap video and post

Besides that I yesterday saw a tweet by Juan Pablo Vielma who is an Associate Professor of Operations Research at MIT and is active in the Julia community as well.

He also mentioned in a comment that he is currently looking into integrating Google OR-tools into JuMP.jl which is the modeling language for mathematical optimization in Julia.

That would be very interesting. Nevertheless I will of course keep working on this as this is for educational reasons (also like how not do some things :D).

Maybe it's also a good time to mention here why I'm writing these posts the way I do again:

  1. I think reading blogs where everything works perfectly because the blogger is a specialist in his career doesn't get you far, similar to why Veritasium asked people questions in his videos not just showing the answer.

  2. Explaining while you're still on a similar level as your reader is far easier as explained also by Wait But Why? who is a layman in the topics before writing about them

  3. I'm writing better code if I have the feeling that I'm watched :D I find bugs before I publish (and after :D ) after thinking about it because I have to write it down in a reasonable manner.

  4. You should not take everything you read as if it is the truth. Definitely not on my blog as there are probably way better methods to write a simplex solver or now a constraint programming solver

  5. Hopefully I can bring you to the stage where you try out coding projects like this by yourself. Often in university or school you often start with: Write a fibonacci or prime function but sometimes it doesn't get much more difficult than that.

Try out some real problems and make them open source!

Okay so what is this post about?

Sometimes you start like I did with the sudoku puzzle and then want to generalize more and sometimes don't quite get your actual goal of the project besides enjoying programming and writing about it ...

I introduced the sum constraint and realized that we now need a right hand side (rhs) which we don't need for the alldifferent constraint. Therefore it defaults to 0. Then later we will have constraints like a == b or a != b where a and b are variables. We would make a new constraint for equal and not_equal or something and would write it like:

CS.add_constraint!(com, CS.equal, a; rhs=b)

maybe but remember that rhs is an integer now so we would need to make it Union{Integer,Variable} or I don't know. It seems like this isn't the way to do.

What would be the ideal way of writing it as the user?

Probably similar to the way it's done it JuMP.jl:

CS.add_constraint!(com, a == b)
CS.add_constraint!(com, a != b)
CS.add_constraint!(com, all_different([a,b]))
CS.add_constraint!(com, sum([a,b]) == 5)

and a current limitation is that we can't have coefficients for a sum constraint which seems to be kind of strict to me as it shouldn't be too complicated to include it into the current implementation. Such that:

CS.add_constraint!(com, 5*a-2*b == 5)

is possible as well and 

CS.add_constraint!(com, 5*a-2*b-5+b==c)

can be transformed as well but the user shouldn't need to care about it.

Now each constraint needs an id and indices such that it gets called if a variable gets changed that the constraint subscribed to.

Each constraint can produce the same output of what changed but different input information.

Later: I think it might be reasonable to save exactly what changed such that we can have a way to reproduce everything in a tree like structure such that optimization later becomes easier i.e. find the sudoku solution where some values add up to the lowest value. For now we only care about the first solution and don't optimize anything (besides performance ;) )

What do we have to do?

The add_constraint! function will only get two inputs now. The model and a constraint of some type. We will have different constraint types for all_different, sum, equal or not_equal. I would call the sum constraint type LinearConstraint as we have some kind of linear combination. It will have a rhs as well as later an operator to support \(\leq\) and \(\geq\) as well.

They all should be a constraint type so we need some type of inheritance. Julia uses so called abstract types for that.

abstract type Constraint
end

It doesn't really hold any information but other types can say they are a subtype and we can write constraint::Constraint instead of constraint::Union{AllDifferentConstraint,LinearConstraint,DirectConstraint} or something like that. We could of course just use constraint instead of any type but it is reasonable to specify what this function expects as an input. Now AllDifferentConstraint is actually the most basic constraint as it only holds the idx, fct, indices and pvals which I think we need for all constraints. Therefore I call it BasicConstraint.

That is easier than the sum constraint so we start with it:

BasicConstraint

abstract type Constraint 
end

mutable struct BasicConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    BasicConstraint() = new()
end

We need the BasicConstraint() = new() as we don't know the index at the beginning when we create it. This has to be filled by the ConstraintSolver.jl but the others will be filled in the corresponding function file like all_different.jl. Same is true for pvals.

I also exclude eq_sum.jl for now and use my current.jl function for testing.

At the moment we have a create_constraint function which fills pvals but maybe we should just call it like that similar as we do with set_in_all_different!. So we create set_pvals!(com::CS.CoM, constraint::Constraint)

Then we have a add_constraint! function:

function add_constraint!(com::CS.CoM, constraint::Constraint)
    constraint.idx = length(com.constraints)+1
    push!(com.constraints, constraint)
    set_pvals!(constraint)
    for (i,ind) in enumerate(constraint.indices)
        push!(com.subscription[ind], constraint.idx)
    end
end

Very basic so far. Now it doesn't work of course. We need to build a constraint. Easy thing first: Reformulate the constraint from the user side:

This were the constraints before:

function add_sudoku_constr!(com, grid)
    for rc=1:9
        #row
        CS.add_constraint!(com, CS.all_different, grid[CartesianIndices((rc:rc,1:9))])
        #col
        CS.add_constraint!(com, CS.all_different, grid[CartesianIndices((1:9,rc:rc))])
    end

    for br=0:2
        for bc=0:2
            CS.add_constraint!(com, CS.all_different, grid[CartesianIndices((br*3+1:(br+1)*3,bc*3+1:(bc+1)*3))])
        end
    end
end

Now it looks like this:

function add_sudoku_constr!(com, grid)
    for rc=1:9
        #row
        variables = grid[CartesianIndices((rc:rc,1:9))]
        CS.add_constraint!(com, CS.all_different([variables...]))
        #col
        variables = grid[CartesianIndices((1:9,rc:rc))]
        CS.add_constraint!(com, CS.all_different([variables...]))
    end

    for br=0:2
        for bc=0:2
            variables = grid[CartesianIndices((br*3+1:(br+1)*3,bc*3+1:(bc+1)*3))]
            CS.add_constraint!(com, CS.all_different([variables...]))
        end
    end
end

I want to have a vector of variables but the uses stores it in a 2D array so [variables...] flattens it to a vector.

which is not defined yet.

Inside all_different.jl:

function all_different(variables::Vector{Variable})
    constraint = BasicConstraint()
    constraint.fct = all_different
    constraint.indices = Int[v.idx for v in variables]
    return constraint
end

Now let's try it out:

include("test/current.jl)
main(;benchmark=false)

Okay first problem:

ERROR: type BasicConstraint has no field in_all_different

Well it makes sense that the all_different constraint doesn't have that so we add a constraint that only if the constraint type has the field in_all_different we make the check.

if :in_all_different in fieldnames(typeof(constraint))

If we only work in Julia 1.2 we could use hasfield but I want to support earlier versions for now.

That's it for BasicConstraint. We will later make a benchmark check whether this is slower somehow.

LinearConstraint

The idea is that we have a function for + between two or more variables and later support coefficients as well such that we can have constraints like 2x+4y == 7. That will be saved as a LinearConstraint and the first part (before ==) is of type LinearVariables.

mutable struct LinearVariables
    indices             :: Vector{Int}
    coeffs              :: Vector{Int}
end

This just saves the id of each variable and the corresponding coefficient. I decided to only support integers. The full constraint type:

mutable struct LinearConstraint <: Constraint
    idx                 :: Int
    fct                 :: Function
    indices             :: Vector{Int}
    pvals               :: Vector{Int}
    coeffs              :: Vector{Int}
    operator            :: Symbol
    rhs                 :: Int
    in_all_different    :: Bool
    LinearConstraint() = new()
end

The operator here is :(==) for the equal case.

I've created a linearcombination.jl which holds the function to create LinearVariables:

Most basic one is:

function Base.:+(x::Variable, y::Variable)
    return LinearVariables([x.idx,y.idx],[1,1])
end

That way you can define a + function in Julia for your own types. This returns a LinearVariables such that we also need a function to sum an already existing LinearVariables with a variable:

function Base.:+(x::LinearVariables, y::Variable)
    lv = LinearVariables(x.indices, x.coeffs)
    push!(lv.indices, y.idx)
    push!(lv.coeffs, 1)
    return lv
end

You probably get the idea. We also need a function to create a LinearVariables when we have 2*x which is done analogously and a function which sums two those types directly.

Then in the eq_sum.jl file we actually have to build the constraint so combine it with == RHS:

function Base.:(==)(x::LinearVariables, y::Int)
    lc = LinearConstraint()
    lc.fct = eq_sum
    lc.indices = x.indices
    lc.coeffs = x.coeffs
    lc.operator = :(==)
    lc.rhs = y
    return lc
end

and I changed:

function eq_sum(com::CS.CoM, constraint::Constraint; logs = true)

to 

function eq_sum(com::CS.CoM, constraint::LinearConstraint; logs = true)

as well as the formulation of add_constraint! in all test cases i.e

CS.add_constraint!(com, sum([com_grid[CartesianIndex(ind)] for ind in s.indices]) == s.result)

Now we have only a support for coefficients that equal 1 which is the reason why I added:

if any(c->c!=1, constraint.coeffs)
    throw(ArgumentError("Currently only a coefficient of 1 is possible in sum constraints"))
end

Okay maybe we should just add that functionality :D

Wait a second I kinda made a mistake in the sum constraint :(

Bugfix of sum constraint

I wrote:

if maxs[i] < pre_maxs[i]
    still_feasible, nremoved = remove_below!(com, search_space[idx], maxs[i])

which reads to a person who doesn't know the code before: If the maximum is lower than previously we remove values below the new maximum which of course is garbage.

It has to be: above the new maximum but I've somehow swapped maximum and minimum before therefore it worked. Weird!

Sorry for the small update here:

Let's check it again: We have [1,2] + [5,9] == 6 this should result in [1]+[5] so:

We have to put the rhs to the left: [1,2] + [5,9] - 6 == 0. Then we compute the full minimum and the full maximum:

full_min = 1+5-6 == 0
full_max = 2+9-6 == 5

Now we remove our current minimum (minimum of the first variable) from the overall minimum and get: 0-1 = -1. Then we check if this is feasible and it is because using our range from [1,2] we can raise it to 0 or higher. Important is that we can raise it to 0!

The next check is: What is the maximum we can reach without the current index.

That is \(5-9 = -4\) which makes sense as the maximum 2 minus the rhs 6 is \(-4\).

If this value plus our maximum is lower than 0. We can't reach the sum even if we use all our max values.

That's the stuff about the feasibility check I've did so far but we don't have to do that every time as it is actually just a check whether full max is \(\geq 0\) and full min is \(\leq 0\). Therefore we only need to compute it once. I really didn't think that through.

Okay to the actual thing: We want to prune values...

We have c_min = 0-1 = -1 and c_max = 5-2 = 3 now actually the c_min gives us information about the maximum value we can use for the current index which is just --1 = 1 because if we go higher than that we can use all minimal values afterwards but still can't reach our goal low value 0.

This means we can remove values above 1.

To the second part we look at variable #2:

c_min = 0-5 = -5 and c_max = 5-9 = -4

c_max means that the maximum we can achieve without using the variable is \(-4\) so the minimum of the variable is \(4\) but we already know that is bigger than \(5\). On the other hand the c_min gives us again that the new maximum is \(5\) which helps.

This fixes everything to 1+5 == 6.

Now it seems to make sense. Sorry for that.

Coefficients

The simple idea is to work with the multiplied variants first so multiply it in the max value and in the min value sums and later divide it.

Whenever we have something like 7*x + y == 50 where y is between 1 and 5 for example we would get y = 1 and 7*x = 49 which means that we now have to divide 49 by 7 to get x. Same is true for remove_below and remove_above and also needs to be done in the feasibility check function of eq_sum.

We also have to check in our simplify function that the coefficients are 1 because otherwise we can't combine our all different constraint with the sum constraint.

Wait a second we need a talk a bit about negative coefficients ;)

So if a coefficient is negative we need to compute min and max the other way around and we need to swap remove_below with remove_above.

I've added some tests for this and also tests for using negative variables to be sure that it works as well.

What else?

It's October which means you can win a Hacktoberfest T-Shirt by digitalocean here.

All you have to do is create four pull requests on GitHub. I would of course appreciate it if you help me with this project if you found a bug or want to add a feature. You can of course contribute to other open source projects ;)

I've created a Patreon page a while ago and I'm very happy to have TWO patrons already!

It costs 1$ per month for getting early access to my new posts. Currently it is 2 days before. If you enjoy this series and want to be up to date consider subscribing on Patreon. If you think it's not worth it you can drop off before the end of this month and you don't pay anything.

Next:

Next post: First recap

I hope that the Kaggle Santa challenge will come soon. My current idea is to share it with you while I try to win :D Last time my team and I achieved quite a good score but in the end it was hard to actually write down what we did and why and I really like the current series. Therefore I might share it on the go. The problem is that I'll not be able to get a good score if I publish my code... as others can simply copy it and improve but I think it sounds good as I really like to read other peoples code when they share it and quite some do (maybe not their best version).

Another problem is that I can't share it with my patrons earlier as I have to share it publicly or only with my team regarding the rules by Kaggle which obviously make sense.

If I'm able to write on this solver during that time I think it's okay and otherwise I'll pause my billing cycle in that time frame on Patreon.

Stay tuned ;)

And as always:

Thanks for reading and special thanks to my first three patrons!

List of patrons paying more than 2$ per month:

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 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 on Patreon for free
Subscribe to RSS
© Ole Kröger. Last modified: June 09, 2020. Website built with Franklin.jl.