Dear list,
I'm looking for some advice on a specific problem. Using RcppDE there is
the possibility to give the optimizer directly an external pointer to
the C++ function it will use as the objective function. I found this
mechanism pretty useful as it may speed up things quite a lot (I have a
problem where the speedup is from 17 minutes to some seconds), so that I
use the same mechanism as RcppDE in our package Rmalschains and in the
Rdonlp2 package, which is available from Rmetrics on Rforge.
The problem that this mechanism has is that it cannot handle additional
parameters to the objective function. Having additional parameters is
often essential, because if you fit a model to data you need the data
available in the target function. I illustrate the problem with an
example I took from the RcppDE tests:
#-----------------------------------------
library(inline)
library(RcppDE)
inc <- 'double rastrigin(SEXP xs) { //here I want to give it an
additional parameter: SEXP additional_parameter
//Do something with the parameter, e.g. use it for result
calculation. Here we just want to print it
//double my_additional_parameter =
Rcpp::as<double>(additional_parameter);
//Rprintf("ap: %f\\n", my_additional_parameter);
Rcpp::NumericVector x(xs);
int n = x.size();
double sum = 20.0;
for (int i=0; i<n; i++) {
sum += x[i]+2 - 10*cos(2*M_PI*x[i]);
}
return(sum);
}
'
src.xptr <- '
typedef double (*funcPtr)(SEXP);
return(XPtr<funcPtr>(new funcPtr(&rastrigin)));
'
create_xptr <- cxxfunction(signature(), body=src.xptr, inc=inc,
plugin="Rcpp")
n <- 10
maxIt <- 100
res <- RcppDE::DEoptim(fn=create_xptr(), lower=rep(-25, n),
upper=rep(25, n),
control=list(NP=10*n, itermax=maxIt, trace=FALSE)) #,
additional_paramater=25)
res$optim
#-----------------------------------------
I currently get around this by having a global singleton object which
holds these parameters. This works but of course is not very nice when
it comes to parallelization. The code is more or less like this:
//----------------------------------------------
class TargetFunction {
private:
static TargetFunction *TargetFunctionSingleton;
std::vector<double> param;
double objval;
public:
void eval(const double* x, int n) {
double sum = 20.0;
for (int i=0; i<n; i++) {
sum += x[i]+2 - 10*cos(2*M_PI*x[i]);
};
//here I can use the parameter now!!
Rprintf("ap: %f\\n", param[0]);
this->objval = sum;
};
void init(std::vector<double> & p_param) {
this->param = p_param;
};
static TargetFunction* getTargetFunctionSingleton() {
if( TargetFunctionSingleton == 0 )
TargetFunctionSingleton = new TargetFunction();
return TargetFunctionSingleton;
};
static void deleteTargetFunctionSingleton(){
if( TargetFunctionSingleton == 0 ) return;
else {
delete TargetFunctionSingleton;
TargetFunctionSingleton = 0;
}
return;
};
double getObjVal() {
return(objval);
};
};
TargetFunction* TargetFunction::TargetFunctionSingleton = 0;
RcppExport SEXP targetFunction(SEXP p_par)
{
Rcpp::NumericVector par(p_par);
TargetFunction* sp = TargetFunction::getTargetFunctionSingleton();
sp->eval(par.begin(), par.size());
return Rcpp::wrap(sp->getObjVal());
}
RcppExport SEXP targetFunctionInit(SEXP p_param) {
TargetFunction::deleteTargetFunctionSingleton();
TargetFunction* sp = TargetFunction::getTargetFunctionSingleton();
std::vector<double> param = Rcpp::as< std::vector<double> >(p_param);
sp->init(param);
return R_NilValue;
}
RcppExport SEXP GetTargetFunctionPtr() {
typedef SEXP (*funcPtr)(SEXP);
return (Rcpp::XPtr<funcPtr>(new funcPtr(&targetFunction)));
}
//-----------------------------------------------------
Now, before doing the optimization, I call targetFunctionInit and set
the additional parameters. Afterwards, everything is as in the example
above, and I have the additional parameters available in the target
function. Now the question is how I could solve this more elegantly, or
more R like. The first thing that comes to mind is to use an R
environment instead of the singleton. However, how can I do this? I
could have a singleton list of objects and then use the address of the R
environment as a hash to find the right object in the list. But this is
probably not really the way R environments should be used, and I wonder
if this will cause any trouble.
Any advise is highly appreciated.
Regards,
Christoph
--
Christoph Bergmeir
e-mail: c.bergm...@decsai.ugr.es
Grupo SCI2S, DiCITS Lab (http://sci2s.ugr.es/DiCITS)
Dpto. de Ciencias de la Computacion e Inteligencia Artificial
E.T.S. Ingenierias de Informatica y Telecomunicacion
Universidad de Granada
18071 - GRANADA (Spain)
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
