Survival Analysis for MIP Solver Running Time¶

Inspired by the recent activity in open-source MIP solvers (HiGHS is becoming stable after the 1.0 release and many users submitting bug reports, SCIP has become open-source software), I wanted to revisit some folklore about running times.

When solving difficult problems, in particular, when running computational experiments, we have to choose a timelimit to tell the solver when to give up solving the problem (that is finding a solution and proving its optimality). If the timelimit is chosen too small, many instances will not be solved within it and therefore we can not compare the running time between different solvers or settings. But if we set the timelimit to long, we might waste computational resources while waiting for a result that never comes.

I remember the rule of thumb:

If the problem is not solved within an hours, it will also not be solved after four hours (or later).

Let's see if we can use some data and survival analysis to confirm this rule.

Data Prep¶

We use running time data from the Mittelmann benchmarks on the MIPLIB set of benchmark instances, a set of representative difficult MIP instances.

From the solver comparison, we only look at the open-source solvers still in active development: HiGHS and SCIP. The timelimit was set as 2 hours (7200 seconds), which we use as our window of observation. The tables contain either the time to optimality (solved instance) or some exceptional status (either timout or an error). We interpret all exceptions as right-censored data.

Kaplan-Meier Estimate¶

We apply the Kaplan-Meier estimate to draw survival curves for each solver.

These show the probability that an instance has survived for some time (ie, has not yet been solved to optimality), together with a 95% confidence interval.

Note that the curves (one for each solver) start at 100% for time 0 and decrease monotonically. The curves stop at 7200s (the time limit set) but have not reached 0%, because of some instances that have not been solved.

We can see a steep decline in the first few minutes, while the curves flattens to a long tail, as predicted by our rule of thumb.

Conditional Survival¶

We can also look at conditional survival, considering only that part of the curves beyond the first hour. That is, we reframe our questions from

What is the probability that an instance has survived (not been solved) after some time?

to

What is the probabilty that an instance is solved after some time of additional waiting, knowing that it wasn't solver after the first hour?

To that end, we merge the data from both solvers and make some adjustments to the resulting survival curve.

The above chart shows the probability than instance is solved after some time (beyond 1h), given that it wasn't solved already. It tells us that an additional hour of computation would only solve about 10% of the remaining instances.

In a recent cookbook, the Healthspan Solution, the authors propose a new way to categorize food: the food triangle.

First of all, only whole food ingredients are considered. Then, rather than grouping foods into one of the macronutrients (fat, protein and carbohydrates), they are sorted by energy density on the one hand (top / down), and by source (animal on the left, plants on the right) on the other:

The recommendation is to eat on the "right side", in particular near the top, in order to get sufficient micronutrients while keeping energy intake in check.

On the other hand, many traditional recipes choose foods from both sides near the bottom, yielding the highest energy density while promoting fat storage due to oxidative priority.

I like this scheme, but I wondered whether most foods could really be assigned one of these sides easily. So this post is about visualizing foods by their energy density and relative source of energy (fat vs carbohydrates).

Energy¶

Our nutrition data lists 142 different foods from different categories with energy (in kJ per 100g) as well as relative weight of the macronutrients.

Since fat has a much higher energy density, we will replace the weight values by the energy contribution per macronutrient.

Food Triangle¶

Our quantified food triangle will be implemented as a scatter plot. For the vertical axis, we use the total energy density (low energy on the top). For the horizontal axis, we take the difference of energy contribution from carbohydrates and fat. So, foods rich in fat tend to be near the left, while foods rich in carbohydrates are on the right. Foods in the middle could have either both fat and carbohydrates, or neither, being rich in protein instead.

To my surprise, most foods do actually line up on either side of the triangle and few are found in the center.

We can also see that all animal foods are on the left, while most plant foods are on the right.

Finally, it's quite obvious that the left side goes a lot "deeper" than the right. Foods rich in fat have a much higher potential for high energy density!

Plant-based Food¶

Let's go into more detail, looking only at plant foods now, split into more categories.

We can see that the different categories form nice clusters in the this view.

Most items on the left side seem to be either nuts or seeds. Let's zoom in on these.

When we remove nuts and seeds, we can get a better picture of the foods near the top of the triangle.

There are still some outliers of higher energy density:

OK, so there are two vegetables high in fat (avocado and olive) and also two types of dried fruits, which might not be considered proper whole foods, I guess.

The other items are either grains or legumes, but not fruits or vegetables.

Near the top we have mostly vegetables and some fruits. The foods lowest in energy density are either cabbages, leafy greens or watery vegetables.

Looking at all vegetables now, we can see that the starchy vegetables are near the bottom. Should peas actually be with the legumes?

All grains and legumes are higher in density than vegetables. It should be noted here that the numbers apply to cooked, not fresh or dried foods. So maybe oats are only at the top here, because they are typically cooked with more water compared to other grains?

Soy beans are quite the outlier here, storing most of their energy as fat.

Animal-based Foods¶

Let's now turn to animal foods, all located on the left side.

There is one outlier with very high energy content:

It's pork belly, with lots of fat. But actually not much more than some nuts.

Near the top, with lower energy density (and high water content), we can find dairy products and shellfish.

Dairy also has significant carbohydrates, but the scallops and mussels seem to store almost all of their energy as protein.

Leaving out the pork belly, the highest energy foods are cheeses. The fishes, meats and organ are rich in fat as well as protein.

Still, I would have expected liver to be farther down. Also, I'm surprised by the wide spread of fat content among different types of fish:

Overall, I was positively surprised by the accuracy of the food triangle schema when looking at the actual numbers.

Energy (Kernel) Density¶

As a final analysis, let's compare the energy density of plants and animals (by count) in a simpler, one-dimensional chart:

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.

Random Instance Data¶

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

Drawing with Makie¶

We define some drawing utitities to be used below.

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

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.

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

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.

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.

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.

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?

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:

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:

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:

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

Let's explore some more time windows:

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.

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

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

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.

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.

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:

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¶

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¶

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:

Computing the Dow Jones Index¶

The Dow Jones index is computed from the prices of its stocks, as a weighted average, where the weight is itself defined through the price of the stock.

We will compute the average price of each stock over the days. Then we will normalize these values by dividing through the total of the average prices. The normalized weights are then multiplied to the daily prices to get the daily value of the index.

In the previous plot, the Dow Jones index is marked in orange. In the following, it will be our target, which we want to approximate by using a small number of the stocks.

Model for Linear Fit¶

Let us start with an optimization model similar to an ordinary linear regression, but using the l1-norm, which is readily formulated as a linear program. In compact vector form, this reads

\begin{align*} \text{minimize} \quad & \lVert w^T P - I \rVert_1 \\ \text{subject to} \quad & w \ge 0 \end{align*}

where $P$ stands for the prices of the individual st, $I$ for our index (the target) and $w$ for the weights we use in our fund. We only allow non-negative weights for use in the portfolio.

We can use a standard linear programming trick and introduce auxiliary variables for the positive and negative parts inside the absolute values and minimize their sum. This is formulated using JuMP:

We could use all days of the year as input for our model, but for our evaluation we only use the days from the first three quarters. That way, we can see how our fit extrapolates through the remaining quarter.

In the scatter chart above, we compare the weights from the Dow Jones index, with the weights found by the LP model. As we can see (the points lie on the diagonal), we recover the actual weights perfectly, even though we are only using a subset of the data.

Model with Sparsity Constraint¶

The last model gave us an exact fit, but used all available stocks.

We are changing it now, and allow it to use only a given number of stocks. For this end, we will introduce additional binary variables to select stocks to become active in our index fund.

The weight of inactive variables must be forced to 0, which we do with a so-called big-M constraint. The constant is chosen in such a way, that any active stock could single-handedly approximate the target index.

We start by using only a single stock and then go up to subsets of 6 stocks out of the 30.

The model seems to be relatively difficult to solve for SCIP, probably because of the big-M constraints, so we set a gap limit of 5%, rather than solving to proven optimality.

We extract the model solutions and concatenate them into a table for visualization.

As we can see, our solutions stick to the target fund quite closely during the training period, but then diverge from it quickly. Allowing more stocks gives us a better fit.

Computing Returns¶

In the previous model, we tried to fit our index fund to the absolute day-to-day closing prices. But the next model will be based on picking representative that are similar, which is defined through correlation of returns.

So let's start by computing daily returns from the closing prices (for the training period).