## Friday, March 10, 2017

### ifthen adventures in GAMS (and R)

GAMS has a function ifthen which can be used in equations. For example:

e.. y =e= ifthen(x<1,(1+x)/2,sqrt(x));

This function implements:

 $y=\begin{cases}\displaystyle\frac{1+x}{2}&\text{if x<1}\\ \sqrt{x}&\text{if x\ge 1}\end{cases}$

This is a perfectly smooth function: it is continuous and differentiable.

However when we use it in GAMS we are in for a surprise:

 variable x,y;equation e;e.. y =e= ifthen(x<1,(1+x)/2,sqrt(x));x.fx = -1;model m /all/;solve m minimizing x using dnlp;

Here we fix $$x=-1$$ so we just want to evaluate at that point. We see:

**** Exec Error at line 4: sqrt: FUNC DOMAIN: x < 0

The underlying reason is that GAMS always evaluates both branches. Obviously we would much prefer it would not do that.

If we replace fixing the variable by:

 x.lo = -1;x.l = 2;

i.e. start at $$x=2$$ and move to the optimal value $$x=-1$$, we see:

 C O N O P T 3   version 3.17A   Copyright (C)   ARKI Consulting and Development A/S                   Bagsvaerdvej 246 A                   DK-2880 Bagsvaerd, Denmark   Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK     0   0        1.4142135624E+00 (Input point)                   Pre-triangular equations:   0                  Post-triangular equations:  1      1   0        0.0000000000E+00 (After pre-processing)     2   0        0.0000000000E+00 (After scaling) ** Feasible solution. Value of objective =    2.00000000000   Iter Phase Ninf     Objective     RGmax    NSB   Step InItr MX OK    11   3        6.5613264487E-03 1.0E+00     1 3.9E-03      F  F    21   3        1.6608134823E-07 1.0E+00     1 6.0E-08      F  F    31   3        3.2205059005E-13 1.0E+00     1 1.5E-11      F  T    41   3        8.5912908646E-17 1.0E+00     1 2.2E-16      F  T    51   3        5.0871607672E-23 1.0E+00     1 3.4E-21      F  T    52   3        5.0871607672E-23 1.0E+00     1 ** Feasible solution. Convergence too slow. The change in objective   has been less than 3.0000E-12 for 20 consecutive iterations

The solver gets stuck at $$x=0$$. This is obviously incorrect. The only hint of the source of the problem is listed here:

 RESOURCE USAGE, LIMIT          0.062      1000.000ITERATION COUNT, LIMIT        52    2000000000EVALUATION ERRORS            803             0

This behavior is really a bug: we should see good feedback where the problem with domain errors originates from. Note that setting the limit of allowed evaluation errors (using option domlim) does not help in this case.

##### R

Interestingly R has something similar. Besides the standard if statement it has a vectorized ifelse function. The standard if is not vectorized which we can see from:

 > curve(if (x<1) (1+x)/2 else sqrt(x),-5,5) Warning message: In if (x < 1) (1 + x)/2 else sqrt(x) :   the condition has length > 1 and only the first element will be used

This only plots $$f(x)=\displaystyle\frac{1+x}{2}$$:

Basically this if statement is interpreted as:

if (x[1]<1) (1+x)/2 else sqrt(x)

With ifelse we get better results:

 > curve(ifelse(x<1,(1+x)/2,sqrt(x)),-5,5)Warning message:In sqrt(x) : NaNs produced

The plot is correct, but from the warning message we see that ifelse is also evaluating the $$\sqrt{x}$$ for all $$x$$. The ifelse function evaluates both branches.

To be precise this is only for vector arguments with some but not all $$x_i<1$$. For scalar arguments (vectors of length one) or more general if all elements of the vector are either less than one or all are $$\ge 1$$ things are ok:

Basically this is due to lazy but vector “all-or-nothing” evaluation.

##### Workaround

Of course the workaround is easy in this particular case. We can use something like:

e.. y =e= ifthen(x<1,(1+x)/2,sqrt(abs(x)));

Now the solver has no problems with this model, and can solve it very quickly (again with starting point $$x=2$$ and optimal solution $$x=-1$$):

 C O N O P T 3   version 3.17A   Copyright (C)   ARKI Consulting and Development A/S                   Bagsvaerdvej 246 A                   DK-2880 Bagsvaerd, Denmark   Iter Phase Ninf   Infeasibility   RGmax    NSB   Step InItr MX OK     0   0        1.4142135624E+00 (Input point)                   Pre-triangular equations:   0                  Post-triangular equations:  1      1   0        0.0000000000E+00 (After pre-processing)     2   0        0.0000000000E+00 (After scaling) ** Feasible solution. Value of objective =    2.00000000000   Iter Phase Ninf     Objective     RGmax    NSB   Step InItr MX OK     4   3       -1.0000000000E+00 0.0E+00     0 ** Optimal solution. There are no superbasic variables.

The same trick can be used in the R code:

curve(ifelse(x<1,(1+x)/2,sqrt(abs(x))),-5,5)

does not give any errors or warnings.