[R] Architecting an optimization with external calls

[Ross Boylan]

On Tue, 2003-11-04 at 13:12, Prof Brian Ripley wrote:
> Look into external pointers.  That is how we have tackled this, e.g. in 
> the ts package.
I got the following to work.  Any comments?  I indicated some areas of
uncertainty in the comments.  I'm also unsure about the differences, if
any, between the DOUBLE_DATA(x)[0] = val style and SET_ELEMENT(x, i,
val) style.

The R code was the best way I could think of to retain the state while
using a function that worked with optim.  The C code remembers the
argument from the previous call and prints it out when invoked.

The driving R code was
loglika <- function(initialparms) {
  opaque <- NULL # holds the data returned by the C function
  innerloglik <- function(parms){
    # If .C can handle externalptr, it would be easier
    result <- .Call("loglik", as.double(parms), opaque)
    opaque <<- result$opaque
    return(result$loglik)
  }
  result <- optim(initialparms, innerloglik)
  .Call("cleanup", opaque)
  opaque <- NULL
}

And the C:
#include <R.h>
#include <Rdefines.h>

SEXP loglik(SEXP params, SEXP data){
  double *pd;
  SEXP returnLik, returnVal, returnNames;

  if (TYPEOF(data) == EXTPTRSXP) {
    pd = (double *) EXTPTR_PTR(data);
    Rprintf("I've been here with %f\n", *pd);
  } else {
    pd = malloc(sizeof(double));
    data = R_MakeExternalPtr(pd, R_NilValue, R_NilValue);
    /* do I need PROTECT(data)?  if so, argument to
     UNPROTECT is variable */
  }
  
  *pd = DOUBLE_DATA(params)[0];

  /* Real computation would go here */

  PROTECT(returnLik = NEW_NUMERIC(1));
  NUMERIC_POINTER(returnLik)[0] = (*pd)*(*pd); /* our log likelihood */
     /* When I just set it = 5.0 I got a segfault.  Maybe because I
	rebuilt the library at the same time.... */

  /* update opaque data with first input parameter */
  Rprintf("Just recorded input parameter %f\n", *pd);

  /* construct return list object */
  PROTECT(returnVal = NEW_LIST(2));
  PROTECT(returnNames = NEW_CHARACTER(2));
  SET_STRING_ELT(returnNames, 0, mkChar("loglik"));
  SET_STRING_ELT(returnNames, 1, mkChar("opaque"));
  SET_NAMES(returnVal, returnNames);
  SET_ELEMENT(returnVal, 0, returnLik);
  SET_ELEMENT(returnVal, 1, data);
  UNPROTECT(3);
  return returnVal;
}

SEXP cleanup(SEXP data){
  if (TYPEOF(data) == EXTPTRSXP) {
    free(EXTPTR_PTR(data));
  } else {
    error("cleanup called without an External Pointer argument");
  }
  return R_NilValue;
}

I was going on a lot of hunches, even with the docs and code together.