NAG Library Callback Functions in R

Introduction
Step 1: Wrapping the Library Routine
Step 2: Wrapping the Callback Function
Step 3: Compiling the Code for Use in R
Step 4: Writing R Code to call the New DLL
Amending the Example to Call the Fortran Library
A More Complex Example (E04UC)
Additional Information

Introduction

Many of the routines in the NAG Library involve callback functions. This documentation gives some notes on how such routines can be called from within R (www.r-project.org). This tutorial was written using R version 2.7.0 in the windows environment. The approach adopted should, however, be applicable to other platforms. Extensive use has been made of the document Writing R Extensions which is available from the R web site and contains additional information on many of the techniques used in these examples.

In order to demonstrate the techniques we consider the specific example of calling the E04 routine E04AB. E04AB searches for the minimum of a function in a given finite interval, of a continuous function of a single variable, using function values only. This routine is available in both the NAG C Libraries, as e04abc, and the NAG Fortran Library as E04ABF. Initially we concentrate on calling e04abc from the C Library, as this is slightly simpler. We later go on to explain the necessary changes to the code in order to call the corresponding Fortran routine.

There are four steps involved in calling a NAG Library routine, which uses a callback function, from within R:

Write a main wrapper function, in C, which:
- Takes R like objects as input.
- Converts them into C (or Fortran) variables.
- Calls the library routine.
- Takes the output from the library routine and converts it back into R like objects.
- Returns the output to the user.
Write a wrapper function, in C, for the callback function which:
- Takes C (or Fortran) variables from the NAG routine as input.
- Converts them into R like objects.
- Calls the callback function from within R.
- Takes the output from the callback function and converts it back into C (or Fortran) variables.
- Returns control to the NAG Library routine.
Compile the C code into a shared library (DLL).
Write some R code that loads the shared library and calls the C wrapper.

As can be seen from above, in order to make use of these instructions you need access to a C compiler irrespective of whether you are using the NAG C or Fortran Libraries. But, see step 3 below for some guidelines on using the mechanism built into the standard R installation, along with a set of freely available tools for this purpose.

Step 1: Wrapping the Library Routine

It is not possible to use any of the NAG Library routines that use callback functions directly from R. The first step in using these functions is to write a wrapper function in C. As this wrapper uses a large number of utility functions and types specific to R (all of which are supplied with the standard installation of R), this wrapper needs to be written in C, irrespective of whether you are eventually going to call the NAG C Library or the NAG Fortran Library.

Writing the wrapper can be broken down into a number of tasks:

Develop the interface for the wrapper.
Other than needing to return an object of type SEXP the interface is a matter of personal taste. In this example we are using the prototype:
```
SEXP we04abc(SEXP funct, SEXP ab,  SEXP e, SEXP mfun, SEXP rho)
```
where the R supplied macro SEXP is a pointer to a structure and indicates an R like object, which includes lists, vectors, scalars and functions, of any type.
The NAG Library routine being wrapped, e04abc, has 10 arguments; three input arguments (e1,e2 and max_fun), four input / output arguments (a, b, comm and fail), two output arguments (x and f) and a pointer to the function being minimized (funct). In prototype for we04abc we have combined e1 and e2 into an R vector of length 2, called e. The value of max_fun will be held in the R scalar mfun. The values of a and b, on entry, are held in the first and second elements of the R vector ab. We do not need equivalents to comm or fail. All of the various output values, that is x, f and a and b on exit from the NAG routine, will be returned in a list. That leaves the arguments funct and rho in the above prototype. The argument funct will hold the callback function and rho will hold an R environment in which funct will be evaluated.

There are two ways of supplying a callback function written in R to a C routine. You can either supply the body of the function, via the R command body, or the function itself, i.e:
- Method 1: As the body of the function:
```
.Call("we04abc",body(funct),as.double(bound),as.double(e),
       as.integer(mfun),new.env())
```
- Method 2: As the function itself:
```
 .Call("we04abc",funct,as.double(bound),as.double(e),
        as.integer(mfun),new.env())
```
where, for example:
```
funct = function(y) {
  sin(y) / y
}
```
The first method has the disadvantage that the names of the variables used in the body of the function must be known when writing the C code, i.e. if the C code assumes that the user will use y as in the above example, but instead they use sin(z) / z, then the routine will not work as expected (but see step 4 for a way around this). The second method has the disadvantage that the C code to implement it is more obscure. In this document we show examples of using both methods.
Convert the R objects into standard C variable types that can be understood by the NAG Library routine.
Prior to converting the R objects supplied by the user we check a couple of assumptions. This can be done using the R supplied functions LENGTH, isLanguage, isFunction and isEnvironment. Therefore to check that the vector ab has length of at least 2, use:
```
if (LENGTH(ab) != 2) {
  error("ab must be a numeric vector of length 2");
}
```
and to check that rho is an environment:
```
if (!isEnvironment(rho)) {
  error("rho must be an environment");
}
```
In the supplied example code we are allowing for both methods of supplying a callback function; as a function or as the body of the function:
```
if (!isLanguage(funct) && !isFunction(funct)) {
  error("funct must be a function or the body of a function");
}
```
The function error is part of the utilities supplied with R and reports errors to the user in an R like manner. Once the initial checks have been carried out, converting the R objects supplied by the user into the usual C variable types is a simple matter of calling some R supplied macros:
```
double a, b;
int max_fun;

a = REAL(ab)[0];
b = REAL(ab)[1];
max_fun = INTEGER(mfun)[0];
```
We do something similar for e. If e is not valid we set the accuracies to zero, which according to the documentation of e04abc will cause default values to be used:
```
e1 = REAL(e)[0];
if (!R_FINITE(e1)) e1 = 0.0;
if (LENGTH(e) > 1) {
  e2 = REAL(e)[1];
  if (!R_FINITE(e2)) e2 = 0.0;
} else {
  e2 = 0.0;
}
```
Where the macro R_FINITE checks for invalid values (including NA, NaN etc). Similarly we check for an invalid value for mfun and use a default value of 100 if one is given:
```
if (max_fun == NA_INTEGER) max_fun = 100;
```
Prepare any relevant user supplied information into a form that can be passed through the NAG Library routine to the callback function wrapper.
In this example, the only information that needs to be passed through the NAG Library routine to the callback function wrapper is the funct and rho arguments. Any additional data that the callback function may require can be passed as global variables from within R. Once the input parameters have been converted we need to prepare the information to pass onto the callback function C wrapper. The mechanism for doing this is via the Nag_Comm communication structure. When calling e04abc from R, the two pieces of information that are required is the (pointer to) the R function (held in the input argument funct) and the (pointer to) the environment in which funct will be evaluated (held in the input argument rho). In addition we will be using an error flag, called ifail, to indicate if the callback function returns an invalid value. We therefore need to allocate sufficient memory to hold these three pieces of information:
```
/* Allocate memory for the communication array, to hold f, rho and ifail */
if (!(p = (void **) malloc(3*sizeof(void *)))) {
  error("Memory allocation error");
}
```
cast the values to void * and store them:
```
/* Convert function pointer funct and rho into Nag_Comm structure */
p[0] = (void *) funct;
p[1] = (void *) rho;
p[2] = (void *) &ifail;
```
and then make sure that the p pointer in the Nag_Comm structure, comm, points to the correct place.
```
comm.p = p;
```
Call the NAG routine.
Which in this example, is just a case of:
```
e04abc(we04abcL, e1, e2, &a, &b, max_fun, &x, &f, &comm, &fail);
```
where we04abcL is a C wrapper to the callback function (see step 2).

Check the error returns from the NAG Library.

For example:

if (fail.code != NE_NOERROR) {
  error(fail.message);
}

Return the results from the NAG function to the user.
As with the prototype for the C wrapper, how the output is returned to the user is a personal choice. The library routine e04abc returns four pieces of information that may be of interest (excluding any possible error exits, which have been handled above), these are x, the estimated value of the minimum, f, the value of the callback function evaluated at x and a and b, the lower and upper bounds for the minimum respectively. In this example we have decided to return these values in a list, with the bounds, a and b, as a vector with two elements. This can be achieved using:
```
/* Get the list of results */
PROTECT(ans = allocVector(VECSXP,3));
/* Return the estimated point of the minimum (x) */
SET_VECTOR_ELT(ans,0,nag_d2r(1,&x));
/* Return the value of F at x (f) */
SET_VECTOR_ELT(ans,1,nag_d2r(1,&f));
/* Return the lower and upper bounds for x */
SET_VECTOR_ELT(ans,2,nag_2d2r(1,&a,1,&b));
```
The R supplied macro PROTECT, stops the garbage collector removing the object inside of it and the R supplied macro SET_VECTOR_ELT adds elements to a list. The routine nag_d2r is a utility function which takes a double precision vector and returns them in the form of a SEXP object and nag_2d2r takes two double precision vectors and appends them into an SEXP object. The code for these two utilities is included with the other example source for this document.
To aid the eventual user of our "e04abc" R routine we also name the elements of the list:
```
/* Get names for the list of results */
PROTECT(names = allocVector(STRSXP, 3));
SET_STRING_ELT(names,0,mkChar("x"));
SET_STRING_ELT(names,1,mkChar("f"));
SET_STRING_ELT(names,2,mkChar("ab"));
setAttrib(ans,R_NamesSymbol,names);
```
Once we have finished we need to unprotect the objects we protected:
```
UNPROTECT(2);
```
The R supplied macro UNPROTECT unprotects the last "n" variables protected. In our example, two variables are unprotected, ans and names. Returning the resulting list to the user is just a matter of returning the variable ans.

Step 2: Wrapping the Callback Function

As described in step 1 there are two ways that the callback function can be supplied to the C wrapper; as an R function, or as the body of an R function (via the R command body). Which of these approaches is adopted will affect how the callback function must be wrapped. Probably the simplest of the two is when the callback function is supplied via the R body command. We will consider this approach first. The function prototype for the callback function wrapper must match that described in the documentation for the NAG routine of interest, in this case:

void NAG_CALL we04abcL(double xc, double *fc, Nag_Comm *comm)

Note that the NAG_CALL at near the beginning of the prototype is required and is defined in the header files supplied with the NAG Library.

Writing the wrapper for the callback function can be broken down into a number of tasks:

Extract the pointers to the callback function and the R environment.

These pointers were packed into the Nag_Comm in step 1 and can be extracted using:

SEXP funct, rho;
void **p;
p = (void **) comm->p;

/* Get the pointer to the R function */
funct = (SEXP) p[0];
/* Get the pointer to the R environment */
rho = (SEXP) p[1];

Make the input arguments to we04abcL available to the callback function.
In this example the callback function only requires one argument, xc, the value at which the callback function is to be evaluated. To make the value of xc available to R, in a variable named x, in the environment defined by rho the following can be used:
```
defineVar(install("x"), nag_d2r(1,&xc), rho);
```
The function nag_d2r is as described in step 1 above. It should be noted that the name used in the install function, x in this case, must correspond to the name used in the body of the R code for the callback function.
Evaluate the callback supplied R function.
The NAG routine expects the value of the callback function, evaluated at xc to be stored in fc, therefore:
```
*fc = REAL(eval(funct,rho))[0];
```
The C function eval (supplied with the R headers etc) works similarly to the R function of the same name.
Check that the returned value is numeric.
When using a callback function written in C (or Fortran) we know that the type of the returned value will be numeric (as the compiler will not allow anything else). However, when the callback function is written in R the returned value could be anything; a character, a list etc. Prior to returning control to the NAG routine we need to check that the callback function has returned a numeric value. The advantage of doing this here, rather than in the calling R code, is that we can ensure that the C wrapper routine we04abc cleans up any allocated memory prior to terminating if an error occurs. Unfortunately, e04abc does not have a quick way of terminating if an error occurs in the callback function. In this case we therefore set a flag in the Nag_Comm argument indicating that an error has occurred and return a default numeric value.
```
if (!R_FINITE(*fc)) {
  ifail = (Integer *) p[2];
  *ifail = 1;
  *fc = 0.0;
}
```

If the second approach to handling the callback function is adopted, i.e. rather than the body of a function we expect the function itself to be supplied, steps (b) and (c), should be replaced by

SEXP R_fcall;

PROTECT(R_fcall = lang2(funct,nag_d2r(1,&xc)));
  
*fc = REAL(eval(R_fcall, rho))[0];
  
UNPROTECT(1);

A full copy of both versions of these wrappers is available in the example source.

Step 3: Compiling the Code for Use in R

The standard R installation comes with a mechanism for compiling C or Fortran code into a shareable library which can then be called by R. The mechanism is accessed via the R CMD SHLIB command line syntax. However, a number of tools, including a compiler, perl etc, common on other platforms is missing from windows. To remedy this, the R group have packaged together these missing tools into a package called Rtools, which is available for download. Its current location can be found from Appendix E of the R Installation and Administration document.

Once these extra tools have been installed the basic command for compiling a shared library (DLL), using the mechanism built into R, is

R CMD SHLIB we04abc.c nag_rwrapper_utils.c

this creates a DLL called we04abc.dll in the current directory. The resulting DLL exports all of the routines in both we04abc.c and nag_rwrapper_utils.c files. In order to use routines from the NAG Library you need to set the environment variable PKG_CFLAGS to point to the location of the NAG header files, for example:

set PKG_CFLAGS=-I. -Ic:/PROGRA~1/NAG/CL08/cldll084zl/include

and then use

R CMD SHLIB we04abc.c nag_rwrapper_utils.c 
   c:/PROGRA~1/NAG/CL08/cldll084zl/lib/CLDLL084Z_nag.lib

to link to the NAG Library at compilation time (assuming that you are linking to NAG Library CLDLL084Z and the installation directory was c:/PROGRA~1/NAG/CL08/cldll084zl/lib/).

If working from a unix like environment, for example via cygwin or a similar emulator, the command to set the environment variable becomes

setenv PKG_CFLAGS "-I. -Ic:/PROGRA~1/NAG/CL08/cldll084zl/include"

The rest of the instructions remain the same.

Step 4: Writing R Code to Call the New DLL

Prior to using our new routine, we must write an R front end for it. In the case where we are using the second method of passing on the callback function (i.e. passing the function directly), this is relatively straight forward:

we04abc = function(funct, bound, e=NA, mfun=NA) {
  ## R code to call e04abc: Calling the C code using the function f directly
  if (!is.loaded("we04abc")) {
    stop("External routine we04abc is not loaded")
  }
  
  ## Call NAG routine
  .Call("we04abc",funct,as.double(bound),as.double(e),as.integer(mfun),
        new.env())
}

It should be noted that, unlike in other tutorials on this web site, rather than using the built in R function .C to call a C routine, we are using .Call, this is because our C wrapper we04abc takes as input and returns R objects.

Finally the function that will be minimised, the callback function, can be written directly in R, for example:

funct = function(y) {
  sin(y) / y
}

Now all we do is load the DLL.

dyn.load("D:/sample_dll/we04abc.dll")

(where "D:/sample_dll" is the location of the DLL) and call the function from R in the usual way:

we04abc(funct, c(3.5,5.0))

It should be noted that a warning of the type:

Warning message:
In inDL(x, as.logical(local), as.logical(now), ...) :
  DLL attempted to change FPU control word from 8001f to 9001f

may appear when the DLL is loaded via dyn.load. This is a common issue when calling DLLs compiled using a visual studio compiler (as per the NAG Library) into R and can be safely ignored.

If we had decided to use the first method of handling the callback function as described in step 1, that is, we assumed that the funct argument to our we04abc C routine was going to be the body of a function, rather than the function itself, we would need to amend our R function as follows:

we04abc2 = function(funct, bound, e=NA, mfun=NA) {
  ## R code to call e04abc: Calling the C code using the body of a function
  ## rather than the function itself

  if (!is.loaded("we04abc")) {
    stop("External routine we04abc is not loaded")
  }
  
  i.funct = function(x) {
    ## Wrap the callback function
    x = funct(x)

    as.double(x)
  }

  ## Call NAG routine
  .Call("we04abc",body(i.funct),as.double(bound),as.double(e),
         as.integer(mfun),new.env())
}

this could then be called as before. The C code in step 2 was written with the assumption that the callback function would use a variable called x to hold the value at which the function is to be evaluated. The advantage of defining the function i.funct in the body of we04abc2, rather than using body(funct) directly in the .Call statement, is that by using x in i.funct and invoking the callback function with funct(x), we ensure that the correct variable names is used whilst allowing the user to use any variable name they want in funct.

Amending the Example to Call the Fortran Library

Calling a routine from the NAG Fortran Library is very similar to calling one from the C Library. The main differences are described below:

There is no Nag_Comm communication structure in the Fortran routine.
Rather than using a structure to pass information from the calling program, through the library routine to the callback function, many of the Fortran routines use an integer array, often called iuser and a double precision array, often called ruser. In order to store the pointers to the callback function funct and the environment rho in these arrays changes must be made in both the main C wrapper (as described in step 1) and the wrapper for the callback function (as described in step 2).
1. In the main C wrapper, replace the C code:
```
/* Allocate memory for the communication array, to hold f, rho 
   and ifail */
if (!(p = (void **) malloc(3*sizeof(void *)))) {
  error("Memory allocation error");
}

/* Convert function pointer funct and rho into Nag_Comm structure */
p[0] = (void *) funct;
p[1] = (void *) rho;
p[2] = (void *) &ifail;
```
  with:
```
long iuser[5];

/* Convert function pointer funct and rho into a couple of integers 
   and store in IUSER */
nag_convert_p2i(funct, &iuser[0]);
nag_convert_p2i(rho, &iuser[2]);
iuser[4] = 0;
```
  where the function nag_convert_p2i takes an object of type SEXP (which is a pointer to a structure) and saves it as two long integers. So in the above example, the value of funct is store in the first two elements of iuser, the value of rho in the third and fourth elements and the fifth element is used to store an error flag.
2. In the wrapper for the callback function, replace the C code:
```
void **p;
p = (void **) comm->p;

/* Get the pointer to the R function */
funct = (SEXP) p[0];
/* Get the pointer to the R environment */
rho = (SEXP) p[1];
```
  with:
```
/* Get the pointer to the R function */
nag_convert_i2p(&iuser[0], &funct);
/* Get the pointer to the R environment */
nag_convert_i2p(&iuser[2], &rho);
```
  where the function nag_convert_i2p is the inverse of nag_convert_p2i, converting two long integers into a pointer.
The communication array ruser is not used in this example.
Different routine name.
The name of the Fortran routine differs to the equivalent C routine name. Often this is just a case of changing the C name to upper case and swapping the last letter from a "c" to an "F". In the case of e04abc, there are two Fortran versions, E04ABF and E04ABA. The second of these, E04ABA, is a thread safe version which passes the two communication arrays iuser and ruser to the callback function. The first, E04ABF, assumes that any data required by the callback function is store in common blocks. As we make use of iuser in this example, we replace a call to e04abc:
```
e04abc(we04abcL, e1, e2, &a, &b, max_fun, &x, 
       &f, &comm, &fail);
```
with one to E04ABA:
```
E04ABA(we04abfL, &e1, &e2, &a, &b, &maxcal, 
       &x, &f, iuser, ruser, &ifail);
```
The is no NagError structure.
Rather than use a structure to return information about possible errors, the NAG Fortran Library uses an integer scalar argument, usually called ifail. On entry, the value of ifail indicates the behaviour of the routine if it detects and error. On exit, the value of ifail indicates what error, if any occurred. When calling a NAG Fortran routine, the parameter ifail should be set to either 1 or -1. A value of 0 on entry will cause the NAG routine to terminate R if an error is encountered. In order to report any errors, we must replace the following C code from step 1:
```
if (fail.code != NE_NOERROR) {
  error(fail.message);
}
```
with:
```
char line[100];
if (ifail != 0) {
  sprintf(line,
          "Error in NAG Library routine E04ABA, IFAIL on exit = %ld\n",
           (long) ifail);
  error(line);
}
```

Other than the points made above, calling a routine from the NAG Fortran Library is pretty much the same as calling one from the C Library. The example source contains examples of doing both, along with the utility functions nag_convert_p2i and nag_convert_i2p.

Additional Information

Information on NAG can be found at our web site, https://support.nag.com/.
Documentation on the NAG Fortran Library is available in a number of formats from here.
Documentation on the NAG C Library is also available in a number of formats from here.
All distributions of the NAG Fortran Library come with a C header file specific to the implementation, which helps you to call the library from C, as well as a document named techdoc.html which describes the use of the header file and the issues involved in mixed language programming. See your installed copy of the NAG Fortran Library for details. In addition, see
- NAG C Header Files FAQ for answers to some common questions.
- Additional resources on calling the NAG Library from R include:
  - Calling NAG Fortran Routines from R
  - NAG Library Callback Functions in R (A More Complex Example)
- Information on R can be found at their website www.r-project.org.
- Additional information on calling C routines into R can be found in the document Writing R Extensions