By: Julia on μβ

Re-posted from: https://matbesancon.xyz/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

By: Mosè Giordano

Re-posted from: http://giordano.github.io/blog/2019-02-21-julia-do-vs-python-with/

I’ve recently started coding in Python more than I’ve used to, and it’s
interesting for me to have a closer look to its features and compare them with
similar Julia’s ones.

The do-block in the Julia programming language and the with statement in
Python can appear to be similar because of some common uses (most notably,
opening a stream and automatically closing it at the end), but are fundamentally
different. In this post I will review them.

with statement in Python

The with
statement has
been introduced in Python with PEP
343 to simplify the control flow of
try/finally
statements.

A good explanation of this feature is given at
http://effbot.org/zone/python-with-statement.htm. In general, you can use
try/finally to handle unmanaged resources like streams:

def controlled_execution(callback):
    set things up
    try:
        callback(thing)
    finally:
        tear things down

def my_function(thing):
    do something with thing

controlled_execution(my_function)

For an object of a class, it is possible to reorganize the above control flow by
using the with statement:

class controlled_execution:
    def __enter__(self):
        set things up
        return thing
    def __exit__(self, type, value, traceback):
        tear things down

with controlled_execution() as thing:
    do something with thing

The __enter__ method prepares the object for later use and returns thing,
which can be then consumed in the body of with. The __exit__ method will be
called in any case at the end of the body of with and can be used to close the
previously opened resources and handle any exception raised within the with
body, if necessary.

Opening a file with with

A common application of the with statement is the open
function:

with open(filename, mode='w') as f:
    do something with f

Just to practically see how the __enter__ and __exit__ methods can be
defined, we implement a very simple class, MyFile, to replicate the behaviour
of the builtin open function when opening a file:

class MyFile:
    def __init__(self, file):
        self.file = file
    def __enter__(self):
        return self.file
    def __exit__(self, type, value, traceback):
        return self.file.close()
    def read(self):
        return self.file.read()

def MyOpen(name, mode='r'):
    return MyFile(open(name, mode))

with MyOpen('file.txt') as f:
    f.read()

do-block in Julia

The
do-block
in Julia addresses a different problem: write in a different way the passing of
a function as argument to another function.

In a simple way, any function func accepting a function as first argument:

function func(f::Function, x, y, z...)
    # do something with arguments x, y, z..., and call f at some point
end

can be called as

func(x, y, z...) do args
    # put here the body of function f(args) that will
    # be passed as first argument to func(f, x, y, z...)
end

For example, the
filter
enables filtering a collection by removing those elements for which the provided
function returns false. We can select the odd numbers between 1 and 100
divisible by 3 with

julia> filter(x -> isinteger(x/3) && isodd(x), 1:100)
17-element Array{Int64,1}:
  3
  9
 15
 21
 27
 33
 39
 45
 51
 57
 63
 69
 75
 81
 87
 93
 99

Here we passed an anonymous
function
as first argument to filter in order to select the desired elements of the
collection passed as second argument. We can rewrite this by using the
do-block in the following way:

julia> filter(1:100) do x
           isinteger(x/3) && isodd(x)
       end
17-element Array{Int64,1}:
  3
  9
 15
 21
 27
 33
 39
 45
 51
 57
 63
 69
 75
 81
 87
 93
 99

The code in the body of the do-block is used to create an anonymous function,
whose arguments are specified after do, that is passed as first argument to
the function before the do.

Opening a file with do

A common use of the do-block is to open a file and automatically close it when
you’re done, just like the with statement in Python. For example:

open("outfile", "w") do io
    write(io, data)
end

In Julia’s Base code, this is achieved by extending the open
method in the
following way:

function open(f::Function, args...; kwargs...)
    io = open(args...; kwargs...)
    try
        f(io)
    finally
        close(io)
    end
end

This method, which accepts a function as first argument, is what makes it
possible to use the do-block in the example above: after effectively opening
the file, the body of the do-block is executed in the try
statement,
and in the end the finally
clause
takes care of safely closing the file, even in case of errors while executing
the function f.

Summary

The with statement in Python is a syntactic sugar for try/finally. This
construct is bound to classes, this means that the __enter__/__exit__
methods are bound to a specific class and they can be defined only once as they
have a fixed signature, so you can’t have multiple __enter__/__exit__
handlers within the same class, or independently from a class. In order for the
with statement to work, both __enter__ and __exit__ must be defined in any
case, however if you don’t really need the try/finally construct you can
still take advantage of the handy with syntax, for example by making the
__exit__ method dummy.

The do-block in Julia is not bound to any type, the only condition to use it
is to have a function taking a function as first argument. The function may
contain a control flow construct, like try/finally, but this is not
mandatory. You can also have different methods for a function to be used with
the do-block, depending on the types and/or the number of the other arguments.