Creation date: 2020-06-17

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.

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
```

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.

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 => { ... }`

.

⚠ Note

I actually do support it such that no extra variable is needed. For a constraint solver this is not hard to achieve ;)

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

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

⚠ Note

I had to add a weird dependency as otherwise tests didn't run through which means I had to increase the minor version number.

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

© Ole Kröger. Last modified: June 27, 2020. Website built with Franklin.jl.