ConstraintSolver.jl Indicator Part 2

A lot of things happened on this blog lately as I've changed the layout completely. If you haven't seen the blog about the blog :D Check out Using Franklin.jl

Now that I don't have two layouts one for Patrons and one for others I can use all the features which only exist in the new layout for everyone and spend more time implementing new ones ;)

Please comment here or reach me otherwise if you have improvements for the new Blog.

Okay enough housekeeping! This blog is again about the constraint solver.

If you haven't read about it before please check out the first post.

I'm basically trying to build a full constraint solver from scratch in Julia, which isn't an easy task. In this series you can read about my journey.

This post is about the indicator constraint again. You might want to check out Indicator Part 1 before continue reading.

This post is about extending JuMP.jl to support constraints of the form:

@constraint(m, b => {[x, y, z] in CS.AllDifferentSet()})

whereas before only

@constraint(m, b => { x + y + z <= 10)})

so affine inner constraints were allowed.

What do we need and how do we know?

Let's just assume we are a user and want to use that constraint:

using JuMP
using ConstraintSolver
const CS = ConstraintSolver

m = Model(CS.Optimizer)
@variable(m, 0 <= x <= 1, Int)
@variable(m, 0 <= y <= 1, Int)
@variable(m, a, Bin)
@constraint(m, a => {[x,y] in CS.AllDifferentSet()})
@objective(m, Max, a)
optimize!(m)

I got this error:

julia> indicator()
ERROR: MethodError: no method matching 
_build_indicator_constraint(::JuMP.var"#_error#68"{Symbol},
::VariableRef, 
::VectorConstraint{VariableRef,ConstraintSolver.AllDifferentSetInternal,VectorShape},
::Type{MathOptInterface.IndicatorSet{MathOptInterface.ACTIVATE_ON_ONE,S}
  where S<:MathOptInterface.AbstractScalarSet})
Closest candidates are:
  _build_indicator_constraint(::Function, ::AbstractVariableRef, ::ScalarConstraint, ::Type{MathOptInterface.IndicatorSet{A,S} where S<:MathOptInterface.AbstractScalarSet}) where A at /home/ole/.julia/dev/JuMP/src/indicator.jl:13

For a constraint which JuMP supports but the solver doesn't the output looks like:

julia> second_order_cone()
ERROR: Constraints of type MathOptInterface.VectorOfVariables-in-MathOptInterface.SecondOrderCone are not supported by the solver and there are no bridges that can reformulate it into supported constraints.

Then one would need to write a MOI.supports_constraint function for that type and actually implement it.

In our case we need to work on the parsing first and change only little on the solver side as we mostly implement the indicator constraint already. In our case there is no real difference of what the inner constraint looks like.

The error message shows us that we need to create a _build_indicator_constraint function and also points us to where that function exists in JuMP for ScalarConstraint.

function _build_indicator_constraint(
    _error::Function, variable::AbstractVariableRef,
    constraint::ScalarConstraint, ::Type{MOI.IndicatorSet{A}}) where A

    set = MOI.IndicatorSet{A}(moi_set(constraint))
    return VectorConstraint([variable, jump_function(constraint)], set)
end

build_indicator_constraint

We need to have

"""
    Support for indicator constraints with a set constraint as the right hand side
"""
function JuMP._build_indicator_constraint(
    _error::Function, variable::JuMP.AbstractVariableRef,
    constraint::JuMP.VectorConstraint, ::Type{MOI.IndicatorSet{A}}) where A

    set = CS.IndicatorSet{A}(variable, MOI.VectorOfVariables(constraint.func), constraint.set, 1+length(constraint.func))
    vov = VariableRef[variable]
    append!(vov, constraint.func)
    return JuMP.VectorConstraint(vov, set)
end

So we need constraint::VectorConstraint here to not only support the all different constraint but all kind of in Set constraints.

We can then extract the variable indices of the inner constraint and create a new set. The new set is needed because of:

struct IndicatorSet{A, S <: AbstractScalarSet} <: AbstractVectorSet
    set::S
    IndicatorSet{A}(set::S) where {A, S <: AbstractScalarSet} = new{A,S}(set)
end

Here S only allows AbstractScalarSet. The set is defined in MOI and you can jump to the code directly when clicking on the link below the code ;) Oh yeah a new feature...

I've just added my own set for this:

struct IndicatorSet{A} <: MOI.AbstractVectorSet
    ind_var :: JuMP.VariableRef
    func    :: MOI.VectorOfVariables
    set     :: MOI.AbstractVectorSet
    dimension::Int
end

I'm not only saving func and set but also ind_var which is the activation variable. Currently I'm not using the variable so I could remove that. The func and set are used for dispatching later and set should be called set because it is the name in the original MOI.IndicatorSet. The dimension which we set to 1+length(constraint.func) is important as JuMP and MOI not work as one might naively think :D I'm still quite sure that I don't know how it works exactly but this is a very rough outline of what I think. Please correct me if I'm wrong in the comments.

How JuMP and MOI (don't) work

One might think (at least I did)

Writing

m = Model(CS.Optimizer)
@variable(m, 0 <= x <= 2, Int)

calls add_variable and add_constraint for lower and upper bound or Interval such that there is a one to one correspondents to my MOI.add_variable and MOI.add_constraint code but this is not how it works.

There are several reasons why this is not the case. One obvious one is that JuMP allows bridges such that when interval is not supported (which actually doesn't get called even if it is supported in this case) it will be reformulated to 0 <= x and x <= 2. Some more of these bridges exist and some are more complicated than this trivial example i.e I don't have to support ACTIVATE_ON_ZERO in the indicator constraint which the user would write like

@constraint(m, !b => { ... })

Because it can be reformulated by adding an extra variable i.e. y such that y is the opposite of b (an extra constraint) and then y => { ... }.

I think what JuMP and MOI roughly do is creating an internal model which gets all variables and constraints and when calling optimize! it gets copied over to the actual optimizer. In this step the order is actually different than the order one would normally think of:

• variables

• constraints

just like the user modeled it...

The so called VectorOfVariables constraints like a in Set() are added first and those variables are added first then the constraint and then the rest of the variables. I'm not too sure why this is the case. I think it has something to do with the fact that some solvers which are not written in Julia support it only in this way. Again please correct me if I'm talking bullshit. ;)

Nevertheless I had the wrong dimension first (somehow I've set it to 0) probably because I didn't know why I need it :/ and that was the reason why variables weren't added to the model because MOI "thought": "Oh this constraint doesn't have variables so I don't need to add them".

Continuation of building the constraint

After adding JuMP._build_indicator_constraint we now have to support that constraint:

function MOI.supports_constraint(
    ::Optimizer,
    func::Type{MOI.VectorOfVariables},
    set::Type{IS},
) where {A, IS<:CS.IndicatorSet{A}}
    return A == MOI.ACTIVATE_ON_ONE || A == MOI.ACTIVATE_ON_ZERO
end

So here I used the same structure as in the other supports_constraint but use CS.IndicatorSet here.

Additionally of course:

function MOI.add_constraint(
    model::Optimizer,
    vars::MOI.VectorOfVariables,
    set::IS,
) where {A, T<:Real,IS<:CS.IndicatorSet{A}}

In it I have a normal:

internals = ConstraintInternals(
    length(com.constraints) + 1,
    vars,
    set,
    Int[v.value for v in vars.variables]
)

as well as for the inner constraint:

inner_internals = ConstraintInternals(
    0,
    MOI.VectorOfVariables(vars.variables[2:end]),
    set.set,
    Int[v.value for v in vars.variables[2:end]]
)
inner_constraint = init_constraint_struct(typeof(set.set), inner_internals)

Here the 2:end is that we don't have the first variable which is the indicator variable be part of the inner constraint. The last line creates the correct struct depending on what inner constraint we have.

i.e I've defined:

function init_constraint_struct(::Type{AllDifferentSetInternal}, internals)
    AllDifferentConstraint(
        internals,
        Int[], # pval_mapping will be filled later
        Int[], # vertex_mapping => later
        Int[], # vertex_mapping_bw => later
        Int[], # di_ei => later
        Int[], # di_ej => later
        MatchingInit(),
        Int[]
    )
end

The last thing we have to do is to update the prune_constraint! and other _constraint functions such that they support CS.IndicatorSet and MOI.IndicatorSet.

That looks like this:

function prune_constraint!(
     com::CS.CoM,
     constraint::IndicatorConstraint,
-    fct::VAF{T},
+    fct::Union{MOI.VectorOfVariables, VAF{T}},
     set::IS;
     logs = true,
-) where {A, T<:Real, ASS<:MOI.AbstractScalarSet, IS<:MOI.IndicatorSet{A, ASS}}
+) where {A, T<:Real, ASS<:MOI.AbstractScalarSet, IS<:Union{IndicatorSet{A}, MOI.IndicatorSet{A, ASS}}}

The rest stayed basically the same. This is available in v0.2.0 as I had some bugs in v0.1.8 :D

Hope you found this helpful and also enjoyed the new layout and functionality like diffs and links to the GitHub code snippets.

Thanks for reading and special thanks to my six patrons!

List of patrons paying more than 4$ per month:

• Site Wang

• Gurvesh Sanghera

• Szymon Bęczkowski

Currently I get 19$ per Month using Patreon and PayPal.

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).

I'll keep you updated on Twitter OpenSourcES as well as my more personal one: Twitter Wikunia_de

Want to be updated? Consider subscribing on Patreon for free