### Table of Contents

# Introduction

Functions play an important role in the GAMS language, especially for non-linear models. Similarly to other programming languages, GAMS provides a number of built-in (intrinsic) functions. However, GAMS is used in an extremely diverse set of application areas and this creates frequent requests for the addition of new and often sophisticated and specialized functions. There is a trade-off between satisfying these requests and avoiding complexity not needed by most users. The GAMS Function Library Facility ( Functions ) provides the means for managing that trade-off. In this Appendix the extrinsic function libraries that are included in the GAMS distribution are described. In addition, we provide some pointers for users wishing to build their own extrinsic function library.

In the tables that follow, the Endogenous Classification (second column) specifies in which models the function can legally appear with endogenous (non-constant) arguments. In order of least to most restrictive, the choices are `any`

, `NLP`

, `DNLP`

or `none`

.

The following conventions are used for the function arguments. Lower case indicates that an endogenous variable is allowed. Upper case indicates that a constant argument is required. The arguments in square brackets can be omitted and default values will be used.

# Fitpack Library

FITPACK by Paul Dierckx ^{1)}
is a FORTRAN based library for one and two dimensional spline interpolation. This library has been repackaged to work with the GAMS Function Library Facility. As it can be seen in the GAMS Test Library model `fitlib01`

the function data needs to be stored in a GDX file `fit.gdx`

containing a three dimensional parameter `fitdata`

. The first argument of that parameter contains the function index, the second argument is the index of the supporting point and the last one needs to be one of `w`

(weight), `x`

(x-value), `y`

(y-value) or `z`

(z-value).

Function | Endogenous Classification | Description |
---|---|---|

`fitFunc(FUNIND,x[,y])` | DNLP | Evalute Spline |

`fitParam(FUNIND,PARAM[,VALUE])` | none | Read or set parameters |

The function `FitParam`

can be used to change certain parameters used for the evaluation:

- 1: Smoothing factor (S)
- 2: Degree of spline in direction x (Kx)
- 3: Degree of spline in direction y (Ky)
- 4: Lower bound of function in direction x (LOx)
- 5: Lower bound of function in direction y (LOy)
- 6: Upper bound of function in direction x (UPx)
- 7: Upper bound of function in direction y (UPy)

This library is made available by the following directive:

$FuncLibIn <InternalLibName> fitfclib

# Piecewise Polynomial Library

This library can be used to evaluate piecewise polynomial functions. The functions which should be evaluated need to be defined and stored in a GDX file like it is done in the GAMS Test Library model [pwplib01]:

* Define two piecewise polynomial functions Table pwpdata(*,*,*) '1st index: function number, 2nd index: segment number, 3rd index: degree' leftBound 0 1 2 1.1 1 2.4 -2.7 0.3 1.2 4 5.6 -4.3 0.5 2.1 0 0 -6.3333 0 2.2 0.3333 1.0370 -12.5554 9.3333 2.3 0.6667 9.7792 -38.7791 29 ; * Write pwp data to gdx file read by external library $gdxout pwp.gdx $unload pwpdata $gdxout

On each row of the Table pwpdata we have

FuncInd.SegInd leftBound Coef0 Coef1 Coef2 ...

`FuncInd`

sets the function index. `SegInd`

defines the index of the segment (or interval) which is decribed here. `LeftBound`

gives the lower bound of the segment. The upper bound is the lower bound on the next row, or infinity if this is the last segment. `CoefX`

defines the Xth degree coefficient of the polynomial corresponding to this segment.

This library is made available by the following directive:

$FuncLibIn <InternalLibName> pwpcclib

**Table 2: Piecewise polynomial functions**

Function | Endogenous Classification | Description |
---|---|---|

`pwpFunc(FUNIND,x)` | DNLP | Piecewise Polynomials |

# Stochastic Library

The stochastic library provides random deviates, probability density functions, cumulative density functions and inverse cumulative density functions for certain distributions. This library is made available by the following directive:

$FuncLibIn <InternalLibName> stodclib

**Table 3: Random number generators**

Function | Description |
---|---|

`SetSeed(SEED)` | defines the seed for random number generator |

Function | Description |
---|---|

`beta(SHAPE_1,SHAPE_2)` | Beta distribution with shape parameters SHAPE_1 and SHAPE_2, see MathWorld |

`cauchy(LOCATION,SCALE)` | Cauchy distribution with location parameter LOCATION and scale parameter SCALE, see MathWorld |

`ChiSquare(DF)` | Chi-squared distribution with degrees of freedom DF, see MathWorld |

`exponential(LAMBDA)` | Exponential distribution with rate of changes LAMBDA, see MathWorld |

`f(DF_1,DF_2)` | F-distribution with degrees of freedom DF_1 and DF_2, see MathWorld |

`gamma(SHAPE,SCALE)` | Gamma distribution with shape parameter SHAPE and scale parameter SCALE, see MathWorld |

`gumbel(LOCATION,SCALE)` | Gumbel distribution with location parameter LOCATION and scale parameter SCALE, see MathWorld |

`invGaussian(MEAN,SHAPE)` | Inverse Gaussian distribution with mean MEAN and scaling parameter SHAPE, see MathWorld |

`laplace(MEAN,SCALE)` | Laplace distribution with mean MEAN and scale parameter SCALE, see MathWorld |

`logistic(LOCATION,SCALE)` | Logistic distribution with location parameter LOCATION and scale parameter SCALE, see MathWorld |

`logNormal(LOCATION,SCALE)` | Log Normal distribution with location parameter LOCATION and scale parameter SCALE, see MathWorld |

`normal(MEAN,STD_DEV)` | Normal distribution with mean MEAN and standard deviation STD_DEV, see MathWorld |

`pareto(SCALE,SHAPE)` | Pareto distribution with scaling parameter SCALE and shape parameter SHAPE, see MathWorld |

`rayleigh(SIGMA)` | Rayleigh distribution with parameter SIGMA, see MathWorld |

`studentT(DF)` | Student's t-distribution with degrees of freedom DF, see MathWorld |

`triangular(LOW,MID,HIGH)` | Triangular distribution between LOW and HIGH, MID is the most probable number, see MathWorld |

`uniform(LOW,HIGH)` | Uniform distribution between LOW and HIGH, see MathWorld |

`weibull(SHAPE,SCALE)` | Weibull distribution with shape parameter SHAPE and scaling parameter SCALE, see MathWorld |

Function | Description |
---|---|

`binomial(N,P)` | Binomial distribution with number of trials N and success probability P in each trial, see MathWorld |

`geometric(P)` | Geometric distribution with success probability P in each trial, see MathWorld |

`hyperGeo(TOTAL,GOOD,TRIALS)` | Hypergeometric distribution with total number of elements TOTAL, number of good elements GOOD and number of trials TRIALS, see MathWorld |

`logarithmic(P-FACTOR)` | Logarithmic distribution with parameter P-FACTOR, also called log-series distribution, see MathWorld |

`negBinomial(FAILURES,P)` | Negative Binomial distribution with the number of failures until the experiment is stopped FAILURES and success probability P in each trial. The generated random number describes the number of successes until we reached the defined number of failures, see MathWorld |

`poisson(LAMBDA)` | Poisson distribution with mean LAMBDA, see MathWorld |

`uniformInt(LOW,HIGH)` | Integer Uniform distribution between LOW and HIGH, see MathWorld |

For each distribution in table Table 3, the library offers four functions in Table 4:

**Table 4: Distribution functions**

Function | Endogenous Classification | Description |
---|---|---|

`d<DistributionName>` | none | generates a random number |

`pdf<DistributionName>` | DNLP (none for discrete distributions) | probability density function |

`cdf<DistributionName>` | DNLP (none for discrete distributions) | cumulative distribution function |

`icdf<DistributionName>` | DNLP (none for discrete distributions) | inverse cumulative distribution function |

The function `d<DistributionName>`

needs the arguments described in table Table 4. The other functions get an additional argument at the first position: The point to evaluate. This parameter can be an endogenous variable. The following table shows all fours functions for the Normal distribution:

**Table 5: Normal distribution functions**

Function | Endogenous Classification | Description |
---|---|---|

`dNormal(MEAN,STD_DEV)` | none | generates a random number with Normal distribution |

`pdfNormal(x,MEAN,STD_DEV)` | DNLP | probability density function for Normal distribution |

`cdfNormal(x,MEAN,STD_DEV)` | DNLP | cumulative distribution function for Normal distribution |

`icdfNormal(x,MEAN,STD_DEV)` | DNLP | inverse cumulative distribution function for Normal distribution |

# LINDO Sampling Library

The LINDO Sampling Library provides samples of random numbers for certain distributions. It is made available by the following directive:

$FuncLibIn <InternalLibName> lsadclib

The function `sampleLS<DistributionName>`

creates a sample of the specified distribution according to the distribution parameters and returns a `HANDLE`

that references the sample as illustrated in the example down below. The Parameter `SAMSIZE`

must be specified and describes the size of the sample while `VARRED`

is optional and provides the possibility to define a variance reduction method (0=none, 1=Latin Hyper Square, 2=Antithetic). If omitted Latin Hyper Square Sampling is used.

**Table 6:LINDO sampling functions**

Function | Description |
---|---|

`getSampleValues(HANDLE)` | Retrieve sampling. |

`induceCorrelation(CORTYPE)` | Induce correlation that has to be set with `setCorrelation` before.`CORTYPE` describes the correlation type and must be one of 0 (PEARSON), 1 (KENDALL) or 2 (SPEARMAN). |

`setCorrelation(SAMPLE1,SAMPLE2,COR)` | Define correlation between two samplings. |

`setSeed(SEED)` | Define the seed for random number generator. |

`setRNG(RNG)` | Define the random number generator to use, possible values are -1 (FREE), 0 (SYSTEM), 1 (LINDO1), 2 (LINDO2), 3 (LIN1), 4 (MULT1), 5 (MULT2) and 6 (MERSENNE). |

Function | Description |
---|---|

`beta(SHAPE_1,SHAPE_2,SAMSIZE[,VARRED])` | Beta distribution specified by two shape parameters. |

`cauchy(LOCATION,SCALE,SAMSIZE[,VARRED])` | Cauchy distribution specified by location and scale parameter. |

`chisquare(DF,SAMSIZE[,VARRED])` | Chi-Squared distribution specified by degrees of freedom. |

`exponential(RATE,SAMSIZE[,VARRED])` | Exponential distribution specified by rate of change. |

`f(DF_1,DF_2,SAMSIZE[,VARRED])` | F distribution specified by degrees of freedom. Note that sampleLSgamma uses another version of the F distribution than dF from the Stochastic Library. |

`gamma(SHAPE,SCALE,SAMSIZE[,VARRED])` | Gamma distribution specified by shape and scale parameter. Note that the use of sampleLSgamma(A,B) is equivalent to the use of dGamma(B,A) from the Stochastic Library. |

`gumbel(LOCATION,SCALE,SAMSIZE[,VARRED])` | Gumbel distribution specified by location and scale parameter. |

`laplace(LOCATION,SCALE,SAMSIZE[,VARRED])` | Laplace distribution specified by location and scale parameter. |

`logistic(LOCATION,SCALE,SAMSIZE[,VARRED])` | Logistic distribution specified by location and scale parameter. |

`lognormal(LOCATION,SCALE,SAMSIZE[,VARRED])` | Log Normal distribution specified by location and scale parameter. |

`normal(MEAN,STD_DEV,SAMSIZE[,VARRED])` | Normal distribution specified by given mean and standard deviation. |

`pareto(SCALE,SHAPE,SAMSIZE[,VARRED])` | Pareto distribution specified by shape and scale parameter. |

`studentt(DF,SAMSIZE[,VARRED])` | Student's t-distribution specified by degrees of freedom. |

`triangular(LOW,MID,HIGH,SAMSIZE[,VARRED])` | Triangular distribution specified by lower and upper limit and mid value. |

`uniform(LOW,HIGH,SAMSIZE[,VARRED])` | Uniform distribution specified by the given bounds. |

`weibull(SCALE,SHAPE,SAMSIZE[,VARRED])` | Weibull distribution specified by scale and shape parameter. Note that the use of sampleLSweibull(A,B) is equivalent to the use of dWeibull(B,A) from the Stochastic Library . |

Function | Description |
---|---|

`binomial(N,P,SAMSIZE[,VARRED])` | Binomial distribution specified by number of trials `N` and success probability `P` in each trial. |

`hypergeo(TOTAL,GOOD,TRIALS,SAMSIZE[,VARRED])` | Hypergeometric distribution specified by total number of elements, number of good elements and number of trials. |

`logarithmic(P-FACTOR,SAMSIZE[,VARRED])` | Logarithmic distribution specified by `P-Factor` . Note that sampleLSlogarithmic uses another version of the logarithmic distribution than dLogarithmic from the Stochastic Library. |

`negbinomial(SUCC,P,SAMSIZE[,VARRED])` | Negative Binomial distribution specified by the number of successes and the probability of success. The generated random number describes the number of failures until we reached the defined number of successes. Note that the use of sampleLSnegbinomial(R,P) is equivalent to the use of dNegBinomial(R,P-1) from the Stochastic Library . |

`poisson(MEAN,SAMSIZE[,VARRED])` | Poisson distribution specified by mean. |

The following example illustrates the use of the sample generator and shows how the commands `setCorrelation`

and `induceCorrelation`

work:

$funclibin lsalib lsadclib function normalSample /lsalib.SampleLSnormal / getSampleVal /lsalib.getSampleValues / setCor /lsalib.setCorrelation / indCor /lsalib.induceCorrelation /; scalar d,h,k; h = normalSample(5,2,12); k = normalSample(5,2,12); set i /i01*i12/; parameter sv(i); loop(i, sv(i) = getSampleVal(k); ); display sv; loop(i, sv(i) = getSampleVal(h); ); display sv; d=setCor(h,k,-1); d=indCor(1); loop(i, sv(i) = getSampleVal(k); ); display sv; loop(i, sv(i) = getSampleVal(h); ); display sv;

The resulting output shows that the values of `sv`

are restructured according to the desired correlation:

---- 18 PARAMETER sv i01 7.610, i02 5.710, i03 3.755, i04 5.306, i05 2.382, i06 2.174 i07 6.537, i08 4.975, i09 8.260, i10 3.067, i11 6.216, i12 4.349 ---- 22 PARAMETER sv i01 7.610, i02 5.710, i03 3.755, i04 5.306, i05 2.382, i06 2.174 i07 6.537, i08 4.975, i09 8.260, i10 3.067, i11 6.216, i12 4.349 ---- 28 PARAMETER sv i01 2.382, i02 4.349, i03 6.216, i04 4.975, i05 7.610, i06 8.260 i07 3.067, i08 5.306, i09 2.174, i10 6.537, i11 3.755, i12 5.710 ---- 32 PARAMETER sv i01 7.610, i02 5.710, i03 3.755, i04 5.306, i05 2.382, i06 2.174 i07 6.537, i08 4.975, i09 8.260, i10 3.067, i11 6.216, i12 4.349

# Build Your Own: Trigonometric Library Example

This library serves as an example of how to code and build an extrinsic function library. The library is included in the GAMS distribution in binary form and also as source code written in C, Delphi, and FORTRAN, respectively, that comes along with the GAMS Test Library models [trilib01], [trilib02], and [trilib03]. Table Table 7 lists the extrinsic functions that are implemented by the libraries.

**Table 7: Trigonometric functions**

Function | Description | End. Classif. |
---|---|---|

`setMode(MODE)` | sets mode globally, could still be overwritten by MODE at (Co)Sine call, possible values are 0=radians and 1=degree | none |

`cosine(x[,MODE])` | returns the cosine of the argument x, default setting: MODE = 0 | NLP |

`sine(x[,MODE])` | returns the sine of the argument x, default setting: MODE = 0 | NLP |

`pi` | value of \(\pi=3.141593...\) | any |

The C implementation of this extrinsic function library can be found in the files `tricclib.c`

and `tricclibql.c`

. Together with the API specification file `extrfunc.h`

, these files document the callbacks that need to implemented by a GAMS extrinsic function library.

Note that the file `tricclibql.c`

, which implements the `querylibrary`

callback, has been generated by the Python helper script `ql.py`

. The purpose of the `querylibrary`

callback is to provide information about the library itself and the extrinsic functions it implements to the GAMS execution system. For example, the information that the `cosine`

function has an endogenous required first argument and an exogenous optional second argument is available from the `querylibrary`

callback. The `ql.py`

file that generated the file `tricclibql.c`

(and also the Delphi and Fortran90 equivalents `tridclibql.inc`

and `triifortlibql.f90`

) does so by processing the specification file `tri.spec`

. See the comments in this file for info on how to write such a specification. The `ql.py`

script is invoked by executing `./ql.py tri.spec`

from the shell.

When implementing a function in a library this function needs to return the function value of the input point. Furthermore, the solver might need derivative information at the input point. Therefore, the implementation of the function needs to be able to return gradient and Hessian values. This can sometimes be inconvenient. GAMS can use the function values at points close to the input point to estimate the gradient and Hessian values using finite differences. However, this is not as accurate as analytic derivatives and requires a number of function evaluations, so the convenience comes at a price. The attribute `MaxDerivative`

in the specification of a function signals GAMS the highest derivatives this function will provide. For higher order derivatives GAMS will use finite differences to approximate the derivative values.

# Build Your Own: Reading the GAMS Parameter File in your Library

This library serves as an example of how to code and build an extrinsic function library that reads the information from the GAMS parameter file. The library is included in the GAMS distribution in binary form and also as source code written in C, that comes along with the GAMS Test Library model [parlib01]. Table Table 8 lists the extrinsic functions that are implemented by the libraries.

**Table 8: GAMS Parameter File reading**

Function | Endogenous Classification | Description |
---|---|---|

`LogOption` | any | Return the value for LogOption found in the parameter file |

The C implementation of this extrinsic function library can be found in the files `parcclib.c`

and `parcclibql.c`

. Together with the API specification file `extrfunc.h`

. The reading of the GAMS Parameter file is done in the function LibInit.

# CPP Library

This library serves both as an example of how to use C++ to obtain gradients and Hessians ``for free'' and as a source of functions based on the multivariate normal distribution. The library is available in compiled form and as C++ source. Test Library model `cpplib00`

exercises the process of building a shared library from C++ source and doing some basic tests, while models [cpplib01] through [cpplib05] are more thorough tests for the CPP library extrinsics shipped with the distribution. These functions are listed and described in Table Table 9.

**Table 9: CPP Library functions**

Function | Endogenous Classification | Description |
---|---|---|

`pdfUVN(x)` | NLP | PDF of uni-variate normal, see MathWorld or R |

`cdfUVN(x)` | NLP | CDF of uni-variate normal, see MathWorld or R |

`pdfBVN(x,y,r)` | NLP | PDF of bivariate normal, see MathWorld or R |

`cdfBVN(x,y,r)` | NLP | CDF of bivariate normal, see MathWorld or R |

`pdfTVN(x,y,z,r21,r31,r32)` | NLP | PDF of trivariate normal, see MathWorldor R |

## Automatic differentiation

Often extrinsic functions are created so that they can be used with endogenous arguments. In such cases, it is necessary that we provide first and second derivatives w.r.t. these arguments in addition to the function values themselves. One way to compute these derivatives is via automatic differentiation techniques (see the Wikipedia article for details). With C++, it is possible to overload the usual arithmetic operators (assignment, addition, multiplication, etc.) so that this automatic differentiation occurs with little or no change to the function-only source code. This is the technique used to compute the derivatives in the CPP Library. Test Library model `cpplib00`

includes all the source code for the CPP Library and illustrates/exercises the steps to build the library from this source. The source is a working self-documentation of how the process of automatic differentiation works.

## Multi-variate Normal Distributions

The CPP Library implements the PDF and CDF for the univariate, bivariate, and trivariate standard normal distributions. We use the standard normal (mean of 0, standard deviation of 1) since intrinsic functions are limited to 20 arguments. The functions for the univariate case are included as convenient examples and should give results (nearly) identical to the `pdfNormal`

and `cdfNormal`

functions from the stochastic extrinsic library LINDO Sampling Library . For the multivariate cases, we based our implementation on the TVPACK from Alan Genz ^{2)} , with modifications to allow for proper computation of derivatives. Again, the number of arguments allowed is a consideration, so we chose to implement functions taking correlation coefficients as arguments, not a covariance matrix. The conversion from the latter case to the former is relatively straightforward. Here we describe it with some `R`

code, making use of the `mnormt`

package that computes the multivariate CDF using similar code from Genz:

# start with a mean mu and variance-covariance matrix S1 x <- c(1.0,3.0,3.0) mu <- c(0.0,1.0,-1.0) S1 <- matrix(c(1.0,1,1.5, 1,4,1.5, 1.5,1.5,9),3,3) v1 <- pmnorm(x=x, mean=mu, varcov=S1) # convert to std normal with 3 correlation coeffs R <- cov2cor(S1) sd <- sqrt(diag(S1)) xn <- (x-mu) / sd v2 <- pmnorm(x=xn, mean=0, varcov=R)

For the bivariate case, there is one correlation coefficient, \(r\). The CDF implementation is not quite accurate to machine precision, having 14 or 15 digits of accuracy. The trivariate case includes 3 correlation coefficients, the 3 off-diagonal elements from the lower triange of \(R\) above. The accuracy of the CDF depends on the inputs: it is higher when the correlation is nearly zero and poorer as the condition number of \(R\) increases. Typically, an accuracy of \(10^{-11}\) is all that can be achieved. In both multivariate cases, you should avoid evaluating at or near points where the correlation matrix is degenerate. At nearly degenerate points, the accuracy of the distribution and density functions suffers, and when the correlation matrix becomes degenerate the distribution becomes so also.

# Extrinsic function vs. External Equation Comparison

For nonlinear models users can extend the solving capability of GAMS by using extrinsic functions and external equations. A function, for example, sin(x) is a intrinsic function because it is defined by GAMS. A function provided by the user is called a extrinsic function and is defined by the user, for example, a user implementation of sin(x) could be called mysin(x). An equality equation, here named e1, can be written as follows in GAMS: e1.. sin(x)=E=1;. However, if the equation evaluation is not done by GAMS then it is called a external equation. The external equation evaluation routine is provided by the user and a solution for a model must satisfy all internal and external equations. A external equation is denoted by equation type =X=.

Example. A possible GAMS syntax of an equation that integrates x from x1 to x2: as a extrinsic function: e1.. integralx(x1,x2) =E= z; as a external equation: e1.. 1*x1 + 2*x2 + 3*z =X= 1;

Recall that the integration can be analytically solved as follows: integralx(x1,x2)=0.5*(sqr(x2)-sqr(x1))

Note that then numbers 1 to 3 for the external equation denote positions while the right hand side (RHS) denote the external equation number.

Some characteristics and clarifications

Characteristic | Extrinsic function | External equation |
---|---|---|

Argument limitation | 10 | No |

Available in statements | Yes | No |

Debugging support | Yes | No |

Returns Hessian to solver | Yes | No |

- An example of a extrinsic function with 10 arguments: integral(x1,x2,x3x,x4,x5,x6,x7,x8,x9,x10)
- An example of a data statement where we use a external function to set a variable level: x1.l=sqr(integralx(0,2));
- Debugging support: support for numerical gradient and hessian
- user provided: integralx.grad(2:0,2) vs. numerically provided by GAMS: integralx.gradn(2:0,2)

- The user can not provide a solver with Hessian information for external equations, however, many solvers can solve without Hessian information, for example,by approximating the Hessian.

**Additional notes**:

- CONOPT verifies by a numerical estimation of the derivatives that the user has provided reasonable derivative information (a unreasonable derivative results in Fatal error).
- Global solvers (ANTIGONE, BARON, Couenne, LindoGlobal, SCIP) analyze the algebra and can therefore not handle extrinsic functions or external equations.

**Examples of use**:

- Extrinsic functions can be used for special purposes, see Fitpack, Piecewise Polynomial, Stochastic Library, LINDO Sampling, Trigonometric libraries in the GAMS user's guide
- External equation provides a unrestricted way to solve custom equations and equation systems, for example, a partial differential equation system.

For more information

- about extrinsic functions see Extrinsic Functions.
- about external equations see External Equations.

^{1)} Paul Dierckx, Curve and Surface Fitting with Splines, Oxford
University Press, 1993, http://www.netlib.org/dierckx/

^{2)} Alan Genz, http://www.math.wsu.edu/faculty/genz/homepage