Extrinsic Functions

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

Table 1: Fitpack functions

Function Endogenous Classification Description
fitFunc(FUNIND,x[,y])DNLP Evalute Spline
fitParam(FUNIND,PARAM[,VALUE])none Read or set parameters

The function FitParamcan 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

Continous distributions

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

Discrete distributions

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

FunctionEndogenous Classification Description
d<DistributionName>nonegenerates 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

FunctionEndogenous 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 setCorrelationbefore.CORTYPEdescribes 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). Continuous distributions 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 . Discrete distributions Function Description binomial(N,P,SAMSIZE[,VARRED]) Binomial distribution specified by number of trials Nand success probability Pin 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

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

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

FunctionEndogenous 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

FunctionEndogenous 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
• 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.