### 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 [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.

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 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:

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));
```