[Rcpp-devel] Best Practice for Optimize Functions in RCPP

[Hao Ye]

Simon,

There are also some minima-finding functions in GSL that you may want to look into. The source for RcppGSL might help with a fully c++ version.

Best,
--
Hao Ye
hye at ucsd.edu

On Dec 23, 2014, at 12:04 AM, Simon Riddell <simonruw at gmail.com> wrote:

> Hello,
> 
> I have been judiciously using RCPP for six months, and have had my numerous questions answered so far from previous threads. I have now run into an issue I cannot get a handle on:
> 
> I have converted a fairly large (recursive) Kalman filter function from R to C++ using RCPP. This function is then called within the R optim() function, subject to constraints, which estimates a set of ~30 parameters. 
> 
> My issue occurs within the C++ Kalman Filter function, where, I use the optimize() function to calculate the yield to maturity rate of various bonds. I do this by calling back to the R environment to access the optimize() function. Below is the R Function I create to be used within the Kalman filter, and below this R function is my method for calling it within the C++ code. To complicate matters further, the R Function calls a C++ Function. To clarify: The Kalman Filter C++ code calls an R function, and this R Function calls an additional separate C++ function. (Code included below)
> 
> As I iterate the Kalman filter it runs perfectly for anywhere from 30 minutes to six hours (and produces correct output when matched to the R code), but it inevitably crashes. From past reading I have done, Dirk has before mentioned that calling an R function many times within C++ is usually not a good idea. As a result I suspect this is my issue (the error codes vary, sometimes mentioning numerical errors, other times recursive errors, other times random RCPP error codes -- I can document and provide them if needed)
> 
> My biggest impasse is I cannot figure out a way to complete this without calling the R optimize() function, and cannot find any RCPP optimize functions to use instead. Thank you for reading.
> 
> 
> R Function (IRNPV.CPP is the C++ function, which R optimizes the rate parameter over):
> 
> optimcpp<-function(CoupDate=1,coupNo=1,coup=1,price=1,rate=1)
> {
> m<-optimize(function(CoupDate,coupNo,coup,price,rate) IRNPV.CPP(CoupDate=CoupDate,coupNo=coupNo,coup=coup,rate)-price)^2,c(-0.05,0.2),tol=1e-20,CoupDate=CoupDate,coupNo=coupNo,coup=coup,price=price)
> m$minimum
> }
> 
> Accessing the R environment within the C++ code:
> CPP.SRC <- '
> using namespace arma;
> Rcpp::Environment global = Rcpp::Environment::global_env();
> Function optimizecpp = global["optimcpp"];
> //Various matrix computations to calculate CD, CN, Co, & Pricex[0]
> optimvec0 = optimizecpp(CD,CN,Co,pricex[0],Ra);
> '
> 
> IRNPV.CPP Function (What the R Optimize() function optimizes 'rate' over -- Very Likely Unnecessary for Purposes of this Question)
> 
> IR.NPV.TIPS.CBF.SRC <- '
> using namespace arma;
> 
> double CD = Rcpp::as<double>(CoupDate);
> double CN = Rcpp::as<double>(coupNo);
> double RN = Rcpp::as<double>(rate);
> double Co = Rcpp::as<double>(coup);
> 
> Rcpp::NumericVector LM;
> mat DiscountFunc;
> double length;
> 
> double price;
> 
> if (CN > 1) {
> 
> length = floor((((CD+0.5*(CN-1))-CD)/0.50))+1;
> 
> for (double i = CD; i <= (CD+0.5*(CN-1))+.05; i += 0.5) {
> 	LM.insert(LM.end(),i);
> }
> 
> DiscountFunc.set_size(LM.size(), 1);
> DiscountFunc.fill(0);
> 
> 
> double k = 0;
> mat::row_iterator q = DiscountFunc.begin_row(0);  
> mat::row_iterator w = DiscountFunc.end_row((LM.size()-1));   
> for(mat::row_iterator i=q; i!=w; ++i)
>      {
>      (*i) = exp(-RN*LM[k]);
>      k = k + 1;
>      }
> 
> price = CD*Co*DiscountFunc[0];
> 
> for (int i=1; i<(LM.size()); ++i) {
> 	price = price+0.5*Co*DiscountFunc[i];
> 	}
> 
> 
> }
> else {
> double DiscountFunc;
> DiscountFunc = exp(-RN*CD);
> price = (1+CD*Co)*DiscountFunc;
> }
> return Rcpp::wrap(price);
> 
> '
> 
> 
> IRNPV.CPP <- cxxfunction(signature(CoupDate="NumericVector", coupNo="NumericVector", coup="NumericVector",rate="NumericVector"),
> 	IR.NPV.TIPS.CBF.SRC, include=CMATH, plugin="RcppArmadillo")
> 
> 
> 
> Thank you,
> Simon
> 
> 
> -- 
> Simon Alexander Riddell
> Economic Research RA
> Federal Reserve Bank
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141223/8f4d6f61/attachment.html>