Quantiles

Calculating quantiles is a standard way to summarize data. In the example below we have some small random data p(i,j,k) where we want to calculate quantiles at the 0, 25, 50, 75, and 100% level for each (i,j):

----      9 PARAMETER p                 k1          k2          k3          k4          k5          k6          k7          k8 i1.j1      18.003      84.483      55.487      30.813      29.929      23.181      35.633      85.771 i1.j2       7.644      50.521      99.814      58.295      99.122      76.463      13.939      64.332 i1.j3      16.792      25.758      67.224      44.100      36.610      35.793      14.018      15.860 i2.j1      59.322      83.258      23.851      66.908      77.810      31.062      11.939      50.736 i2.j2      16.857      87.374      27.246      29.296      59.802      72.549      63.197      46.916 i2.j3      41.917      12.652      32.107       5.609      34.516      19.028      64.927      56.514 ----    143 PARAMETER quantiles                 0%         25%         50%         75%        100% i1.j1      18.003      28.242      33.223      62.736      85.771 i1.j2       7.644      41.375      61.313      82.128      99.814 i1.j3      14.018      16.559      30.775      38.483      67.224 i2.j1      11.939      29.259      55.029      69.633      83.258 i2.j2      16.857      28.783      53.359      65.535      87.374 i2.j3       5.609      17.434      33.312      45.566      64.927

In GAMS this is not so easy. Here is one way:

set i /i1*i2/     j /j1*j3/     k /k1*k8/     q /'0%', '25%', '50%', '75%', '100%'/ ; parameter p(i,j,k); p(i,j,k) = uniform(1,100); display p; parameter    v(k) 'vector to be sorted'    r(k) 'sorting rank'    perc0(q) /'0%' EPS, '25%' 25, '50%' 50, '75%' 75, '100%' 100 /    perc(q) 'percentiles'    quantiles(i,j,q) ; $libinclude rank loop((i,j),   v(k) = p(i,j,k);   perc(q) = perc0(q); $ libinclude rank v k r perc   quantiles(i,j,q) = perc(q); ); display quantiles;

We use here the rank utility by Tom Rutherford. A disadvantage is that it is designed to work on vectors. That means we need to loop to extract a k-vector. Under the hood a GDX file is written for each vector, a program called gdxrank is called, and another GDX file is written with the ranks of each vector element. For very large data sets this turns out to be very slow.

Note that in the data statement for perc0 we use EPS instead of 0. This on purpose: EPS will survive the trip to gdxrank, while a pure zero is not stored and not transmitted. If we would have used 0 instead of EPS we would have seen:

----    142 PARAMETER quantiles                25%         50%         75%        100% i1.j1      28.242      33.223      62.736      85.771 i1.j2      41.375      61.313      82.128      99.814 i1.j3      16.559      30.775      38.483      67.224 i2.j1      29.259      55.029      69.633      83.258 i2.j2      28.783      53.359      65.535      87.374 i2.j3      17.434      33.312      45.566      64.927

How would this look in R?

R has the quantile function to calculate percentiles of a vector. We can use the by function to apply this on a slice of the data. This can look like:

> df     i  j  k     value 1  i1 j1 k1 18.002966 2  i1 j1 k2 84.483404 3  i1 j1 k3 55.487160 4  i1 j1 k4 30.812652 5  i1 j1 k5 29.929000 6  i1 j1 k6 23.181234 7  i1 j1 k7 35.633220 8  i1 j1 k8 85.770764 9  i1 j2 k1  7.644259 10 i1 j2 k2 50.520856 11 i1 j2 k3 99.813645 12 i1 j2 k4 58.294604 13 i1 j2 k5 99.122171 14 i1 j2 k6 76.462796 15 i1 j2 k7 13.938556 16 i1 j2 k8 64.332157 17 i1 j3 k1 16.792269 18 i1 j3 k2 25.757973 19 i1 j3 k3 67.223932 20 i1 j3 k4 44.100282 21 i1 j3 k5 36.610326 22 i1 j3 k6 35.792695 23 i1 j3 k7 14.017667 24 i1 j3 k8 15.860077 25 i2 j1 k1 59.322251 26 i2 j1 k2 83.258388 27 i2 j1 k3 23.850758 28 i2 j1 k4 66.907712 29 i2 j1 k5 77.809903 30 i2 j1 k6 31.062189 31 i2 j1 k7 11.938737 32 i2 j1 k8 50.736102 33 i2 j2 k1 16.857103 34 i2 j2 k2 87.373769 35 i2 j2 k3 27.246340 36 i2 j2 k4 29.295618 37 i2 j2 k5 59.801636 38 i2 j2 k6 72.549188 39 i2 j2 k7 63.196619 40 i2 j2 k8 46.915989 41 i2 j3 k1 41.917392 42 i2 j3 k2 12.651840 43 i2 j3 k3 32.107014 44 i2 j3 k4  5.608600 45 i2 j3 k5 34.516477 46 i2 j3 k6 19.027860 47 i2 j3 k7 64.926986 48 i2 j3 k8 56.513809 > q <- by(df$value,list(df$i,df$j),quantile) > q : i1 : j1       0%      25%      50%      75%     100% 18.00297 28.24206 33.22294 62.73622 85.77076 ------------------------------------------------------------------------ : i2 : j1       0%      25%      50%      75%     100% 11.93874 29.25933 55.02918 69.63326 83.25839 ------------------------------------------------------------------------ : i1 : j2        0%       25%       50%       75%      100% 7.644259 41.375281 61.313381 82.127640 99.813645 ------------------------------------------------------------------------ : i2 : j2       0%      25%      50%      75%     100% 16.85710 28.78330 53.35881 65.53476 87.37377 ------------------------------------------------------------------------ : i1 : j3       0%      25%      50%      75%     100% 14.01767 16.55922 30.77533 38.48282 67.22393 ------------------------------------------------------------------------ : i2 : j3       0%      25%      50%      75%     100% 5.60860 17.43385 33.31175 45.56650 64.92699 >

We can compute all quantiles in a oneliner script. Although the results are correct, it would be nice if we can organize the results in a proper data frame for further processing. Here is my first attempt:

> cbind(unique(df[,1:2]), +       matrix(unlist(q),ncol=5,byrow=T,dimnames=list(NULL,names(q[[1,1]]))) + )     i  j        0%      25%      50%      75%     100% 1  i1 j1 18.002966 28.24206 33.22294 62.73622 85.77076 9  i1 j2 11.938737 29.25933 55.02918 69.63326 83.25839 17 i1 j3  7.644259 41.37528 61.31338 82.12764 99.81365 25 i2 j1 16.857103 28.78330 53.35881 65.53476 87.37377 33 i2 j2 14.017667 16.55922 30.77533 38.48282 67.22393 41 i2 j3  5.608600 17.43385 33.31175 45.56650 64.92699

The unique function populates the leftmost columns i and j. The unlist operation will make one long vector of the quantiles. We make the matrix of the correct size by using ncol=5. Finally, as the column names are lost in this operation, we re-establish them by dimnames.

Unfortunately this is wrong as unique populates the columns i and j in a different order than by does. The easy fix is as follows:

> q <- by(df$value,list(df$j,df$i),quantile) > qdf <- cbind(unique(df[,1:2]), +          matrix(unlist(q),ncol=5,byrow=T,dimnames=list(NULL,names(q[[1,1]]))) +       ) > qdf     i  j        0%      25%      50%      75%     100% 1  i1 j1 18.002966 28.24206 33.22294 62.73622 85.77076 9  i1 j2  7.644259 41.37528 61.31338 82.12764 99.81365 17 i1 j3 14.017667 16.55922 30.77533 38.48282 67.22393 25 i2 j1 11.938737 29.25933 55.02918 69.63326 83.25839 33 i2 j2 16.857103 28.78330 53.35881 65.53476 87.37377 41 i2 j3  5.608600 17.43385 33.31175 45.56650 64.92699

The change is in the order of list(df$j,df$i) passed on to by. Now unique and q have the same order.

Performance

If we change the data to:

set i /i1*i50/     j /j1*j50/     k /k1*k50/     q /'0%', '25%', '50%', '75%', '100%'/ ;

we see the following timings:

• GAMS rank loop: 161 seconds
• R calculate quantiles: 1.5 seconds
• R convert results into a data frame: 0.08 seconds

Element renaming

In GAMS set elements (used for indexing) are strings. This means it is not necessary to use non-descriptive numbers as in:

set i /1*50/;

Actually, even with numbered elements, it is good practice to prefix such elements, as in:

set i /i1*i50/;

This will help when inspecting large, multi-dimensional data-structures.

As an example consider the code fragment:

set

 i /1*8/

 j /1*4/

 k /1*6/

;

parameter p(i,j,k);

p(i,j,k) = uniform(0,1);

I see these type of numbered sets quite often. Here is why this may be a bad idea: especially after pivoting a bit, it is really difficult to see what is going on:

If we prefix the set elements with a name, we achieve much more clarity:

When I receive a model with numbered set elements the first thing I try to do is to prefix the labels.

Of course, in practical models, we often can use meaningful names, but sometimes we just have numbers as IDs. I have also seen cases where the ID is a number or a very long description. In that case, the number can be a better candidate for set element.

Exporting set elements

In GAMS a parameter is really a sparse (multi-dimensional) matrix. This corresponds quite nicely to a table in a RDBMS (relational database) or to a dataframe in R or Python. A matrix in R or Python is dense, so those data structures do not work as well for large, sparse parameters, although R has nice facilities to use matrices with row and column names.

Example

Let’s generate some data in GAMS and export to an SQLite database:

set

 i 'observations' /i1*i50/

 s 'scenarios' /scen1*scen3/

;

parameter results(s,i);

results('scen1',i) = normal(0,1);

results('scen2',i) = normal(-1+4*ord(i)/card(i),0.4);

results('scen3',i) = uniform(4,5);

execute_unload "results.gdx",results;

execute "gdx2sqlite -i results.gdx -o results.db";

We can try to plot this data in R as follows:

library("RSQLite")
library("ggplot2")
 
# connect to the sqlite file
sqlite = dbDriver("SQLite")
db = dbConnect(sqlite, "results.db")

# retrieve data
df = dbGetQuery(db,"select * from results")
head(df)
 
ggplot(data=df,aes(x=i,y=value,color=s)) + geom_line()
 
dbDisconnect(db)

The output of head(df) indicates we read the data correctly:

## head(df)

     s  i      value

1 scen1 i1 -0.3133429

2 scen1 i2  0.3276748

3 scen1 i3  0.4635588

4 scen1 i4 -1.8299478

5 scen1 i5 -0.7316124

6 scen1 i6 -0.9715983

Unfortunately, things don’t work exactly as expected. We see an error message:

geom_path: Each group consists of only one observation. Do you need to adjust

the group aesthetic?

And the plot does not look correct:

The reason is the variable df$i.If we convert this back to integers we are ok. To do this conversion we actually have to do two things:

1. Drop the first character ‘i’ from df$i.
2. Convert df$i from string to numeric (type casting).

Interestingly we can fix this in several stages:

Partial GAMS solution

We can not really rename set elements in GAMS. However we can create a new parameter, a mapping set and use a sum statement to map between sets:

set

 i 'observations' /i1*i50/

 s 'scenarios' /scen1*scen3/

 i2 'observations as number' /1*50/

 map(i2,i) '1-1 mapping set'

;

map(i2,i)$(ord(i2)=ord(i)) = yes;

parameter results(s,i);

results('scen1',i) = normal(0,1);

results('scen2',i) = normal(-1+4*ord(i)/card(i),0.4);

results('scen3',i) = uniform(4,5);

parameter results2(s,i2);

results2(s,i2) = sum(map(i2,i),results(s,i));

Now we can export the parameter results2 instead of results. Note that index i2 will still be a string when it arrives in R. The type casting to a numeric type has to be done in R or in SQL (see next sections how this can be done).

Conversion in R

The conversion to a numeric type can be done in R in just one line

df$i <- as.numeric(substring(df$i,2))

Now things look better:

Conversion in SQL

We can also perform the conversion in SQL while reading the data from the SQLite database file:

select s,
      cast(substr(i,2) as numeric) as i,
      value
from results

Scripting

This whole thing can be easily scripted in GAMS.

The CDF of the Gamma distribution in GAMS

For a model (1) I needed to calculate the cumulative distribution function \(F(x)\) of the Gamma distribution. There is no uniform consensus:  the Gamma distribution has different parametrizations. I use the \(k, \theta\) parametrization from (2). This would yield (3):

\[F(x) =  \gamma( \frac{x}{\theta},k)\]

where \(\gamma(x,a)\) is the incomplete Gamma function defined by:

\[\gamma(x,a) = \frac{1}{\Gamma(a)} \int_0^x t^{a-1}e^{-t} dt\]

The function \(\Gamma(x)\) is the (complete) Gamma function:

\[\Gamma(x) = \int_0^{\infty} t^{x-1} e^{-t} dt \]

Unfortunately there a few alternative definitions of the incomplete Gamma function. In (2) the following definition is used:

\[\gamma(s,x) = \int_0^x t^{s-1} e^{-t} dt\]

So watch out: different definitions are being used.

As  a result we can implement the CDF of the Gamma distribution with parameter \(k\) and \(\theta\) as follows (2):

GammaReg(x/theta,k)

There is also a more obscure way to calculate this using an extrinsic function. The syntax is unfortunately needlessly complicated and the documentation is really, really bad (4).

scalars   k / 0.363 /   theta / 27.863 /   x / 1.5 /   p ; p = gammareg(x/theta,k); display p; $FuncLibIn stodclib stodclib function cdfGamma / stodclib.cdfGamma /; p = cdfGamma(x,k,theta); display p;

References

1. Modeling Flood Damages, http://yetanothermathprogrammingconsultant.blogspot.com/2017/07/modeling-flood-damages.html
2. Gamma Distribution, https://en.wikipedia.org/wiki/Gamma_distribution
3. New Special Functions in GAMS, http://www.amsterdamoptimization.com/pdf/specfun.pdf
4. Stochastic Library, https://www.gams.com/latest/docs/userguides/userguide/_u_g__extrinsic_functions.html#UG_ExtrinsicFunctions_StochasticLibrary

no loops please

In (1) a piece of matlab code is presented:

There may be better way to formulate this in Matlab.  In many languages, including Matlab, R and GAMS loops are considered bad if there are alternatives. To be honest, I don’t see immediately how this particular fragment can be vectorized in Matlab.

The suggested  GAMS version is:

I strongly disagree. I believe this can be formulated without loops as:

 c(i,j,k)$(ord(k)<=plan(j,i)) = sub(j,i); 

Also notice that MATLAB will store this in a fully allocated or dense matrix (unless explicitly using a sparse matrix) while GAMS always stores things in a sparse format.

References

1. Generate Data using for loop in GAMS. https://stackoverflow.com/questions/45167648/generate-data-using-for-loop-in-gams

Ubuntu Bash on Windows

Windows 10 has a beta feature: Bash.exe with an Ubuntu Linux subsystem (see (1) and (2)), also known as WSL (Windows Subsystem for Linux). This will allow you to run Linux command line tools from within Windows. This is not the same as running Linux in a Virtual Machine (VM) using virtualization software like VirtualBox or VMware. It is also different from Cygwin or MINGW which is a windows port of GNU tools. In the case of  WSL, real, unmodified Ubuntu Linux binaries are executed.





There is a single Linux instance, so when invoking bash.exe several times, you talk to the same instance:

The main issue I noticed is that some networking commands (like ping) requires bash.exe to be started as administrator.

References

1. Jack Hammons, Bash on Ubuntu on Windows, https://msdn.microsoft.com/commandline/wsl
2. Learning about Bash on Windows Subsystem for Linux, https://blogs.msdn.microsoft.com/commandline/learn-about-bash-on-windows-subsystem-for-linux/
3. VirtualBox, https://www.virtualbox.org/
4. VMware, https://www.vmware.com/
5. Cygwin, https://www.cygwin.com/
6. MINGW, Minimalist GNU for Windows, http://www.mingw.org/

Modeling flood damages

While working on a investment planning model to combat damages resulting from flooding, I received the results from a rainfall model that calculates damages as a result of excessive rain and flooding. The picture the engineers produced looks like:

 

This picture has the scenario on the x-axis (sorted by damage) and the damages on the y-axis. This picture is very much like a load-duration curve in power generation.

For a more “statistical” picture we can use standard histogram (after binning the data):

Gamma distribution

We can use standard techniques to fit a distribution. When considering a Gamma distribution (1), a simple approach is the method of moments. The mean and the variance of the Gamma distribution with parameters \(k\) and \(\theta\) are given by:

\[\begin{align} &E[X] = k \cdot \theta\\ & Var[X] = k \cdot \theta^2 \end{align} \]

Using sample mean \(\mu\) and standard deviation \(\sigma\), we can solve:

\[\begin{align} & k \cdot \theta = \mu\\ & \sqrt{k} \cdot \theta = \sigma \end{align} \]

This can be easily solved numerically and it actually seems to work:

Weibull distribution

An alternative distribution that is sometimes suggested is the Weibull distribution (2).  The method of moments estimator for the Weibull distribution with parameters \(\lambda\) and \(k\) can be found by solving the system:

\[\begin{align} & \lambda \Gamma(1+1/k) = \mu\\ &\lambda \sqrt{\Gamma(1+2/k)+\left(\Gamma(1+1/k)\right)^2} = \sigma \end{align}\]

I was unable to get a solution from this: solvers failed miserably on this.

An alternative approach would be to use an MLE (Maximum Likelihood Estimation) technique. This yields a system of equations:

\[\begin{align} & \lambda^k = \frac{1}{n}\sum_{i=1}^n x_i^k\\ &k^{-1} = \frac{\sum_{i=1}^n x_i^k \ln x_i}{\sum_{i=1}^n x_i^k} – \frac{1}{n}\sum_{i=1}^n \ln x_i \end{align}\]

Note that we can solve the second equation first to solve for \(k\) and then calculate \(\lambda\) using the first equation. A solver like CONOPT will do this automatically by recognizing the triangular structure of this system. Note that our data set contains many \(x_i=0\). These are replaced by \(x_i=0.001\) so we can take the logarithm.

This gives is a very similar picture:

References

1. Gamma distribution, https://en.wikipedia.org/wiki/Gamma_distribution
2. Weibull distribution, https://en.wikipedia.org/wiki/Weibull_distribution

Rectangles: no-overlap constraints

The question of how we can formulate constraints that enforce rectangles not to overlap comes up regularly (1). There are basically four cases to consider when considering two rectangles \(i\) and \(j\):

In the picture above, we indicate by \((x_i,y_i)\) the left-lower corner of a rectangle (these are typically decision variables). The height and the width are \(h_i\) and \(w_i\) (these are typically constants). From the above, we can formulate the “no-overlap”  constraints as:

\[ \bbox[lightcyan,10px,border:3px solid darkblue] { \begin{align} &x_i+w_i \le x_j  & \text{or} \\ &x_j+w_j \le x_i  & \text{or} \\ &y_i+h_i \le y_j  & \text{or} \\ &y_j+h_j \le y_i \end{align} }\]

The “or” condition can be modeled with the help of binary variables \(\delta\) and big-M constraints:

\[ \bbox[lightcyan,10px,border:3px solid darkblue] { \begin{align} &x_i+w_i \le x_j + M_1 \delta_{i,j,1}\\ &x_j+w_j \le x_i + M_2 \delta_{i,j,2}\\ &y_i+h_i \le y_j + M_3 \delta_{i,j,3}\\ &y_j+h_j \le y_i+ M_4 \delta_{i,j,4}\\ &\sum_k \delta_{i,j,k} \le 3\\ & \delta_{i,j,k} \in \{0,1\} \end{align} }\]

The binary variables make sure at least one of the constraints is active (not relaxed).

If we have multiple rectangles we need to compare all combinations, i.e. we generate these constraints for  \(\forall i,j\). Except: we can not compare a rectangle with itself. A first approach would be to generate the above constraints for all \(i\ne j\). This is correct, but we can apply an optimization. Note that we actually compare things twice: if rectangles \(i\) and \(j\) are non-overlapping then we don’t need to check the pair of rectangles \(j\) and \(i\). That means we only need constraints and variables \(\delta_{i,j,k}\) for \(i<j\).

It is important to make constants \(M_k\) as small as possible. In general we have a good idea about good values for these \(M_{k}\)’s. E.g. consider the case where we need to pack boxes in a container (2):

If the container has length and width \(H\) and \(W\) then we can easily see: \(M_1=M_2=W\) and \(M_3=M_4=H\).

References

1. How to write this constraint in LPSolve?, https://stackoverflow.com/questions/44858353/how-to-write-this-constraint-in-lp-solver-logical-and-and-doesnot-exist
2. Filling up a container with boxes, http://yetanothermathprogrammingconsultant.blogspot.com/2016/06/filling-up-container-with-boxes.html

GLPK was used in the proof of the 300 year old Kepler conjecture

A lot of authors.

Apparently 43,078 LPs had to be solved.

References

1. Thomas Hales e.a., A Formal Proof of the Kepler Conjecture, Forum of mathematics, Pi (2017), vol 5, e2, https://www.cambridge.org/core/services/aop-cambridge-core/content/view/78FBD5E1A3D1BCCB8E0D5B0C463C9FBC/S2050508617000014a.pdf/formal_proof_of_the_kepler_conjecture.pdf
2. Original message appeared in the Help-GLPK mailing list,  https://lists.gnu.org/archive/html/help-glpk/2017-07/msg00001.html