Extrinsic Functions

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 [43] 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.

Function Description End. Classif. for Continuous Distributions End. Classif. for Discrete Distributions
d<DistributionName>(parameters) Generate a random deviate (sample from the distribution) none none
pdf<DistributionName>(x,parameters) Probability density function DNLP none
cdf<DistributionName>(x,parameters) Cumulative distribution function DNLP none
icdf<DistributionName>(x,parameters) Inverse cumulative distribution function DNLP none

For example, the functions for the Normal distribution are

Function Description End. Classif.
dNormal(MEAN,STD_DEV) Samples a random number from the Normal distribution none
pdfNormal(x,MEAN,STD_DEV) Probability density function for Normal distribution DNLP
cdfNormal(x,MEAN,STD_DEV) Cumulative distribution function for Normal distribution DNLP
icdfNormal(x,MEAN,STD_DEV) Inverse cumulative distribution function for Normal distribution DNLP

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 id x in the system. Returns 0 on success and UNDF in case of error.
  • Lock(x): Request to aquire or lock a mutex with id x: it returns immediately if the mutex is available, otherwise it waits until the mutex becomes available. Returns 0 on success and UNDF 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 id x. Returns 0 on success and UNDF in case of error. If processes are waiting to acquire a lock on x, one of them will acquire the lock. Note that the lock should be held by the unlocking process: behavior is undefined otherwise.
  • TryLock(x): Like Lock, but if the mutex is not available, return 1 immediately instead of waiting.
  • TimedLock(x,y): Like Lock, but instead of waiting forever, return 1 after waiting unsuccessfully for y milliseconds.
  • Delete(x): Erases a mutex with id x from the system. Returns 0 on success and UNDF 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:

# 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)
GamsInteractive c
GamsWorkspace x

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:

  1. 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.
  2. 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));