Constraint Solver: How to use JuMP for your own solver?

Creation date: 2020-01-10

Tags: julia, optimization, constraint-programming

First of all: Happy new year to everyone!

This is kind of part 15 of the series about: How to build a Constraint Programming solver? But it is also interesting for you if you don't care about constraint programming. You can learn how to use JuMP.jl as the modelling layer for your own solver.

Anyway these are the other posts so far in the ongoing series:

All of that can only happen because of the generous support of some people. I'm very happy that I can cover my server expenses now because of my Patrons. To keep this going please consider supporting this blog on Patreon. If you do I will create a better/faster blog experience and you will get posts like this earlier! Everything will be free forever but if you have some cash to spare -> This might be an option ;)

Special limited offer... I either just lost your attention or I have it now :D I decided to offer a 1:1 (binary) talk where you can ask me about my work and just chat with me if you're willing to join the Patron 4$ club. ;) This offer is until the 10th of February.

If you want to just keep reading -> Here you go ;)

Maybe I should have done that from the ground up but wasn't sure how to and probably it was easier for the blog as well to start without a big code base. Therefore here it is now: Using JuMP.jl as the modelling layer for the constraint solver I'm creating for the past couple of months.

UserInterface change

Before I explain how it's done I want to show how the user interface changes first:

How it was before

com = CS.init() # creating a constraint solver model
com_grid = Array{CS.Variable, 2}(undef, 9, 9)
for (ind,val) in enumerate(grid)
    if val == 0
        com_grid[ind] = add_var!(com, 1, 9)
    else
        com_grid[ind] = add_var!(com, 1, 9; fix=val)
    end
end

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

which is a bit of the sudoku part. Now it looks like this:

m = Model(with_optimizer(CS.Optimizer))
@variable(m, 1 <= x[1:9,1:9] <= 9, Int)
# set variables
for r=1:9, c=1:9
    if grid[r,c] != 0
        @constraint(m, x[r,c] == grid[r,c])
    end
end
for rc = 1:9
    @constraint(m, x[rc,:] in CS.AllDifferentSet(9))
    @constraint(m, x[:,rc] in CS.AllDifferentSet(9))
end

I think it looks a lot cleaner now but more importantly it is the standard interface for mathematical optimization solvers written in Julia. This also means that if there will exist other constraint solvers in the future one can easily swap the solver. There is currently a catch because of CS.AllDifferentSet(9) but if there is more support for constraint solvers in the future this will be easier.

What needs to be done?

It is useful to create a directory in src/ called MOI_wrapper for all these interface things. You might want to check the MathOptInterface.jl documentation as well, which is kind of the lower level version of JuMP.

I've created five files in this directory

I'll go over the basics a bit like the general setup and then some complications I had for things which are constraint solver specific. Additionally in the end there will be a short section on a better structure of pruning which just removes unnecessary code blocks and removes some pruning steps which are not needed.

General stuff

We need some Optimizer struct:

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{CoM, Nothing}
    variable_info::Vector{Variable}
    # which variable index, (:leq,:geq,:eq,:Int,:Bin), and lower and upper bound
    var_constraints::Vector{Tuple{Int64,Symbol,Int64,Int64}} 
    status::MOI.TerminationStatusCode
    options::SolverOptions
end

that holds the complete ConstraintSolverModel as well as some info about the variables, the status and the solver options as well as a function which gets called when the solver gets added to the model.

function Optimizer(;options...) 
    com = CS.init()
    options = combine_options(options)
    return Optimizer(
        com, 
        [], 
        [],
        MOI.OPTIMIZE_NOT_CALLED,
        options
    )
end

then there is the function MOI.optimize!(model::Optimizer) which gets called when we want to optimize/solve the model. Similar to our previous solve! call. It will now be optimize!(m).

A bit about the @variable stuff:

In src/MOI_wrapper/variables.jl we need to define first what is supported. For example Integer and Binary variables or only real variables. This is normally for stuff like solving a linear program or mixed integer linear solvers.

This gets written down like:

const SVF = MOI.SingleVariable
MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.LessThan{Float64}}) = true

and the corresponding implementation for it:

function MOI.add_constraint(model::Optimizer, v::SVF, lt::MOI.LessThan{Float64})
    # adding info to variable info & var_constraints
    return MOI.ConstraintIndex{SVF, MOI.LessThan{Float64}}(constraint_index)
end

we also need a function for a simple @variable(m, x) (without constraints)

function MOI.add_variable(model::Optimizer)
    # adding to variable_info
    return MOI.VariableIndex(index)
end

It turns out that we also need a add_variables function as otherwise @variable(m, x[1:9]) would not work:

MOI.add_variables(model::Optimizer, n::Int) = [MOI.add_variable(model) for i in 1:n]

Some minor comments: There is a special case for interval variables like @variable(m, 1 <= x <= 8) which needs to be implemented if 

MOI.supports_constraint(::Optimizer, ::Type{SVF}, ::Type{MOI.Interval{Float64}}) = true

is set but it is also supported if we don't have the supports_constraint for Interval but instead only for LessThan and GreaterThan then it gets bridged.

Similar things are done for the normal constraints like a+b==2 which we use in the killer sudoku part

const = MOI.ScalarAffineFunction{Float64}
MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{MOI.EqualTo{Float64}}) = true

We then formulate everything into our LinearConstraint. Using JuMP/MOI we also don't have to care about special cases like when someone writes a+2-4+b == 12. Reformulation is done automatically.

To access the values of variables or the objective value as well as the status we need to implement:

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
	return model.status
end

function MOI.get(model::Optimizer, ::MOI.ObjectiveValue)
    if model.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    return model.inner.best_sol
end

function MOI.get(model::Optimizer, ::MOI.VariablePrimal, vi::MOI.VariableIndex)
    if model.status == MOI.OPTIMIZE_NOT_CALLED
        @error "optimize! not called"
    end
    check_inbounds(model, vi)
    return CS.value(model.inner.search_space[vi.value])
end

Please feel free to ask any question about further implementation and check the code changes on this PR #51.

Now I'll tell a bit about specific problems I had to implement stuff for the constraint solver.

Specific problems

We have to check that each variable is bounded and is either integer or boolean. This is a special case as normally unbounded and real values are the easier way but not for the constraint solver. This is not a standard thing so before we optimize the model I run a check on each variable and if this constraint is not satisfied an error is thrown.

Additionally there are two constraints which are normally not supported by other solvers namely the all different constraint and !=.

The all different constraint is not that complicated as there is a way to create constraints of the form: @constraint(m, [x,y,z] in Set) for that we just need to define a set which needs a dimension.

struct AllDifferentSet <: MOI.AbstractVectorSet
    dimension :: Int64
end

and then convert that into our BasicConstraint struct in:

MOI.supports_constraint(::Optimizer, ::Type{MOI.VectorOfVariables}, ::Type{AllDifferentSet}) = true
   
function MOI.add_constraint(model::Optimizer, vars::MOI.VectorOfVariables, set::AllDifferentSet)
    ...
end

I expected it to be harder which was the main reason why I invented my own user interface :/ Well learned something again and hope you too ;)

Okay what about the != constraint. Of course one can write [x,y] in CS.AllDifferentSet(2) instead of x != y but that is extremely unsatisfying.

Fortunately there is a way in JuMP to also define any kind of symbol constraint in the form vars symbol vars/constant using sensetoset.

struct NotEqualSet{T} <: MOI.AbstractScalarSet 
    value :: T
end

sense_to_set(::Function, ::Val{:!=}) = NotEqualSet(0.0)

MOIU.shift_constant(set::NotEqualSet, value) = NotEqualSet(set.value + value)

and then we need to support x - b != 0 so a single affine function:

MOI.supports_constraint(::Optimizer, ::Type{SAF}, ::Type{NotEqualSet{Float64}}) = true

In the corresponding constrain function we need to check that there are no constraints of the form a - b + c != 0 or a + b != 0 or some kind of other form which we currently don't support.

After I managed to implement that I realized that the objective for the graph coloring problem can't be expressed as well. What we did before is to minimize over all variables. (I use we all the time and here it sounds like I'm less stupid because "we" didn't realize it's a too specialized idea. Sorry...)

Mathieu Besançon pointed out to me that smart people normally create an extra variable for that and then add the constraint @constraint(m, extra_var .>= all_vars) and then simply minimize that one extra variable. If you have some questions about how to implement JuMP interfaces or anything related to JuMP you should check out the Gitter developer chat. I mean you can ask me as well but those people can actually help you ;)

I'm sure I missed something to write about so feel free to comment if you have questions regarding my implementation.

There it is: Implementing the solver options isn't that hard as well: options.jl which includes the combine_options function seen above. Basically another struct which holds the options and then load defaults and replace them if the user wants to change something.

Pruning function

Currently I have a lot of code like:

for constraint in com.constraints
    com.info.pre_backtrack_calls += 1
    feasible = constraint.fct(com, constraint)
    if !feasible
        feasible = false
        break
    end
end

but I already have a prune! function so I should use it. I've added some parameters to it:

function prune!(com::CS.CoM; pre_backtrack=false, all=false, only_once=false, initial_check=false)

to have control whether all constraints should be called or only on the changes and whether a constraint should be called even though all variables of the constraint are fixed already. This needs to be done to check that the start variables actually satisfy the constraints.

Additionally it was clear to me that the order of the constraint is somewhat relevant and I didn't want it to be a factor as I think it's unexpected behavior if it's slower if you change the order of constraints. JuMP orders them in their own fashion anyway so I wanted to be able to compare it more to the current master. That is the reason why order the constraints by the number of open values left (which I did before in the prune! function which I didn't use everywhere...) but now I also order by a hash of each constraint if the number of open values is the same.

Hope you learned how to use JuMP for your own project and maybe a bit about how to integrate your own special constraints into it which might not seem supported by now.

A new post is out: Sum constraint with real coefficients

What's next?

Some ideas but feel free to add yours ;)

If you enjoy the blog in general and at least some of the constraint solver posts I would be very flattered if you consider supporting this blog. There are several ways you can do that:

Thanks for reading and special thanks to my four patrons!

List of patrons paying more than 2$ per month:

Currently I get 9$ per Month using Patreon and PayPal when I reach my next goal of 50$ per month I'll create a faster blog which will also have a better layout:

If you have ideas => Please contact me!

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 and receiving a mail whenever a new post comes out.

Powered by Buttondown.


Subscribe to RSS
© Ole Kröger. Last modified: June 09, 2020. Website built with Franklin.jl.