Table of Contents
Introduction
Mathematical functions play an important role in the GAMS language, especially for nonlinear models. Like other programming languages, GAMS provides a number of built-in or intrinsic functions. 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 provides the means for managing this trade-off, since it allows users to import functions from an external library into a GAMS model. However, these external libraries can currently only provide functionality for the evaluation of functions (incl. their first and second derivatives) in a point. Solvers that need to analyze the algebraic structure of the model instance are therefore not able to work with extrinsic functions. This includes the class of deterministic global solvers, see column "Global" in this table, while, for example, stochastic global solvers can work with extrinsic functions.
In this chapter we will demonstrate how to access functions from an extrinsic function library in a GAMS model and we will describe the extrinsic function libraries that are included in the GAMS distribution. In addition, we will provide some pointers for users who wish to build their own extrinsic function library.
Using Function Libraries
Function libraries are made available to a model with the following compiler directive:
$funcLibIn <InternalLibName> <ExternalLibName>
Here InternalLibName
is a handle that will be used to refer to the library inside the model source code, ExternalLibName
is the file name of the shared library that implements the extrinsic functions. If no path is given, GAMS will look for the library in the directory extrinsic_functions
in the GAMS standard locations and in the GAMS system directory. To access a library that does not reside in these standard places, the external name should include a relative or absolute path to the location of the library. GAMS will then search for the specified library using the mechanisms specific to the host operating system. When processing the directive $funcLibIn, GAMS will validate the library, make the included functions available for use and add a table of the included functions to the listing file.
- Note
- The function library facility gives users complete control over naming so that potential name conflicts between libraries can be avoided.
Before the individual functions may be used, they have to be declared in the following way:
Function <InternalFuncName> /<InternalLibName>.<FuncName>/;
Here InternalFuncName
is the name of the individual function that will be used in the GAMS code. The user may choose this internal name freely and thus avoid potential naming conflicts. InternalLibName
is the name of the function library as defined by the $funcLibIn directive and FuncName
is the name of the individual function in the external function library. Once functions have been declared in this way they may be used like intrinsic functions.
Consider the following simple example:
$funcLibIn myLib tricclib
Function myCos /myLib.Cosine/
mySin /myLib.Sine/
myPi /myLib.Pi/;
Scalar d;
d = myCos(myPi/3);
display d;
Note that in the first line the external trigonometric library tricclib
is activated and the internal name myLib
is specified for it. Then the functions are declared. Observe that Cosine
, Sine
and Pi
are functions in the trigonometric library. After the library has been loaded and the functions have been declared, the functions may be used as usual. The trigonometric library is discussed in section Example: Trigonometric Library below.
Libraries that are included in the GAMS Distribution
In this section we will present the libraries that are included in the GAMS distribution: the Fitpack Library, the Piecewise Polynomial Library, the Stochastic Library and the LINDO Sampling Library. Note that the LINDO Sampling Library is only available for LINDO license holders. In addition, we will provide details on the Mutex Library.
In the tables that follow, the "Endogenous Classification" (last column) specifies in which models the function may legally appear. In order of least to most restrictive, the choices are any, DNLP, NLP, none. See section Classification of Models for details on model types in GAMS. Note well that functions classified as any are only permitted with exogenous (constant) arguments.
A word on the notation in the tables below: for function arguments, lower case indicates that an endogenous variable is allowed. For details on endogenous variables, see section Functions in Equation Definitions. Upper case function arguments indicate that a constant is required. Arguments in square brackets may be omitted: the default values used in such cases are specified in the function description provided.
The Fitpack Library
FITPACK by Paul Dierckx [42] is a Fortran-based library for one-dimensional and two-dimensional spline interpolations. This library has been repackaged to work with the GAMS Function Library Facility. The model [FITLIB01] from the GAMS Test Library is an example of the use of the library FITPACK inside GAMS.
Note that the supporting points to which the function will be fit need to be stored in a three-dimensional parameter, fitdata
, in a GDX file fit.gdx
. The first dimension is a function index, the second dimension is the index of the supporting point and the third dimension takes one of the following four values: "w"
(weight), "x"
(x-value), "y"
(y-value) or "z"
(z-value).
The FITPACK library is made available with the following directive:
$funcLibIn <InternalLibName> fitfclib
It provides the following functions:
Function | Description | End. Classif. |
---|---|---|
fitFunc(FUNIND,x[,y]) | Evaluate spline | DNLP |
fitParam(FUNIND,PARAM[,VALUE]) | Read or set parameters | none |
The function FitParam
may be used to change certain parameters that are used for the evaluation. The following values are defined:
- 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)
The Piecewise Polynomial Library
The Piecewise Polynomial Library may be used to evaluate piecewise polynomial functions. An example is given in the model [PWPLIB01] in the GAMS Test Library. Note that the functions that are to be evaluated need to be defined and stored in a GDX file. The following code snippet serves as illustration:
* 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, which will be read by external library (pwpcclib)
$gdxout pwp.gdx
$unload pwpdata
$gdxout
Observe that on each row of the table pwpdata
we have the following entries:
FuncInd.SegInd leftBound Coef0 Coef1 Coef2
Here FuncInd
sets a function index and SegInd
defines the index of the segment (or interval) which is described. Further, LeftBound
gives the lower bound of the segment. The upper bound will be taken from the lower bound on the following segment, or set to infinity in case it is the last segment. Finally, CoefX
defines the X
-th degree coefficient of the polynomial corresponding to this segment.
The Piecewise Polynomial Library is made available with the following directive:
$funcLibIn <InternalLibName> pwpcclib
It provides the following function:
Function | Description | End. Classif. |
---|---|---|
pwpFunc(FUNIND,x) | Piecewise polynomials | DNLP |
The Stochastic Library
The Stochastic Library provides random deviates, probability density functions, cumulative density functions and inverse cumulative density functions for certain continuous and discrete distributions. This library is made available with the following directive:
$funcLibIn <InternalLibName> stodclib
The continuous distributions that are available with this library are the following:
Distribution | 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 change 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) | Lognormal 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 , where 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 |
Further, the following discrete distributions are available:
Distribution | 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 FAILURES being the number of failures until the experiment is stopped and success probability P in each trial. The generated random number describes the number of successes until the defined number of failures is reached, see MathWorld |
poisson(LAMBDA) | Poisson distribution with mean LAMBDA , see MathWorld |
uniformInt(LOW,HIGH) | Integer Uniform distribution between LOW and HIGH , see MathWorld |
Note that for each distribution the library offers the following four functions, where DistributionName
is the name of the distribution as listed in the tables above, parameters
are the parameters associated with each distribution, and x
is the point at which the function is to be evaluated. Note that x
may be an endogenous variable.
For example, the functions for the Normal distribution are
Finally, the seed for the various random number generators can be set by using the following function:
Function | Description | End. Classif. |
---|---|---|
SetSeed(SEED) | Defines the seed for random number generator | none |
In the following example, a sample of size 20 is generated from the Normal, Binomial, Cauchy, and Lognormal distributions each:
$funcLibIn stolib stodclib
Functions randnorm /stolib.dnormal /
randbin /stolib.dbinomial /
randcauchy /stolib.dcauchy /
randlognorm /stolib.dlognormal /;
Set i / i1*i20 /;
Set j / norm, binomial, cauchy, lognorm /;
Parameter randx(i,j) "distribution sample";
randx(i,"norm") = randnorm(5,2);
randx(i,"binomial") = randbin(10,0.5);
randx(i,"cauchy") = randcauchy(5,1);
randx(i,"lognorm") = randlognorm(1.2,0.3);
display randx;
In the example, first the stochastic library is made available in GAMS, then the functions that will be used from the library are declared, giving them names under which to refer to them in the GAMS model.
The output generated by the display statement is the following:
---- 15 PARAMETER randx distribution sample norm binomial cauchy lognorm i1 4.373 4.000 5.520 3.132 i2 5.655 6.000 6.813 4.192 i3 5.927 6.000 5.426 2.801 i4 1.340 4.000 5.898 3.689 i5 3.537 3.000 -3.069 2.746 i6 3.057 5.000 5.518 3.430 i7 4.212 3.000 0.136 1.577 i8 6.869 7.000 5.068 3.857 i9 3.481 4.000 -13.856 3.977 i10 5.001 4.000 5.274 2.261 i11 3.182 5.000 4.383 2.686 i12 5.688 6.000 2.914 2.102 i13 3.675 6.000 4.693 3.105 i14 4.028 5.000 1.813 2.418 i15 8.767 5.000 12.190 2.182 i16 3.558 3.000 -112.644 2.092 i17 2.402 4.000 4.078 2.523 i18 2.249 2.000 0.996 2.777 i19 5.639 4.000 3.931 2.678 i20 7.374 4.000 4.671 3.159
The 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
Observe that a LINDO license is required to use this library.
The following table list the LINDO sampling functions.
Function | Description | End. Classif. |
---|---|---|
sampleLS<DistributionName>(parameters) | Creates a sample using the distribution DistributionName , according to the distribution parameters , and returns a HANDLE that references the sample, as illustrated in the example below. | none |
getSampleValues(HANDLE) | Retrieves sampling created by the function sampleLS. See example below. | none |
induceCorrelation(CORTYPE) | Induces the correlation that was set with the function setCorrelation before. CORTYPE describes the correlation type: 0 (Pearson), 1 (Kendall) or 2 (Spearman). See example below. | none |
setCorrelation(SAMPLE1,SAMPLE2,COR) | Defines correlation between two samplings. See example below. | none |
setSeed(SEED) | Specifies the seed for the random number generator. | none |
setRNG(RNG) | Specifies the random number generator that will be used. Possible values are -1 (FREE), 0 (SYSTEM), 1 (LINDO1), 2 (LINDO2), 3 (LIN1), 4 (MULT1), 5 (MULT2), and 6 (MERSENNE). | none |
The following tables list the available continuous and discrete distributions, respectively. Note that the parameter SAMSIZE
must be specified and describes the size of the sample. However, the parameter VARRED
is optional and facilitates choosing a variance reduction method. The values are 0
(meaning "none"), 1
(meaning "Latin Hyper Square") and 2
(meaning "Antithetic"). The default is Latin Hyper Square sampling, it will be used if no variance reduction method is specified.
Continuous Distribution | Description |
---|---|
beta(SHAPE_1,SHAPE_2,SAMSIZE[,VARRED]) | Beta distribution specified by two shape parameters. |
cauchy(LOCATION,SCALE,SAMSIZE[,VARRED]) | Cauchy distribution specified by the location and the 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 the function sampleLSf uses another version of the F distribution than the function dF from the Stochastic Library. |
gamma(SHAPE,SCALE,SAMSIZE[,VARRED]) | Gamma distribution specified by shape and scale parameter. Note that the function sampleLSgamma(A,B) is equivalent to the function 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 function sampleLSweibull(A,B) is equivalent to the function dWeibull(B,A) from the Stochastic Library. |
Discrete Distribution | 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 the function 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 the defined number of successes is reached. Note that the function sampleLSnegbinomial(R,P) is equivalent to the function 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 the effect of the functions setCorrelation
and induceCorrelation
:
$funcLibIn lsalib lsadclib
Functions normalSample / lsalib.SampleLSnormal /
getSampleVal / lsalib.getSampleValues /
setCor / lsalib.setCorrelation /
indCor / lsalib.induceCorrelation /;
Scalars d "dummy"
h "handle for first sample"
k "handle for second sample" ;
Set i "sample index" / i01*i12 /;
Parameters sv_h(i) "sample values for handle h"
sv_k(i) "sample values for handle k";
* generate two handles for 12 samples from normal distribution with mean 5 and std.dev. 2 each
h = normalSample(5,2,12);
k = normalSample(5,2,12);
* retrieve sample values from Lindo library
loop(i, sv_h(i) = getSampleVal(h) );
loop(i, sv_k(i) = getSampleVal(k) );
display sv_h, sv_k;
* set and induce a correlation between samples h and k
d = setCor(h,k,-1);
d = indCor(1);
* retrieve sample values again from correlated distribution
loop(i, sv_h(i) = getSampleVal(h) );
loop(i, sv_k(i) = getSampleVal(k) );
display sv_h, sv_k;
The resulting output shows that the values of sv_k
are reordered according to the desired correlation:
---- 25 PARAMETER sv_h sample values for handle h i01 2.079, i02 6.454, i03 4.437, i04 2.747, i05 5.339, i06 4.059, i07 6.311, i08 7.512, i09 8.280, i10 3.380, i11 4.596, i12 5.752 ---- 25 PARAMETER sv_k sample values for handle k i01 5.509, i02 3.021, i03 7.550, i04 6.002, i05 4.227, i06 0.704, i07 3.890, i08 9.474, i09 5.084, i10 4.592, i11 3.311, i12 6.442 ---- 35 PARAMETER sv_h sample values for handle h i01 2.079, i02 6.454, i03 4.437, i04 2.747, i05 5.339, i06 4.059, i07 6.311, i08 7.512, i09 8.280, i10 3.380, i11 4.596, i12 5.752 ---- 35 PARAMETER sv_k sample values for handle k i01 7.550, i02 3.021, i03 9.474, i04 6.442, i05 4.592, i06 0.704, i07 3.890, i08 6.002, i09 3.311, i10 5.084, i11 4.227, i12 5.509
The Mutex Library
The Mutex Library allows us to work with a process-level mutex - a program object that can be used to synchronize multiple processes. For example, a mutex can be used to coordinate mutually exclusive access to a resource shared between processes and, thus, prevent multiple processes from accessing the shared resource at the same time. A mutex object is created with a unique name and can be acquired or locked by only one process at a time. Each process waits to acquire the lock for (i.e. ownership of) the mutex object before executing the code that accesses the shared resource. After accessing the shared resource, the process releases or unlocks the mutex object and other processes waiting for the mutex can acquire the lock and continue execution.
The Mutex Library may be used to prevent simultaneous access to shared resources when running GAMS programs concurrently and the errors and hard-to-explain behavior that result. Such behavior is not unusual or surprising: GAMS even has some language features to spawn programs asynchronously (e.g. Asynchronous Execution), and these concurrent GAMS programs often share files. Consider the following example where we have two GAMS programs A and B concurrently running a common code that writes a GDX file containing a scalar x
:
scalar x;
repeat
execute_load 'x.gdx', x;
x = x + 1;
execute_unload 'x.gdx', x;
until x > 100;
To ensure correct execution of the code, we want to ensure that programs A and B have mutually exclusive access to the GDX file. In addition, we want the load-increment-unload sequence to occur atomically (i.e. no other process reads or writes the GDX file during this sequence). The extrinsic function library mtxcc
provides us with some functions to accomplish this and similar tasks:
Create(x)
: Creates a mutex with idx
in the system. Returns0
on success andUNDF
in case of error.Lock(x)
: Request to aquire or lock a mutex with idx
: it returns immediately if the mutex is available, otherwise it waits until the mutex becomes available. Returns0
on success andUNDF
in case of error. Any process waiting in this call will wait forever if the mutex is never unlocked.Unlock(x)
: Unlocks or releases a previously locked mutex with idx
. Returns0
on success andUNDF
in case of error. If processes are waiting to acquire a lock onx
, one of them will acquire the lock. Note that the lock should be held by the unlocking process: behavior is undefined otherwise.TryLock(x)
: LikeLock
, but if the mutex is not available, return1
immediately instead of waiting.TimedLock(x,y)
: LikeLock
, but instead of waiting forever, return1
after waiting unsuccessfully fory
milliseconds.Delete(x)
: Erases a mutex with idx
from the system. Returns0
on success andUNDF
in case of error. Note that a process does not need to own the mutex to erase it. If a mutex is deleted and there are remaining processes queued up waiting for the unlock, these processes will wait forever.
Continuing our previous example, we need to make the functions Lock
and Unlock
available so we can use them to protect the entire loop body by enclosing it in Lock
/Unlock
calls.
$funcLibIn mtxlib mtxcclib
Function lock / mtxlib.Lock /
unlock / mtxlib.Unlock /;
scalar x;
repeat
abort$lock(12345) 'Problems locking mutex';
execute_load 'x.gdx', x;
x = x + 1;
execute_unload 'x.gdx', x;
abort$unlock(12345) 'Problems unlocking mutex';
until x > 100;
Before we can use the mutex, we first have to create it. Also if both programs are done, we need to make sure that the named mutex is erased at the end of the programs, otherwise named mutexes collect in your system and consume resources (memory or disk depending on the implementation). Because the deletion of a named mutex might be forgotten or a program exits while holding the lock to the mutex or another concurrent system uses the same id, it is a good idea to find a unique mutex name, create it once, and erase it at the end as shown in the following complete example:
* generate a unique mutex ID
$eval mtxID round(frac(jnow)*24*60*60*1000)
$onEcho > x.gms
$funcLibIn mtxlib mtxcclib
Function lock / mtxlib.Lock /
unlock / mtxlib.Unlock /;
scalar x;
repeat
abort$lock(%mtxID%) 'Problems locking mutex';
execute_load 'x.gdx', x;
x = x + 1;
put_utility 'log' / 'Program %prgID% increased x to ' x:0:0;
execute_unload 'x.gdx', x;
abort$unlock(%mtxID%) 'Problems unlocking mutex';
until x >= 100;
$offEcho
* Initialize x.gdx
scalar x /1/;
$gdxOut x.gdx
$unload x
$gdxOut
$onEcho > create.gms
$funcLibIn mtxlib mtxcclib
Function create / mtxlib.Create /;
abort$create(%mtxID%) 'problems creating mutex';
$offEcho
* Create the mutex
$call.checkErrorLevel gams create.gms lo=2
* Asynchronously spawn A and B
$call.async 'gams x.gms --prgID=A fileStem=A'
$eval jhA JobHandle
$call.async 'gams x.gms --prgID=B fileStem=B'
$eval jhB JobHandle
* Wait until both jobs are back
$set statA 0
$set statB 0
$label l1
$eval x sleep(1)
$if not %statA% == 2 $eval statA JobStatus(%jhA%)
$if not %statB% == 2 $eval statB JobStatus(%jhB%)
$if not %statA%%statB% == 22 $goto l1
$onEcho > erase.gms
$funcLibIn mtxlib mtxcclib
Function delete / mtxlib.Delete /;
abort$delete(%mtxID%) 'problems erasing mutex';
$offEcho
* Erase the mutex
$call.checkErrorLevel gams erase.gms lo=2
The log indicates which program updated x
and the GDX file x.gdx
at what stage:
Program B increased x to 20 Program B increased x to 21 Program A increased x to 22 Program A increased x to 23
Models dicegrid and asyncincbi show more realistic uses for the mtxcc
library.
Build Your Own Library
This section discusses the creation of a custom extrinsic function library. Before attempting to implement such a library, we suggest to study the example libraries for which source code and test models are available. These libraries are studied below.
- Attention
- Building extrinsic function libraries requires the knowledge of a regular programming language (like C/C++, Fortran, ...) and experience with handling compilers and linkers to build dynamically linked libraries. In some situation, it may be easier to use the simpler GAMS macros or batinclude files to define own functions.
- Note
- Extrinsic functions are limited to 20 scalar arguments and return a scalar value.
An extrinsic function library consists of a specification part and a number of callbacks to evaluate the defined functions at an input point.
The specification part is implemented by a callback querylibrary
. It returns information about the library itself, available functions, their arguments, endogenous classification, etc. to the GAMS execution system. C, Delphi, or Fortran source code for this callback can be generated automatically by using the Python helper script ql.py
. The script processes a specification file *.spec
, which is specified as first argument. The format of this file is documented in the file tri.spec
. Both ql.py
and tri.spec
are contained in the source of the trigonometric library examples in the GAMS test library, obtainable via
$ testlib trilib01 $ unzip trisource.zip
If an extrinsic function will be used within equations of a GAMS model, next to the function value evaluation callback, also callbacks that compute first and second derivatives with respect to all endogenous arguments at an input point should be provided. Occasionally, this can be inconvenient. Observe that GAMS can use the function values at points close to the input point to estimate the derivate values using finite differences. However, this method is not as accurate as analytic derivatives and requires a number of function evaluations, thus 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. However, a better alternative is often the use of automatic differentiation when implementing the function evaluation. This is demonstrated in the CPP Library Example, see section Automatic Differentiation for more details.
GAMS offers some support to check the implementation of of derivatives for extrinsic functions via the function suffixes grad
, gradn
, hess
and hessn
. These function suffixes are defined for intrinsic and extrinsic functions. For example, for an extrinsic function userfunc
, the gradient evaluation that the user implemented may be called with userfunc.grad
. Further, an approximation of the gradient by finite differences is available by calling userfunc.gradn
. Comparing the results of these two calls can often help to check the implementation of the gradient. The same principle applies for the Hessian and the function suffixes .hess
and .hessn
. The GAMS options FDDelta and FDOpt can be used to influence the finite difference calculations. For more details, see model [DERIVTST] in the GAMS Model Library.
Example: Trigonometric Library
The Trigonometric 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. In addition, the source code in C, Delphi, and Fortran is available in the include files of the models [TRILIB01], [TRILIB02] and [TRILIB03] respectively. The library implements the following extrinsic functions:
Function | Description | End. Classif. |
---|---|---|
setMode(MODE) | Sets mode globally. Possible values are 0 for radian and 1 for degree. May be overwritten by the optional argument MODE in the functions cosine and sine . | none |
cosine(x[,MODE]) | Returns the cosine of the argument x . Note that the argument MODE is optional, default setting: MODE = 0 . | NLP |
sine(x[,MODE]) | Returns the sine of the argument x . Note that the argument MODE is optional, 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 be implemented by a GAMS extrinsic function library. The file tricclibql.c
implements the querylibrary
callback, which provides information about the library itself and the extrinsic functions it implements. For example, the information that the function cosine
has an endogenous required first argument and an exogenous optional second argument is available from the querylibrary
callback. The file tricclibql.c
(and also the Delphi and Fortran90 equivalents tridclibql.inc
and triifortlibql.f90
) has been generated by the script ql.py
by processing the specification file tri.spec
.
Example: The CPP Library
The CPP 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], [CPPLIB02], [CPPLIB03], [CPPLIB04], and [CPPLIB05] are more thorough tests for the CPP library extrinsics shipped with the distribution. These functions are listed and described in the following table. Note that in keeping with the language conventions of statistics, PDF is shorthand for "probability density function" and CDF is shorthand for "cumulative distribution function".
Function | Description | End. Classif. |
---|---|---|
pdfUVN(x) | PDF of uni-variate Normal distribution, see MathWorld or R | NLP |
cdfUVN(x) | CDF of uni-variate Normal distribution, see MathWorld or R | NLP |
pdfBVN(x,y,r) | PDF of bivariate Normal distribution, see MathWorld or R | NLP |
cdfBVN(x,y,r) | CDF of bivariate Normal distribution, see MathWorld or R | NLP |
pdfTVN(x,y,z,r21,r31,r32) | PDF of trivariate Normal distribution, see MathWorldor R | NLP |
Automatic Differentiation
Often, extrinsic functions are created in order to be used with endogenous arguments. In such cases it is necessary to provide first and second derivatives with respect to these arguments in addition to the function values themselves. One way to compute these derivatives is via automatic differentiation techniques (see the article in Wikipedia for details).
Note that with C++ it is possible to overload the usual arithmetic operators (assignment, addition, multiplication, etc.) so that 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. Observe that the Test Library model [CPPLIB00] includes all the source code for the CPP Library and illustrates the steps needed 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
As shown in the table above, 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 functions pdfNormal and cdfNormal from the stochastic library. For the multivariate cases, the implementation is based on TVPACK from Alan Genz, with some modifications to allow for proper computation of derivatives. Note that we chose to implement the functions taking correlation coefficients as arguments, not a covariance matrix. The conversion from a covariance matrix to correlation coefficients is straightforward. The following R
code describes this conversion. Here the package mnormt
is used that computes the multivariate CDF with similar code to Genz:
Note that for the bivariate case there is one correlation coefficient r
. The CDF implementation is not quite accurate to machine precision, it has 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 lower as the condition number of R
increases. Typically, an accuracy of \(10^{-11}\) is all that can be achieved. In both multivariate cases, we recommend to 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. As the correlation matrix becomes degenerate, the distribution becomes degenerate too.
Remark: Stateful Libraries
While GAMS intrinsic functions are stateless, users may implement stateful
extrinsic functions, i.e., functions that have some memory. There are two ways to achieve stateful extrinsic functions:
- Library initialization (
LibInit
): at initialization time, the function library reads some data that is required to evaluate the provided functions. An example is the Piecewise Polynomial Library. - Previous function calls: function calls that alter the execution of successive function calls. An example is the function
SetMode
from the Trigonometric Library.
- Attention
- Altering the state of an extrinsic function library by function evaluation calls is problematic, since different parts of the GAMS system potentially use different instances of the function library.
For example, consider that the function SetMode
of the Trigonometric Library is called via SetMode(1)
before a solve statement. Unless option solvelink is set to 5, the solver will run in a separate process with a new instance of the function library and therefore will use the default mode
, which is zero. Further, if solvelink
is set to zero, the GAMS process will terminate in order to execute the solve statement and will restart a new GAMS process after the solve. The restarted GAMS process will load a fresh instance of the extrinsic function library, which has no memory of the value of mode
from before the solve statement. This problem is demonstrated in the GAMS Test Library model [TRILIB04].
Extrinsic Functions vs. External Equations
In addition to extrinsic functions, GAMS offers another facility to include additional mathematical functions in GAMS: external equations. These equations are denoted by equation type =x=
. A feasible solution for a model instance must satisfy all internal and external equations. External equations are introduced and discussed in chapter External Equations.
Similar to extrinsic functions, it is the users responsibility to provide routines that evaluates the external equation. Further, both facilities are especially pertinent to nonlinear models.
An overview of some characteristics of extrinsic functions and external equations is given in the following table:
Characteristic | Extrinsic Function | External Equation |
---|---|---|
Maximum number of arguments | 20 | No limit |
Available in statements | Yes | No |
Debugging support | Yes | No |
Returns Hessian to solver | Yes | No |
Consider the following example, which demonstrates how an equation \( z = \int_{x_1}^{x_2} x\, dx \) may be formulated using either an extrinsic function and an external equation:
* using an extrinsic function
e1.. integralx(x1,x2) =e= z;
* using an external equation
e1.. 1*x1 + 2*x2 + 3*z =x= 1;
Note that the function integralx
is assumed to be a user-defined extrinsic function, probably implemented as \( \texttt{integralx}(x_1,x_2) = \frac{1}{2}(x_2^2-x_1^2)\).
Note further, that the integers 1 to 3 in the external equation are mapped to indices for variables that are used in external functions inside an external module. The right-hand side of the external equation denotes the external equation number. It is assumed that the user-defined external equation implements the function \((x_1,x_2,z)\mapsto \int_{x_1}^{x_2} x\, dx - z\).
As a simple example for the use of an extrinsic function in a statement, consider the following assignment, which sets the level of variable y
to the value of \((\int_{0}^{2} x\, dx)^2\):
y.l = sqr(integralx(0,2));