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

Creation date: 2020-01-10

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? • We need to create a JuMP optimizer • Define which variable types we support • Define which constraint types are supported • Define which objectives are supported • Functions to get the status, values and objective value • Solver/Optimizer options 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 • MOI_wrapper.jl -> the general optimizer stuff • variables.jl -> Things for @variable • constraints.jl -> Things for @constraint • objective.jl -> Things for @objective • results.jl -> Getting back the results of the solver 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 ;) • Hopcroft Karp bipartite matching • reuse memory in bipartite matching • Looking at the search tree for bigger graph coloring • Support for additional constraints • Bigger benchmark set (i.e. bigger sudokus such that solving time gets increased) 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: • Star the project on GitHub • Follow me on Twitter • Spread the blog on twitter and mention me ;) • Or if you want to support it financially please consider a donation via Patreon, GitHub or old school PayPal • If you subscribe until the 10th of February 2020 to spend about 4$ per month: We will have a 1:1 call where you can ask me anything about the blog!

Thanks for reading and special thanks to my four patrons!

List of patrons paying more than 2$per month: • Site Wang • Gurvesh Sanghera 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:

• Code which I refer to on the side

• Code diffs for changes

• Easier navigation with Series like Simplex and Constraint Programming

• And more