Search code examples
c++rrcpp

Using R functions within a C++ function


I'm trying to convert R code implementing the golden section method to C++ code. Here is the R code:

goldensectionR <- function(f, dXl, dXr, dXm, dTol = 1e-9, ...) {
  
  dFr = f(dXr, ...)
  dFl = f(dXl, ...)
  dFm = f(dXm, ...)
  
  dRho = (1.0 + sqrt(5))/2.0
  
  if (dFl > dFm | dFr > dFm) {
    stop("Inital conditions are not satisfied")  
  }
  
  while (abs(dXr - dXl) > dTol) {
    
    if (dXr - dXm > dXm - dXl) {
      dY = dXm + (dXr - dXm)/(1.0 + dRho)
      dFy = f(dY, ...)
      if (dFy >= dFm) {
        dXl = dXm
        dXm = dY
      } else {
        dXr = dY  
      }
    }  else {
      dY = dXm - (dXm - dXl)/(1.0 + dRho)
      dFy = f(dY, ...)
      if (dFy >= dFm) {
        dXr = dXm
        dXm =dY
      } else {
        dXl = dY  
      }
    }   
    dFm = f(dXm, ...)  
  }
  return(dXm)
  
}

Here is my attempt so far to recreate this code in C++:

//[[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <Rcpp.h>

using namespace arma;
using namespace Rcpp;

// [[Rcpp::export]]
double goldensection(Function f, double dXl, double dXr, double dXm, double dTol = 1e-9) {
  
  double dFr = f(dXr);
  double dFl = f(dXl);
  double dFm = f(dXm);
  
  double dRho = (1.0 + sqrt(5))/2.0;
  
  if (dFl > dFm || dFr > dFm) {
    stop("Inital conditions are not satisfied");
  }
  
  while (fabs(dXr - dXl) > dTol) {
    
    if (dXr - dXm > dXm - dXl) {
      double dY = dXm + (dXr - dXm)/(1.0 + dRho);
      double dFy =  f(dY);
      if (dFy >= dFm) {
        dXl = dXm;
        dXm = dY;
      } else {
        dXr = dY; 
      }
    }  else {
      double dY = dXm - (dXm - dXl)/(1.0 + dRho);
      double dFy =  f(dY);
      if (dFy >= dFm) {
        dXr = dXm;
        dXm =dY;
      } else {
        dXl = dY;  
      }
    }   
    dFm = f(dXm);  
  }
  return(dXm);
    
}

The f off course refers to a function defined in R. When trying to import this code into R using sourceCpp I get in the places where I have used f in the function the error message:

cannot convert 'SEXP' {aka 'SEXPREC*'} to 'double' in initialization

So I'm obviously not using the R function correctly in the C++ function. How do you do this correctly?


Solution

  • The problem is that Rcpp doesn't know that the function will return a double, so what it makes it return instead is a SEXP - which is basically a wrapper that could stand in for lists, numbers, strings or other things. If you're sure your SExpr will always be a real you can use the asReal function to cast the SEXP to double, e.g.

    double dFy = asReal(f(dY));