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