Author Archives: Julia on μβ

Static lists in Julia

Re-posted from: https://matbesancon.github.io/post/2019-03-30-static-list/

This post explores the possibility to build static lists in Julia, meaning
lists for which the size is known at compile-time. This is inspired by
a post
on a Scala equivalent but will take different roads to see more than a plain port.
Of course, this implementation is not that handy nor efficient but
is mostly meant to push the limits of the type system,
especially a trick of using recursive types as values
(replacing a dependent type system).
Some other references:

The list operations are inspired by the implementation in DataStructures.jl
StaticArrays.jl is a good inspiration for static data structures in Julia

First thoughts: value type parameter

Julia allows developers to define type parameters.
In the case of a list, the most obvious one may be the
type of data it contains:

abstract type MyList{T} end

Some types are however parametrized on other things, if we look at the
definition of AbstractArray for example:

  AbstractArray{T,N}
Supertype for N-dimensional arrays (or array-like types) with elements of type T.

The two type parameters are another type T and integer N for the
dimensionality (tensor rank). The only constraint for a value to be
an acceptable type parameter is to be composed of plain bits, complying
with isbitstype.

This looks great, we could define our StaticList
directly using integers.

"""
A static list of type `T` and length `L`
"""
abstract type StaticList{T,L} end

struct Nil{T} <: StaticList{T,0} end

StaticList{T}() where T = Nil{T}()
StaticList(v::T) where T = Cons(v, Nil{T}())

struct Cons{T,L} <: StaticList{T,L+1}
    h::T
    t::StaticList{T,L}
    function Cons(v::T, t::StaticList{T,L}) where {T,L}
        new{T,L}(v,t)
    end
end

# Usage:
# Cons(3, Nil{Int}()) is of type StaticList{Int,1}
# Cons(4, Cons(3, Nil{Int}())) is of type StaticList{Int,2}

If you try to evaluate this code, you will get an error:

ERROR: MethodError: no method matching +(::TypeVar, ::Int64)

Pretty explicit, you cannot perform any computation on values used as type
parameters. With more complex operations, this could make the compiler hang,
crash or at least perform poorly (we would be forcing the compiler to execute
this code at compile-time).

One way there might be around this is macros or replacing sub-typing with
another mechanism. For the macro-based approach,
ComputedFieldTypes.jl
does exactly that. More discussion on computed type parameters in
[1] and [2].

Edit: using integer type parameters can be achieved using ComputedFieldTypes.jl as such:

julia> using ComputedFieldTypes

julia> abstract type StaticList{T,L} end

julia> struct Nil{T} <: StaticList{T,0} end

julia> @computed struct Cons{T,L} <: StaticList{T,L}
           h::T
           t::StaticList{T,L-1}
           function Cons(v::T, t::StaticList{T,L0}) where {T,L0}
               L = L0+1
               new{T,L}(v,t)
           end
       end

julia> Cons(3, Nil{Int}())
Cons{Int64,1,0}(3, Nil{Int64}())

julia> Cons(4, Cons(3, Nil{Int}()))
Cons{Int64,2,1}(4, Cons{Int64,1,0}(3, Nil{Int64}()))

This might be the neatest option for building the StaticList.

Recursive natural numbers

We can use the same technique as in the Scala post, representing natural
number using recursive types.

ZeroLength is a special singleton type
Next{L} represents the number following the one represented by L

We can modify our previous example:

"""
A type parameter for List length, the numerical length can be retrieved
using `length(l::Length)`
"""
abstract type Length end
struct ZeroLength <: Length end
struct Next{L<:Length} <: Length end

"""
A linked list of size known at compile-time
"""
abstract type StaticList{T,L<:Length} end

struct Nil{T} <: StaticList{T,ZeroLength} end

StaticList{T}() where T = Nil{T}()
StaticList(v::T) where T = Cons(v, Nil{T}())

struct Cons{T,L<:Length} <: StaticList{T,Next{L}}
    h::T
    t::StaticList{T,L}
    function Cons(v::T, t::StaticList{T,L}) where {T,L<:Length}
        new{T,L}(v,t)
    end
end

"""
By default, the type of the Nil is ignored if different
from the type of first value
"""
Cons(v::T,::Type{Nil{T1}}) where {T,T1} = Cons(v, Nil{T}())

We can then define basic information for a list, its length:

Base.length(::Type{ZeroLength}) = 0
Base.length(::Type{Next{L}}) where {L} = 1 + length(L)

Base.eltype(::StaticList{T,L}) where {T,L} = T
Base.length(l::StaticList{T,L}) where {T,L} = length(L)

One thing should catch your attention in this block,
we use a recursive definition of length for the Length type,
which means we can blow our compiler. However, both of the definitions
are static, in the sense that they don’t use type information, so
the final call should reduce to spitting out the length cached at compile-time.
You can confirm this is the case by checking the produced assembly instructions with @code_native.
We respected our contract of a list with size known at compile-time.

Implementing a list-y behaviour

This part is heavily inspired by the DataStructures.jl list implementation,
as such we will not re-define methods with semantically similar but
implement them for our list type. Doing so for your own package
allows user to switch implementation for the same generic code.

The first operation is being able to join a head with an existing list:

DataStructures.cons(v::T,l::StaticList{T,L}) where {T,L} = Cons(v,l)

"""
Allows for `cons(v,Nil)`. Note that the `Nil` type is ignored.
"""
DataStructures.cons(v::T,::Type{Nil}) where {T} = StaticList(v)

(::Colon)(v::T,l::StaticList{T,L}) where {T,L} = DataStructures.cons(v, l)
(::Colon)(v::T,::Type{Nil}) where {T,L} = DataStructures.cons(v, Nil{T})

Implementing the odd ::Colon methods allows for a very neat syntax:

l0 = StaticList{Int}()
l1 = 1:l0
l2 = 2:l1

Unlike the Scala post, we are not using the :: operator which
is reserved for typing expressions in Julia.
We can add a basic head and tail methods, which allow querying
list elements without touching the inner structure. This
will be useful later on.

DataStructures.head(l::Cons{T,L}) where {T,L} = l.h
DataStructures.tail(l::Cons{T,L}) where {T,L} = l.t

Testing list equality can be done recursively, dispatching on the three
possible cases:

==(l1::StaticList, l2::StaticList) = false

function ==(l1::L1,l2::L2) where {T1,L,T2,L1<:Cons{T1,L},L2<:Cons{T2,L}}
    l1.h == l2.h && l1.t == l2.t
end

"""
Two `Nil` are always considered equal, no matter the type
"""
==(::Nil,::Nil) = true

We can now define basic higher-order functions, such as zip below,
and implement the iteration interface.

function Base.zip(l1::Nil{T1},l2::StaticList{T2,L2}) where {T1,T2,L2}
    Nil{Tuple{T1,T2}}
end

function Base.zip(l1::Cons{T1,L1},l2::Cons{T2,L2}) where {T1,L1,T2,L2}
    v = (l1.h, l2.h)
    Cons(v,zip(l1.t,l2.t))
end

Base.iterate(l::StaticList, ::Nil) = nothing
function Base.iterate(l::StaticList, state::Cons = l)
    (state.h, state.t)
end

Iterating over our lists is fairly straight-forward, and will be more efficient than
the recursive implementations of the higher-order functions, we still kept it for
equality checking, more a matter of keeping a functional style in line with the Scala post.

The case of list reversal is fairly straightforward: iterate and accumulate
the list in a new one.

function Base.reverse(l::StaticList{T,L}) where {T,L}
    l2 = Nil{T}
    for h in l
        l2 = Cons(h, l2)
    end
    l2
end

We define the cat operation between multiple lists.

function Base.cat(l1::StaticList{T,L},l2::StaticList{T,L}) where {T,L}
    l = l2
    for e in reverse(l1)
        l = Cons(e, l)
    end
    l
end

The reverse is necessary to keep the order of the two lists.

Special-valued lists

Now that we have a basic static list implementation, we can spice things up.
StaticList is just an abstract type in our case, not an algebraic data type
as in common functional implementations, meaning we can define other sub-types.

Imagine a numeric list, with a series of zeros or ones somewhere.
Instead of storing all of them, we can find a smart way of representing them.
Let us define a static list of ones:

struct OnesStaticList{T<:Number,L<:Length} end

Base.iterate(l::OnesStaticList, ::Type{ZeroLength}) = nothing
function Base.iterate(l::OnesStaticList{T,L}, state::Type{Next{L1}} = L) where $
    (one(T), L1)
end

This list corresponds to the 1 value of type T, repeated for all elements.
In a similar fashion, one can define a ZeroList:

struct ZerosStaticList{T<:Number,L<:Length} end

Base.iterate(l::ZerosStaticList, ::Type{ZeroLength}) = nothing
function Base.iterate(l::ZerosStaticList{T,L}, state::Type{Next{L1}} = L) where$
    (zero(T), L1)
end

One thing to note is that these lists are terminal, in the sense that they cannot
be part of a greater list. To fix this, we can add a tail to these as follows:

struct ZerosStaticList{T<:Number,L<:Length,TL<:StaticList{T,<:Length}}
	t::TL
end

Base.iterate(l::ZerosStaticList, ::Type{ZeroLength}) = l.t
function Base.iterate(l::ZerosStaticList{T,L}, state::Type{Next{L1}} = L) where$
    (zero(T), L1)
end

The t field of the list contains the tail after the series of zeros,
we can thus build a much simpler representation in case of long constant series.
In a similar fashion, one could define a constant list of N elements, storing
the value just once.

Multi-typed lists

There is one last extension we can think of with this data structure.
Since we have a recursive length parameter, why not add it a type at each new node?

abstract type TLength end
struct TZeroLength <: TLength end
struct TNext{T,L<:TLength} <: TLength end

abstract type TStaticList{L<:TLength} end

struct TNil <: TStaticList{TZeroLength} end

struct TCons{T, L<:TLength} <: TStaticList{TNext{T,L}}
    h::T
    t::TStaticList{L}
    function TCons(v::T, t::TStaticList{L}) where {T,L<:TLength}
        new{T,L}(v,t)
    end
end

With such construct, all nodes can be of a different type T, without
removing the type information from the compiler.

julia> TCons(3,TNil())
TCons{Int64,TZeroLength}(3, TNil())

julia> TCons("ha", TCons(3,TNil()))
TCons{String,TNext{Int64,TZeroLength}}("ha", TCons{Int64,TZeroLength}(3, TNil()))

One interesting thing to note here is that the type takes the same
structure as the list itself:

Type: either a T and a TLength containing the rest of the type, or TNil
Data: either a value of a given type and the rest of the list, or empty list

Conclusion

The Julia type system and compiler allow for sophisticated specifications
when designing data structures, which gives it a feel of compiled languages.
This however should not be abused, in our little toy example, the type parameter
grows in complexity as the list does, which means the compiler has to carry out
some computation.

If you want some further compile-time tricks, Andy Ferris’s
workshop at JuliaCon 2018 details how to perform compile-time computations
between bits and then bytes.

If you have any idea how to implement StaticList using integer parameters instead
of custom struct I would be glad to exchange. Porting this to
use ComputedFieldTypes.jl might be a fun
experiment.

Feel free to reach out any way you prefer, Twitter,
email to exchange or discuss this post.

Sources

Header image source: https://pxhere.com/en/photo/742575

[1] A proposal on Julia “Defer calculation of field types until type parameters are known”, julia/issues/18466
[2] Discussion on compile-time computations on Discourse

Static lists in Julia

By: Julia on μβ

Re-posted from: https://matbesancon.xyz/post/2019-03-30-static-list/

The list operations are inspired by the implementation in DataStructures.jl
StaticArrays.jl is a good inspiration for static data structures in Julia

First thoughts: value type parameter

Julia allows developers to define type parameters.
In the case of a list, the most obvious one may be the
type of data it contains:

abstract type MyList{T} end

Some types are however parametrized on other things, if we look at the
definition of AbstractArray for example:

  AbstractArray{T,N}
Supertype for N-dimensional arrays (or array-like types) with elements of type T.

This looks great, we could define our StaticList
directly using integers.

"""
A static list of type `T` and length `L`
"""
abstract type StaticList{T,L} end

struct Nil{T} <: StaticList{T,0} end

StaticList{T}() where T = Nil{T}()
StaticList(v::T) where T = Cons(v, Nil{T}())

struct Cons{T,L} <: StaticList{T,L+1}
    h::T
    t::StaticList{T,L}
    function Cons(v::T, t::StaticList{T,L}) where {T,L}
        new{T,L}(v,t)
    end
end

# Usage:
# Cons(3, Nil{Int}()) is of type StaticList{Int,1}
# Cons(4, Cons(3, Nil{Int}())) is of type StaticList{Int,2}

If you try to evaluate this code, you will get an error:

ERROR: MethodError: no method matching +(::TypeVar, ::Int64)

Edit: using integer type parameters can be achieved using ComputedFieldTypes.jl as such:

julia> using ComputedFieldTypes

julia> abstract type StaticList{T,L} end

julia> struct Nil{T} <: StaticList{T,0} end

julia> @computed struct Cons{T,L} <: StaticList{T,L}
           h::T
           t::StaticList{T,L-1}
           function Cons(v::T, t::StaticList{T,L0}) where {T,L0}
               L = L0+1
               new{T,L}(v,t)
           end
       end

julia> Cons(3, Nil{Int}())
Cons{Int64,1,0}(3, Nil{Int64}())

julia> Cons(4, Cons(3, Nil{Int}()))
Cons{Int64,2,1}(4, Cons{Int64,1,0}(3, Nil{Int64}()))

This might be the neatest option for building the StaticList.

Recursive natural numbers

We can use the same technique as in the Scala post, representing natural
number using recursive types.

ZeroLength is a special singleton type
Next{L} represents the number following the one represented by L

We can modify our previous example:

"""
A type parameter for List length, the numerical length can be retrieved
using `length(l::Length)`
"""
abstract type Length end
struct ZeroLength <: Length end
struct Next{L<:Length} <: Length end

"""
A linked list of size known at compile-time
"""
abstract type StaticList{T,L<:Length} end

struct Nil{T} <: StaticList{T,ZeroLength} end

StaticList{T}() where T = Nil{T}()
StaticList(v::T) where T = Cons(v, Nil{T}())

struct Cons{T,L<:Length} <: StaticList{T,Next{L}}
    h::T
    t::StaticList{T,L}
    function Cons(v::T, t::StaticList{T,L}) where {T,L<:Length}
        new{T,L}(v,t)
    end
end

"""
By default, the type of the Nil is ignored if different
from the type of first value
"""
Cons(v::T,::Type{Nil{T1}}) where {T,T1} = Cons(v, Nil{T}())

We can then define basic information for a list, its length:

Base.length(::Type{ZeroLength}) = 0
Base.length(::Type{Next{L}}) where {L} = 1 + length(L)

Base.eltype(::StaticList{T,L}) where {T,L} = T
Base.length(l::StaticList{T,L}) where {T,L} = length(L)

Implementing a list-y behaviour

The first operation is being able to join a head with an existing list:

DataStructures.cons(v::T,l::StaticList{T,L}) where {T,L} = Cons(v,l)

"""
Allows for `cons(v,Nil)`. Note that the `Nil` type is ignored.
"""
DataStructures.cons(v::T,::Type{Nil}) where {T} = StaticList(v)

(::Colon)(v::T,l::StaticList{T,L}) where {T,L} = DataStructures.cons(v, l)
(::Colon)(v::T,::Type{Nil}) where {T,L} = DataStructures.cons(v, Nil{T})

Implementing the odd ::Colon methods allows for a very neat syntax:

l0 = StaticList{Int}()
l1 = 1:l0
l2 = 2:l1

DataStructures.head(l::Cons{T,L}) where {T,L} = l.h
DataStructures.tail(l::Cons{T,L}) where {T,L} = l.t

Testing list equality can be done recursively, dispatching on the three
possible cases:

==(l1::StaticList, l2::StaticList) = false

function ==(l1::L1,l2::L2) where {T1,L,T2,L1<:Cons{T1,L},L2<:Cons{T2,L}}
    l1.h == l2.h && l1.t == l2.t
end

"""
Two `Nil` are always considered equal, no matter the type
"""
==(::Nil,::Nil) = true

We can now define basic higher-order functions, such as zip below,
and implement the iteration interface.

function Base.zip(l1::Nil{T1},l2::StaticList{T2,L2}) where {T1,T2,L2}
    Nil{Tuple{T1,T2}}
end

function Base.zip(l1::Cons{T1,L1},l2::Cons{T2,L2}) where {T1,L1,T2,L2}
    v = (l1.h, l2.h)
    Cons(v,zip(l1.t,l2.t))
end

Base.iterate(l::StaticList, ::Nil) = nothing
function Base.iterate(l::StaticList, state::Cons = l)
    (state.h, state.t)
end

The case of list reversal is fairly straightforward: iterate and accumulate
the list in a new one.

function Base.reverse(l::StaticList{T,L}) where {T,L}
    l2 = Nil{T}
    for h in l
        l2 = Cons(h, l2)
    end
    l2
end

We define the cat operation between multiple lists.

function Base.cat(l1::StaticList{T,L},l2::StaticList{T,L}) where {T,L}
    l = l2
    for e in reverse(l1)
        l = Cons(e, l)
    end
    l
end

The reverse is necessary to keep the order of the two lists.

Special-valued lists

Imagine a numeric list, with a series of zeros or ones somewhere.
Instead of storing all of them, we can find a smart way of representing them.
Let us define a static list of ones:

struct OnesStaticList{T<:Number,L<:Length} end

Base.iterate(l::OnesStaticList, ::Type{ZeroLength}) = nothing
function Base.iterate(l::OnesStaticList{T,L}, state::Type{Next{L1}} = L) where $
    (one(T), L1)
end

This list corresponds to the 1 value of type T, repeated for all elements.
In a similar fashion, one can define a ZeroList:

struct ZerosStaticList{T<:Number,L<:Length} end

Base.iterate(l::ZerosStaticList, ::Type{ZeroLength}) = nothing
function Base.iterate(l::ZerosStaticList{T,L}, state::Type{Next{L1}} = L) where$
    (zero(T), L1)
end

One thing to note is that these lists are terminal, in the sense that they cannot
be part of a greater list. To fix this, we can add a tail to these as follows:

struct ZerosStaticList{T<:Number,L<:Length,TL<:StaticList{T,<:Length}}
	t::TL
end

Base.iterate(l::ZerosStaticList, ::Type{ZeroLength}) = l.t
function Base.iterate(l::ZerosStaticList{T,L}, state::Type{Next{L1}} = L) where$
    (zero(T), L1)
end

Multi-typed lists

There is one last extension we can think of with this data structure.
Since we have a recursive length parameter, why not add it a type at each new node?

abstract type TLength end
struct TZeroLength <: TLength end
struct TNext{T,L<:TLength} <: TLength end

abstract type TStaticList{L<:TLength} end

struct TNil <: TStaticList{TZeroLength} end

struct TCons{T, L<:TLength} <: TStaticList{TNext{T,L}}
    h::T
    t::TStaticList{L}
    function TCons(v::T, t::TStaticList{L}) where {T,L<:TLength}
        new{T,L}(v,t)
    end
end

With such construct, all nodes can be of a different type T, without
removing the type information from the compiler.

julia> TCons(3,TNil())
TCons{Int64,TZeroLength}(3, TNil())

julia> TCons("ha", TCons(3,TNil()))
TCons{String,TNext{Int64,TZeroLength}}("ha", TCons{Int64,TZeroLength}(3, TNil()))

One interesting thing to note here is that the type takes the same
structure as the list itself:

Type: either a T and a TLength containing the rest of the type, or TNil
Data: either a value of a given type and the rest of the list, or empty list

Conclusion

If you want some further compile-time tricks, Andy Ferris’s
workshop at JuliaCon 2018 details how to perform compile-time computations
between bits and then bytes.

If you have any idea how to implement StaticList using integer parameters instead
of custom struct I would be glad to exchange. Porting this to
use ComputedFieldTypes.jl might be a fun
experiment.

Feel free to reach out any way you prefer, Twitter,
email to exchange or discuss this post.

Sources

Header image source: https://pxhere.com/en/photo/742575

[1] A proposal on Julia “Defer calculation of field types until type parameters are known”, julia/issues/18466
[2] Discussion on compile-time computations on Discourse

Multiple dispatch – an example for mathematical optimizers

By: Julia on μβ

Re-posted from: https://matbesancon.github.io/post/2019-02-24-multiple-dispatch-optimizers/

In a recent pull request on a personal project, I spent some time designing
an intuitive API for a specific problem. After reaching a satisfying result,
I realized this would never have been possible without one of the central
mechanisms of the Julia language: multiple dispatch. Feel free to read the
Julia docs on the topic
or what Wikipedia has to say
about it.

This post is a walkthrough for multiple dispatch for a case in mathematical
optimization. The first part will introduce the problem context and requires
some notion in mathematical optimization, if this stuff is scary, feel free to
skip to the rest directly.

Refresher on if-then-else constraints

I promised an example oriented towards mathematical optimization, here it is:
it is common to model constraints with two variables $(x, y)$,
$x$ continuous and $y$ binary stating:

$y = 0 \Rightarrow x = 0$
If $y = 1$, there is no specific constraint on $x$

Some examples of models with such constraint:

Facility location: if a wharehouse is not opened, $y = 0$, then the quantity
served by this point has to be $x = 0$, otherwise, the quantity can go up to
the wharehouse capacity.
Unit commitment (a classic problem for power systems): if a power plant
has not been activated for a given hour, then it cannot supply any power,
otherwise, it can supply up to its capacity.
Complementarity constraints: if a dual variable $\lambda$ is 0,
then the corresponding constraint is not active (in non-degenerate cases,
the slack variable is non-zero)

Logical constraints with such if-then-else structure cannot be handled by
established optimization solvers, at least not in an efficient way. There are
two usual ways to implement this, “big-M” type constraints and special-ordered
sets of type 1 SOS1.

A SOS1 constraint specifies that out of a set of variables or expressions,
at most one of them can be non-zero. In our case, the if-then-else constraint
can be modeled as:
$$SOS1(x,\, 1-y)$$

Most solvers handling integer variables can use these $SOS1$ constraints
within a branch-and-bound procedure.

The other formulation is using an upper-bound on the $x$ variable, usually
written $M$, hence the name:

$$x \leq M \cdot y $$

If $y=0$, $x$ can be at most 0, otherwise it is bounded by $M$. If $M$
is sufficiently big, the constraint becomes inactive.
However, smaller $M$ values yield tighter formulations, solved more efficiently.
See Paul Rubin’s
detailed blog post on the subject. If we want bounds as tight as possible, it
is always preferable to choose one bound per constraint, instead of one unique
$M$ for them all, which means we need a majorant of all individual $M$.

As a rule of thumb, big-M constraints are pretty efficient if $M$ is tight,
but if we have no idea about it, SOS1 constraints may be more interesting,
see [1] for recent numerical experiments applied to bilevel problems.

Modeling if-then-else constraints

Now that the context is set, our task is to model if-then-else constraints
in the best possible way, in a modeling package for instance. We want the user
to specify something as:

function handle_ifthenelse(x, y, method, params)
    # build the constraint with method using params
end

Without a dispatch feature baked within the language, we will end up doing
it ourselves, for instance in:

function handle_ifthenelse(x, y, method, params)
    if typeof(method) == SOS1Method
        # model as SOS1Method
    elseif typeof(method) == BigMMethod
        # handle as big M with params
    else
        throw(MethodError("Method unknown"))
    end
end

NB: if you have to do that in Julia, there is a isa(x, T) function
verifying if x is a T in a more concise way, this is verifying sub-typing
instead of type equality, which is much more flexible.

The function is way longer than necessary, and will have to be modified every
time. In a more idiomatic way, what we can do is:

struct SOS1Method end
struct BigMMethod end

function handle_ifthenelse(x, y, ::SOS1Method)
    # model as SOS1Method
end

function handle_ifthenelse(x, y, ::BigMMethod, params)
    # handle as big M with params
end

Much better here, three things to notice:

This may look similar to pattern matching in function arguments if you are
familiar with languages as Elixir. However, the method to use can be determined
using static dispatch, i.e. at compile-time.
We don’t need to carry around params in the case of the SOS1 method,
since we don’t use them, so we can adapt the method signature to pass only
what is needed.
This code is much easier to document, each method can be documented on
its own type, and the reader can refer to the method directly.

Cherry on top, any user can define their own technique by importing our function
and defining a new behavior:

import OtherPackage # where the function is defined

struct MyNewMethod end

function handle_ifthenelse(x, y, ::MyNewMethod)
    # define a new method for ifthenelse, much more efficient
end

Handling big M in an elegant way

We have seen how to dispatch on the technique, but still we are missing one
point: handling the params in big-M formulations. If you have pairs of $(x_j,y_j)$,
then users may want:

$$ x_j \leq M_j \cdot y_j\,\, \forall j $$

Or:
$$ x_j \leq M \cdot y_j\,\, \forall j $$

The first formulation requires a vector of M values, and the second one
requires a scalar. One default option would be to adapt to the most general one:
if several M values are given, build a vector, if there is only one, repeat it
for each $j$. One way to do it using dynamic typing:

struct BigMMethod end

function handle_ifthenelse(x, y, ::BigMMethod, M::Union{Real,AbstractVector{<:Real}})
    if M isa Real
        # handle with one unique M
    else
        # it is a vector
        # handle with each M[j]
    end
end

Note that we can constrain the type of M to be either a scalar or a Vector
using Union type. Still, this type verification can be done using dispatch,
and we can handle the multiple cases:

struct BigMMethod end

"""
Use one unique big M value
"""
function handle_ifthenelse(x, y, ::BigMMethod, M::Real)
    # handle with one unique M
end

"""
Use a vector of big M value
"""
function handle_ifthenelse(x, y, ::BigMMethod, Mvec::AbstractVector)
    # handle with each Mvec[j]
end

This solution is fine, and resolving most things at compile-time.
Also, note that we are defining one signature as a convenience way redirecting
to another.

Polishing our design: enriched types

The last solution is great, we are dispatching on our algorithm and parameter
types. However, in a realistic research or development work, many more
decisions are taken such as algorithms options, number types, various parameters.
We will likely end up with something similar to:

function do_science(x, y, z,
                    ::Alg1, params_alg_1,
                    ::Alg2, params_alg_2,
                    ::Alg3, # algortithm 3 does not need parameters
                    ::Alg4, params_alg_4)
    # do something with params_alg_1 for Alg1
    # do something with params_alg_2 for Alg2
    # ...
end

Requiring users to pass all arguments and types in the correct order.
A long chain of positional arguments like this end makes for error-prone
and cumbersome interfaces. Can we change this? We created all our types as
empty structures struct A end and use it just to dispatch. Instead,
we could store adapted parameters within the corresponding type:

struct Alg1
    coefficient::Float64
    direction::Vector{Float64}
end

# define other types

function do_science(x, y, z, a1::Alg1, a2::Alg2, ::Alg3, a4::Alg4)
    # do something with params_alg_1 for Alg1
    # a1.coefficient, a1.direction...
    # do something with Alg2
    # ...
end

Getting back to our initial use case of BigMMethod, we need to store
the $M$ value(s) in the structure:

struct BigMMethod
    M::Union{Float64, Vector{Float64}}
end

This seems fine, however, the Julia compiler cannot know the type of the M
field at compile-time, instead, we can use a type parameter here:

struct BigMMethod{MT<:Union{Real, AbstractVector{<:Real}}}
    M::MT
    BigMMethod(M::MT) where {MT} = new{MT}(M)
end

When constructing the BigMMethod with this definition, it can be specialized
on MT, the type of M, two examples of valid definitions are:

BigMMethod(3.0)
# result: BigMMethod{Float64}(3.0)

BigMMethod(3)
# result: BigMMethod{Int}(3)

BigMMethod([3.0, 5.0])
# result BigMMethod{Vector{Float64}}([3.0, 5.0])

The advantage is we can now specialize the handle_ifthenelse
signature on the type parameter of M, as below:

"""
Use one unique big M value
"""
function handle_ifthenelse(x, y, bm::BigMMethod{<:Real})
    # handle with one unique M bm.M
end

"""
Use a vector of big M value
"""
function handle_ifthenelse(x, y, bm::BigMMethod{<:AbstractVector})
    # handle with each bm.M[j]
end

The advantage is a strictly identical signature, whatever the method and
its parameters, users will always call it with:
handle_ifthenelse(x, y, bm::BigMMethod{<:AbstractVector})

Conclusion: avoiding a clarity-flexibility trade-off

In this simple but commonly encountered example, we leveraged multiple dispatch,
the ability to choose a function implementation depending on the type of its
arguments. This helped us define a homogeneous interface for specifying a type
of constraint, specializing on the method (SOS1 or big M) and on the data
available (one M or a vector of M values).

Performance bonus, this design is providing the Julia compiler with strong type
information while remaining flexible for the user. In Julia terminology,
this property is called type stability.
We would not have benefitted from this property if we had used reflection-based
design (with typeof() and isa).

This idea of using big-M as an example did not come up in the abstract but is
a simplification of the design used in the
BilevelOptimization.jl
package. Remember I mentioned complementarity constraints, it is exactly this
use case.

If you are interested in more examples of multiple dispatch and hands-on
use cases for the Julia type system, check out
these
two
articles.
Feel free to reach out any way you prefer, Twitter,
email.

Edit 1: thanks BYP for sharp proofreading and constructive critics.

Edit 2: Thanks Mathieu Tanneau for pointing out the alternative solution of
indicator constraints instead of big M, as documented in Gurobi, CPLEX.

Edit 3: For more info on big M constraints and underlying issues, you can read
Thiago Serra’s post, which includes nice visualizations of the problem space.

Sources:

[1] Henrik Carøe Bylling’s thesis, KU, http://web.math.ku.dk/noter/filer/phd19hb.pdf

juliabloggers.com

A Julia Language Blog Aggregator

Author Archives: Julia on μβ

Static lists in Julia

Table of Contents

First thoughts: value type parameter

Recursive natural numbers

Implementing a list-y behaviour

Special-valued lists

Multi-typed lists

Conclusion

Sources

Static lists in Julia

First thoughts: value type parameter

Recursive natural numbers

Implementing a list-y behaviour

Special-valued lists

Multi-typed lists

Conclusion

Sources

Multiple dispatch – an example for mathematical optimizers

Table of Contents

Refresher on if-then-else constraints

Modeling if-then-else constraints

Handling big M in an elegant way

Polishing our design: enriched types

Conclusion: avoiding a clarity-flexibility trade-off