Refined Delaunay Triangulations for Flow Problems

I am trying to solve a network design problem for gas pipeline networks (see PipeLayout). In the original problem statement, only the coordinates of the source and sink nodes are given. Pipes can be connected in any way, also using Steiner nodes at any location.

But I also consider a discretized variation of the problem, where the pipes are selected from a given ground structure. Initially, I thought of using a regular (square) grid as a ground structure, but this brings problems related to solution symmetry (among other things).

So, instead I will now try to find a (not so regular) triangulation of the area that starts with the given source & sink nodes, but may add additional points for refinement.

Mesh generation with triangles is a well-studied problems, with mature implementations (see Triangle), but is usually done in order to solve PDE systems or similar. So these tools define a good quality triangulation in terms of triangles with large angles or small areas.

For our application, the triangle cells themselves are not relevant, but rather the edges around them. These are used to connect the nodes, so they should provide good candidates for shortest paths and spanning trees. Therefore, we will evaluate several triangulations based on the minimum tree that spans all terminal nodes.

Implementation notes.

As usual, we will experiment with Julia code. The triangulation code is part of the PipeLayout package, and based on TriangleMesh. Visualizations are made with Makie and Cairo. Optimization models are formulated with JuMP and solved with SCIP.

In [1]:
# Activate custom environment
using Pkg
Pkg.activate("../../environments/refined-delaunay-for-flow-problems/");

# Julia standard library
using LinearAlgebra
using Random
using SparseArrays
using Statistics

# Third-party packages
using AbstractPlotting
using CairoMakie
using JuMP
using LightGraphs
using SCIP

# My own code (not registered)
using PipeLayout
using PipeLayout: delaunay_triangulation, make_pslg, TriangleSwitches, triangulate
using PipeLayout: antiparallel_digraph, euclidean_steiner_tree

Random Instance Data

We generate a single set of 7 nodes that is used throughout the experiment.

In [2]:
Random.seed!(0);
const N = 7
const WIDTH = 500
const HEIGHT = 300

x = round.(0.95 * WIDTH * rand(N))
y = round.(0.95 * HEIGHT * rand(N))
points = [x y];

Drawing with Makie

We define some drawing utitities to be used below.

In [3]:
function make_scene(width=WIDTH, height=HEIGHT)
    return Scene(resolution=(width, height), show_axis=false, scale_plot=false)
end

function draw_points!(points; markersize=4, color=:black)
    scatter!(points[:, 1], points[:, 2], markersize=markersize, color=color)
end

function draw_edges!(points, edges; color=:gray)
    linesegments!(points[edges'[:], :], color=color)
end

function draw_triangulation(points, edges)
    make_scene()
    draw_edges!(points, edges)
    draw_points!(points)   
end

draw_triangulation(trimesh) = draw_triangulation(trimesh.point', trimesh.edge');

Triangulation with TriangleMesh

TriangleMesh wraps the Triangle C library for Julia. They provide, in particular, an algorithm for Constrained Delaunay Triangulation. This takes as input a planar polygon and outputs a triangulation of the nodes, while ensuring that the existing edges are part of the triangulation (but may be subdivided).

In [4]:
# Create Delaunay Triangulation from points
deltri = delaunay_triangulation(points)

# Create Polygon from triangulation (as input for further refinements)
delpoly = make_pslg(deltri)

draw_triangulation(deltri)
Out[4]:
In [5]:
# Create triangulation of the raw points, specifying a minimum angle constraint
draw_triangulation(triangulate(points,
        TriangleSwitches(minimum_angle=20, conforming_delaunay=true)))
Out[5]:
In [6]:
# Create a CDT, keeping the edges of the original triangulation 
draw_triangulation(triangulate(delpoly,
        TriangleSwitches(minimum_angle=20, conforming_delaunay=true)))
Out[6]:

The simple Delaunay Triangulation makes sense visually and even has nice theoretical properties about approximating Euclidean distances with shortest paths!

The second triangle mesh introduces few additional nodes in order to get more regular triangles. But in turn, it also loses some of the "direct" connections between the terminal nodes.

Finally, the CDT keeps the direct edges of the Delaunay Triangulation, but refines the triangles, using many more additional nodes and edges.

Steiner Tree Model

We will now define an optimization model to determine the shortest subgraph that spans the given terminal nodes (Steiner tree in graphs). The model uses a multi-commodity flow formulation.

With the help of this model, we can see which edges of the triangulation are actually useful for our application, as well as quantitatively compare the triangulations.

In [7]:
function edge_lengths(points, edges)
    diff = points[edges[:, 1], :] - points[edges[:, 2], :]
    return dropdims(mapslices(norm, diff, dims=2), dims=2)
end

edge_lengths(trimesh) = edge_lengths(trimesh.point', trimesh.edge')

function edge_length_map(trimesh)
    edges = trimesh.edge'
    lengths = edge_lengths(trimesh)

    length_map = Dict{Tuple{Int64, Int64}, Float64}()
    for e = 1:size(edges, 1)
        s, t = edges[e, :]
        length_map[s, t] = lengths[e]
        length_map[t, s] = lengths[e]
    end
    return length_map
end

function steiner_tree(trimesh, terminals)
    # Compute length by Euclidean distance of nodes.
    lengths = edge_length_map(trimesh)
    
    # Build digraph with all antiparallel arcs, for nonnegative flow.
    graph = antiparallel_digraph(trimesh)
    nodes = collect(1:nv(graph))
    arcs = collect(keys(lengths))
    
    length(terminals) >= 2 || error("Need at least 2 terminals.")
    all(terminals .<= nv(graph)) || error("Terminals out of range.")
    root = terminals[1]
    sinks = terminals[2:end]
   
    demand(v, s) = 1.0*(v == s) - 1.0*(v == root)
    
    # Using arc length for fixed capacity cost and multi-commodity flow.
    model = JuMP.direct_model(SCIP.Optimizer(display_verblevel=0))
    @variable(model, select[a in arcs], Bin, container=SparseAxisArray)
    @variable(model, flow[a in arcs,s in sinks]  0, container=SparseAxisArray)
    @constraint(model, balance[v in nodes, s in sinks],
        sum(flow[(n, v), s] - flow[(v, n), s] for n in neighbors(graph, v))
        == demand(v, s))
    @constraint(model, capacity[a in arcs, s in sinks], flow[a, s] <= select[a])
    @objective(model, Min, sum(lengths[a] * select[a] for a in arcs))
    
    optimize!(model)
    
    # Return objective, as well as selected edges.
    return objective_value(model), value.(select)
end
Out[7]:
steiner_tree (generic function with 1 method)

Another utility will take a triangulation as input, compute the optimal subtree and visualize both.

In [8]:
function draw_tree(trimesh, terminals=1:N)
    points, edges = trimesh.point', trimesh.edge'
    obj, select = steiner_tree(trimesh, terminals)
    @show obj
    selected = [select[(s,t)] + select[(t,s)] for (s,t) in eachrow(edges)]
    active_edges = edges[selected .> 0.0, :]
    
    draw_triangulation(trimesh)
    linesegments!(points[active_edges'[:], :], color=:plum, linewidth=3)
    draw_points!(points[terminals, :], markersize=5, color=:teal)
end
Out[8]:
draw_tree (generic function with 2 methods)
In [9]:
draw_tree(deltri)
obj = 693.4867635250124
Out[9]:
In [10]:
tri = triangulate(points, TriangleSwitches(minimum_angle=20, conforming_delaunay=true))
draw_tree(tri)
obj = 700.5151721904267
Out[10]:
In [11]:
tri = triangulate(points, TriangleSwitches(maximum_area=sqrt(WIDTH*HEIGHT), conforming_delaunay=true))
draw_tree(tri)
obj = 697.572096498523
Out[11]:
In [12]:
tri = triangulate(delpoly, TriangleSwitches(maximum_area=sqrt(WIDTH*HEIGHT), conforming_delaunay=true))
draw_tree(tri)
obj = 693.4867635250124
Out[12]:

The Delaunay Triangulation allows for a solution of length 693.5. Other point cloud triangulations based on minimum angle or maximum area constraints yield poorer solutions, of length 700.5 and 687.5, while using a larger graph.

The CDT recovers the solution of the Delaunay Triangulation, but does not use any of the additional nodes or edges introduced in the subdivision.

Euclidean Steiner Tree

We can "cheat" and solve the EST problem directly (based on GeoSteiner). This will gives a lower bound on the length (what we can hope to achieve), and maybe a polygon to be triangulated.

In [13]:
est = euclidean_steiner_tree(points)

make_scene()
draw_edges!(est...)
draw_points!(est[1][1:N, :])
Out[13]:

We can see that 3 additional Steiner nodes are inserted.

Let's now build a triangulation of this tree and compute the length.

Then, we show a point cloud triangulation, based on the terminals and the Steiner nodes of the EST.

In [14]:
estpoly = make_pslg(est...)
draw_tree(triangulate(estpoly, TriangleSwitches(
            minimum_angle=20, maximum_area=sqrt(WIDTH*HEIGHT), convex_hull=true)))
obj = 664.018124937279
Out[14]:
In [15]:
draw_tree(triangulate(est[1], TriangleSwitches(minimum_angle=20, maximum_area=sqrt(WIDTH*HEIGHT))))
obj = 696.9337210443265
Out[15]:

The EST has a length of 664, which is significantly shorter than the 693.5 of the Delaunay Triangulation, so there is still some room for improvement.

Using only the additional Steiner nodes does not give a convincing solution.

Conclusion

Overall, the experiment is a little disappointing. Using any point cloud triangulation does not result in short trees, even if the triangulation becomes finer.

The Delaunay Triangulation is a strong candidate, given that it does not introduce any new nodes.

Improving on the Delaunay Triangulation with refinement based on CDT has not been fruitful on this example data, as the additional edges are never used for the tree.

Intuitively, to find short trees, one would like to connect the terminal nodes with direct, line-of-sight edges. But the irregular triangulation often runs orthogonal to that desired line, which would require frequent detours.

Finally, the actual solution of the EST can be used as the basis of a triangulation, which then contains the relevant edges and also looks regular and "natural". But from a methodology point of view, this means using the final result as input of the solution process.

Outlook

Luckily, the actual application has additional features, beyond the minimum length spanning tree. So the EST could still be used effectively in a preprocessing step. Maybe some pertubation of the node positions can be applied to yield triangles of more equal size, as in force-directed_graph_drawing.

Sum Up Rounding for Ramp Down of Medication Intake

The paper Solving Mixed-integer Control Problems by Sum Up Rounding With Guaranteed Integer Gap by Sager et al introduces a rounding strategy for control functions over time with respect to accumulated rounding error.

In this post, we want to apply Sum Up Rounding to the timing of medicine intake, in particular to the case where we would like to ramp down the dose, but it is not practical to take fractional amounts. So, rather than taking half a pill daily, we would take a whole pill every other day.

The central question is:

On which days should we take our medicine if we wanted to ramp down the dose linearly over a given period, say, 4 weeks?

In [1]:
using Random

using DataFrames
using VegaLite
In [2]:
"Sum Up Rounding of a sequence of continuous values between 0 and 1."
function sum_up_round(x)
    all(0.0 .<= x .<= 1.0) || error("Input must lie in [0, 1]!")
    
    result = zeros(size(x))
    sum = 0.0
    for (index, value) in enumerate(x)
        sum += value
        rounded = round(sum)
        result[index] = rounded
        sum -= rounded
    end

    @assert -1.0 <= sum <= 1.0
    @assert all(0.0 .<= result .<= 1.0)
    
    return result
end
Out[2]:
sum_up_round

The implementation is straight-forward: We iterate over the sequence of continuous control values, round at each time step, and add the remainder (rounding error) to the running sum.

Let's apply it to some sequences to see it in effect! First, random numbers:

In [3]:
Random.seed!(1)

random = rand(10)
rounded = sum_up_round(random)
@show random rounded

data = vcat(DataFrame(day=1:10, dose=random, mode="continuous"),
            DataFrame(day=1:10, dose=rounded, mode="discrete"))
@vlplot(
    mark={:bar, opacity=0.7},
    height=80, width=500, data=data,
    x="day:o", y={:dose, stack=nothing}, row=:mode
)
random = [0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666]
rounded = [0.0, 1.0, -0.0, -0.0, 0.0, 1.0, 1.0, 1.0, -0.0, 1.0]
Out[3]:
modecontinuous0.00.51.0dosediscrete0.00.51.0dose12345678910day

As we can see in the above plot, we have a value of 1 on day 2 in the result, even though the original control has a value below 0.5, because the method takes into account also the values from previous results (that were rounded down).

Let's consider a wave pattern:

In [4]:
days = 1:80
wave = (cos.(2π/length(days) .* days) .+ 1) ./ 2
rounded = sum_up_round(wave)

data = vcat(DataFrame(day=days, dose=wave, mode="continuous"),
            DataFrame(day=days, dose=rounded, mode="discrete"))
@vlplot(
    mark={:bar, opacity=0.7},
    height=80, width=500, data=data,
    x="day:o", y={:dose, stack=nothing}, row=:mode
)
Out[4]:
modecontinuous0.00.51.0dosediscrete0.00.51.0dose1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980day

The result is not fully symmetric, but we can see a higher frequencies of 1s where the original control is relative large. Not also the value at day 39, where the original control is almost 0, but we have to correct for the rounded down values for the preceding 12 days!

Let's now answer our actual question, about the Sum Up Rounding of a ramped down medicine intake:

In [5]:
"Linear ramp down covering a given number of days."
ramp_down(length) = range(1.0, 0.0, length=length)
Out[5]:
ramp_down
In [6]:
days = 1:80
ramp = ramp_down(length(days))
rounded = sum_up_round(ramp)

data = vcat(DataFrame(day=days, dose=ramp, mode="continuous"),
            DataFrame(day=days, dose=rounded, mode="discrete"))
@vlplot(
    mark={:bar, opacity=0.7},
    height=80, width=500, data=data,
    x="day:o", y={:dose, stack=nothing}, row=:mode
)
Out[6]:
modecontinuous0.00.51.0dosediscrete0.00.51.0dose1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980day

As expected, the frequency of 1 outputs is decreasing, but not monotonically!

Let's explore some more time windows:

In [7]:
parts = DataFrame[]
for w in 1:8
    n = 7*w
    ramp = collect(ramp_down(n))
    push!(parts, DataFrame(weeks=w, day=1:n, mode="continuous", dose=ramp))
    rounded = sum_up_round(ramp)
    push!(parts, DataFrame(weeks=w, day=1:n, mode="binary", dose=rounded))
end
data = vcat(parts...)

@vlplot(
    mark={:bar, opacity=0.7},
    height=40, width=500, data=data,
    x="day:o", y={:dose, stack=nothing}, row="weeks:n", color=:mode
)
Out[7]:
weeks101dose201dose301dose401dose501dose601dose701dose801dose1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556daybinarycontinuousmode

Finally, let's look at the case where we want to ramp down to a lower dose, say 1/3, instead of 0, and then stay there for another week.

In [8]:
parts = DataFrame[]
high, low = 1, 1/3
for w in 3:8
    n = 7*w
    ramp = vcat(range(high, low, length=n - 7), repeat([low], 7)) 
    push!(parts, DataFrame(weeks=w, day=1:n, mode="continuous", dose=ramp))
    rounded = sum_up_round(ramp)
    push!(parts, DataFrame(weeks=w, day=1:n, mode="binary", dose=rounded))
end
data = vcat(parts...)

@vlplot(
    mark={:bar, opacity=0.7},
    height=40, width=500, data=data,
    x="day:o", y={:dose, stack=nothing}, row="weeks:n", color=:mode
)
Out[8]:
weeks301dose401dose501dose601dose701dose801dose1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556daybinarycontinuousmode

Please don't mistake the above experiments for medical advice!

Managing Exceptions with ResultTypes

Managing Exceptions with ResultTypes

In this post, we want to briefly review some examples of exception handling, present the ResultTypes package as an alternative and finally, show how these mechanisms can cooperate nicely.

Exceptions in Julia

Julia supports Exceptions with mechanisms well-known from other languages:

A function can throw an exception when facing invalid input or a situation that is in some otherwise unexpected. The caller can then handle the exception by wrapping the function call in a try/catch block. If not handled, the exception will travel upwards on the call stack all the way to the REPL (or even stop the process) and print a stack trace. Julia offers several built-in exception types, one can also define custom Exception types, or simply call the error function as a shortcut to throw an ErrorException with a given message.

HTTP example: GET

Let's use a simple example function that wraps HTTP.get and checks the URL for the transport protocol to demonstrate usage of exceptions.

In [1]:
using HTTP
In [2]:
"GET given URL, if using HTTPS, else throw exception."
function https_get1(url)
    if !startswith(url, "https://")
        throw(DomainError("$url: only HTTPS requests allowed!"))
    end
    return HTTP.get(url)
end
Out[2]:
https_get1
In [3]:
# Test it with a happy case, printing the result.
https_get1("https://httpbin.org/status/200")
Out[3]:
HTTP.Messages.Response:
"""
HTTP/1.1 200 OK
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:22 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

"""
In [4]:
# Test it with an error case, showing the stack trace and exception. 
https_get1("http://httpbin.org/status/200")
DomainError with http://httpbin.org/status/200: only HTTPS requests allowed!:


Stacktrace:
 [1] https_get1(::String) at ./In[2]:4
 [2] top-level scope at In[4]:1

In our case, we have an idea of how to fix the error causing issue, by changing the URL. Let's reiterate on our https_get function by handling the exception:

In [5]:
"GET with URL, enforce HTTPS via exception handling."
function https_get2(url)
    try
        return https_get1(url)
    catch e
        if isa(e, DomainError)
            return https_get1(replace(url, "http" => "https", count=1))
        else
            rethrow(e)
        end
    end
end
Out[5]:
https_get2
In [6]:
# Try previous error case again, which is now fixed:
https_get2("http://httpbin.org/status/200")
Out[6]:
HTTP.Messages.Response:
"""
HTTP/1.1 200 OK
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:24 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

"""
In [7]:
# Try another error case, where an exception is thrown by the `HTTP` package,
# which we don't handle. Here, we are again shown the complete (because of rethrow)
# and long stack trace, full of HTTP.jl internals.
https_get2("https://httpbin.org/status/404")
HTTP.ExceptionRequest.StatusError(404, HTTP.Messages.Response:
"""
HTTP/1.1 404 Not Found
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:25 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

""")

Stacktrace:
 [1] #request#1 at /home/rs/.julia/packages/HTTP/GN0Te/src/ExceptionRequest.jl:22 [inlined]
 [2] (::getfield(HTTP, Symbol("#kw##request")))(::NamedTuple{(:iofunction,),Tuple{Nothing}}, ::typeof(HTTP.request), ::Type{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}, ::HTTP.URIs.URI, ::HTTP.Messages.Request, ::Array{UInt8,1}) at ./none:0
 [3] (::getfield(Base, Symbol("###48#49#50")){ExponentialBackOff,getfield(HTTP.RetryRequest, Symbol("##2#3")){Bool,HTTP.Messages.Request},typeof(HTTP.request)})(::Base.Iterators.Pairs{Symbol,Nothing,Tuple{Symbol},NamedTuple{(:iofunction,),Tuple{Nothing}}}, ::Function, ::Type, ::Vararg{Any,N} where N) at ./error.jl:231
 [4] ##48#51 at ./none:0 [inlined]
 [5] #request#1 at /home/rs/.julia/packages/HTTP/GN0Te/src/RetryRequest.jl:44 [inlined]
 [6] #request at ./none:0 [inlined]
 [7] #request#1(::VersionNumber, ::String, ::Nothing, ::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Type{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/MessageRequest.jl:47
 [8] request at /home/rs/.julia/packages/HTTP/GN0Te/src/MessageRequest.jl:28 [inlined]
 [9] #request#1(::Int64, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Type{HTTP.RedirectRequest.RedirectLayer{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/RedirectRequest.jl:24
 [10] request(::Type{HTTP.RedirectRequest.RedirectLayer{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/RedirectRequest.jl:21
 [11] #request#5(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:300
 [12] #request#6 at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:300 [inlined]
 [13] request at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:310 [inlined] (repeats 2 times)
 [14] #get#13 at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:382 [inlined]
 [15] get at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:382 [inlined]
 [16] https_get1(::String) at ./In[2]:6
 [17] https_get2(::String) at ./In[5]:4
 [18] top-level scope at In[7]:1

ResultTypes

There is an alternative pattern of error handling that does not use exceptions, but communicates about success or failures with structural return values from functions. Rather than simply returning the actual result of the function, that value is wrapped in a parameterized type that could also contain an error description, for example in the form of an exception value.

The package ResultTypes provides a Julia implementation of this pattern and the README shows a usage example and benchmarks with performance benefits.

Another possible advantage is that it nudges developers to be more explicit about errors that could happen with function calls and might lead to dealing with the errors in a more local context. Further, the results (including errors) are just values (nothing special/magical) and can be dealt with programmatically. These points are also discussed in a blog post about error handling in Golang.

GET example with ResultTypes

Let's try to see what our example GET wrapper would look like if we used ResultTypes rather than throwing an exception.

In [8]:
using ResultTypes
In [9]:
"GET given URL, if using HTTPS, else return error value."
function https_get3(url)::Result{HTTP.Messages.Response, DomainError}
    if !startswith(url, "https://")
        return DomainError(url, "Insecure protocol!")
    end
    return HTTP.get(url)
end
Out[9]:
https_get3
In [10]:
# Try happy case again, this time with wrapped result:
https_get3("https://httpbin.org/status/200")
Out[10]:
Result(HTTP.Messages.Response:
"""
HTTP/1.1 200 OK
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:25 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

""")
In [11]:
# Try bad case again, resulting in error value.
https_get3("http://httpbin.org/status/200")
Out[11]:
ErrorResult(HTTP.Messages.Response, DomainError("http://httpbin.org/status/200", "Insecure protocol!"))
In [12]:
# Try third-party exception, which is still thrown :-\
https_get3("https://httpbin.org/status/404")
HTTP.ExceptionRequest.StatusError(404, HTTP.Messages.Response:
"""
HTTP/1.1 404 Not Found
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:26 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

""")

Stacktrace:
 [1] #request#1 at /home/rs/.julia/packages/HTTP/GN0Te/src/ExceptionRequest.jl:22 [inlined]
 [2] (::getfield(HTTP, Symbol("#kw##request")))(::NamedTuple{(:iofunction,),Tuple{Nothing}}, ::typeof(HTTP.request), ::Type{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}, ::HTTP.URIs.URI, ::HTTP.Messages.Request, ::Array{UInt8,1}) at ./none:0
 [3] (::getfield(Base, Symbol("###48#49#50")){ExponentialBackOff,getfield(HTTP.RetryRequest, Symbol("##2#3")){Bool,HTTP.Messages.Request},typeof(HTTP.request)})(::Base.Iterators.Pairs{Symbol,Nothing,Tuple{Symbol},NamedTuple{(:iofunction,),Tuple{Nothing}}}, ::Function, ::Type, ::Vararg{Any,N} where N) at ./error.jl:231
 [4] ##48#51 at ./none:0 [inlined]
 [5] #request#1 at /home/rs/.julia/packages/HTTP/GN0Te/src/RetryRequest.jl:44 [inlined]
 [6] #request at ./none:0 [inlined]
 [7] #request#1(::VersionNumber, ::String, ::Nothing, ::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Type{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/MessageRequest.jl:47
 [8] request at /home/rs/.julia/packages/HTTP/GN0Te/src/MessageRequest.jl:28 [inlined]
 [9] #request#1(::Int64, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::Type{HTTP.RedirectRequest.RedirectLayer{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/RedirectRequest.jl:24
 [10] request(::Type{HTTP.RedirectRequest.RedirectLayer{HTTP.MessageRequest.MessageLayer{HTTP.RetryRequest.RetryLayer{HTTP.ExceptionRequest.ExceptionLayer{HTTP.ConnectionRequest.ConnectionPoolLayer{HTTP.StreamRequest.StreamLayer}}}}}}, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/RedirectRequest.jl:21
 [11] #request#5(::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::Function, ::String, ::HTTP.URIs.URI, ::Array{Pair{SubString{String},SubString{String}},1}, ::Array{UInt8,1}) at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:300
 [12] request at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:300 [inlined]
 [13] #request#6 at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:314 [inlined]
 [14] request at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:310 [inlined] (repeats 2 times)
 [15] #get#13 at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:382 [inlined]
 [16] get at /home/rs/.julia/packages/HTTP/GN0Te/src/HTTP.jl:382 [inlined]
 [17] https_get3(::String) at ./In[9]:6
 [18] top-level scope at In[12]:1

As we have seen in the last example, even if we decide to use ResultTypes in our code, we are not safe from exceptions being thrown in the calls below.

Does that limit the scope of the pattern and its implementation with ResultTypes? Can we only use it internally within our libraries, but still deal with exceptions bubbling up from other code?

What should we do about return values in user-facing functions in our library API? We can not expect everybody to learn about and deal with two different patterns of error handling just to use our library.

Systematic Error-Handling in C++

I was reminded of ResultTypes recently, when I researched error handling patterns in C++, and stumbled upon a talk titled "Systematic Error-Handling in C++" by Andrei Alexandrescu (video, slides).

In the first half of the presentation, he motivates and sketches the implementation of an Expected<T> type which is equivalent to what is done in ResultTypes. In addition to the advantages shown above, he also mentions dealing with errors across threads (multiple simultaneous exceptions being thrown), but I'm not sure how this applies to Julia.

Most interesting to me was slide 27 (Icing) with this definition of fromCode:

template <class F>
static Expected fromCode(F fun) {
    try {
        returnExpected(fun());
    } catch(...) {
        return fromException();
    }
}

auto r = Expected<string>::fromCode([] {...});

This provides a bridge between the worlds of exception-throwing and result-returning. Thrown exceptions are captured and instead returned as error values. This applies in particular to the case where we call functions from third-party libraries.

Wrapping Callees

We repeat our example function, but attempt to capture all exceptions (including those from HTTP) and work with result values exclusively.

In [13]:
"GET given URL using HTTPS or return error value."
function https_get4(url)
    if !startswith(url, "https://")
        return ErrorResult(DomainError(url, "Insecure protocol!"))
    end
    try
        # Happy path: wrap result value.
        return Result(HTTP.get(url))
    catch e
        # Turn all exceptions into error values.
        return ErrorResult(e)
    end
end
Out[13]:
https_get4
In [14]:
# Test our own error case.
https_get4("http://httpbin.org/status/200")
Out[14]:
ErrorResult(Any, DomainError("http://httpbin.org/status/200", "Insecure protocol!"))
In [15]:
# Test error case in called function, also resulting in error value,
# no longer an exception with stack trace etc.
https_get4("https://httpbin.org/status/400")
Out[15]:
ErrorResult(Any, HTTP.ExceptionRequest.StatusError(400, HTTP.Messages.Response:
"""
HTTP/1.1 400 Bad Request
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:27 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

"""))

Wrapping for Callers

Let us assume that https_get4 was our library-internal utility function, but we would like to expose a version to our library users. They are expecting exceptions, so we throw them in case of error values:

In [16]:
"GET from URL with HTTPS or throw exception."
function https_get5(url)
    result = https_get4(url)
    if ResultTypes.iserror(result)
        # Error case, get exception and throw it.
        throw(unwrap_error(result))
    end
    # Happy case, return the raw result.
    return unwrap(result)
end
Out[16]:
https_get5
In [17]:
# Test working case, returning raw result.
https_get5("https://httpbin.org/status/200")
Out[17]:
HTTP.Messages.Response:
"""
HTTP/1.1 200 OK
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:28 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

"""
In [18]:
# Test with our error (using HTTP), now thrown as exception.
https_get5("http://httpbin.org/status/200")
DomainError with http://httpbin.org/status/200:
Insecure protocol!

Stacktrace:
 [1] https_get5(::String) at ./In[16]:6
 [2] top-level scope at In[18]:1
In [19]:
# Test with third-party error, also thrown as exception.
https_get5("https://httpbin.org/status/500")
HTTP.ExceptionRequest.StatusError(500, HTTP.Messages.Response:
"""
HTTP/1.1 500 Internal Server Error
Connection: keep-alive
Server: gunicorn/19.9.0
Date: Sat, 26 Jan 2019 13:04:39 GMT
Content-Type: text/html; charset=utf-8
Access-Control-Allow-Origin: *
Access-Control-Allow-Credentials: true
Content-Length: 0
Via: 1.1 vegur

""")

Stacktrace:
 [1] https_get5(::String) at ./In[16]:6
 [2] top-level scope at In[19]:1

Notice that the stack trace in the last example is shallow. That is, it starts from the throw statement in our own functions and does no longer contain the levels below from where the exception originates.

This can be seen as positive or negative, but in any case, I don't know how it could be changed. Use of rethrow is not allowed here, because there is no try/catch block.

Conclusion

We have seen how we can bridge between the patterns of throwing exceptions and return error values, easily and in both directions. Maybe this will convince some developers to use ResultValues in their own packages?

So far, I have only detected its use in Dispatcher.jl. Please share your experiences with error handling in Julia!

Index Funds with Mixed-Integer-Programming

Index Funds with Mixed-Integer-Programming

We will analyze daily price data for stocks in the Dow Jones index and then try to build an accurate index fund using a small numbers of stocks therein.

Similar material was already used in a presentation at PyData Berlin 2017. See the "Tour of popular packages" notebook. Back then, we worked with Julia 0.6 and used the packages DataFrames, Plots and JuMP. Now, we work with Julia 1.0 and use packages from the Queryverse, VegaLite and IndexedTables for data prep and visualization. Also, I added an alternative model.

Loading the Data

In [1]:
using Queryverse
using IndexedTables: ndsparse

The data consists of daily closing prices from 2016 of the 30 stocks that participate in the Dow Jones and are given in a CSV file:

In [2]:
;head -n4 ../files/dowjones2016.csv
Date,Symbol,Price
2016-01-04,AAPL,105.349997999999999
2016-01-04,AXP,67.589995999999999
2016-01-04,BA,140.500000000000000
In [3]:
price = load("../files/dowjones2016.csv") |> ndsparse
price |> @take(3)
Out[3]:
Date Symbol Price
2016-01-04 "AAPL" 105.35
2016-01-04 "AXP" 67.59
2016-01-04 "BA" 140.5
In [4]:
price |> @vlplot(:line, x=:Date, y=:Price, color={"Symbol:o", scale={scheme="viridis"}},
                 width=600, height=450)
Out[4]: