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

jashander at ucdavis.edu jashander at ucdavis.edu
Tue Dec 23 08:14:44 CET 2014


Simon, 


If you're using Boost, // [[Rcpp::depends(BH)]] might be very helpful (if tools/minima.hpp is in the BH package). See http://gallery.rcpp.org/articles/a-first-boost-example/






 - Jaime

On Mon, Dec 22, 2014 at 10:37 PM, Simon Riddell <simonruw at gmail.com>
wrote:

> Brief Update,
> Avraham's advice helped me get a better idea of where to start. What I am
> now trying to do is learn from this post, which explains how to use
> external libraries in RCPP (
> http://stackoverflow.com/questions/13995266/using-3rd-party-header-files-with-rcpp).
> And secondly, try to use that to implement this optimize function (
> http://www.boost.org/doc/libs/1_57_0/libs/math/doc/html/math_toolkit/internals1/minima.html
> ).
> If anyone has additional links or input I would be grateful to receive it,
> I do not expect anyone to hold my hand through this.
> Thank you again for your time and help,
> Simon
> On Mon, Dec 22, 2014 at 9:25 PM, Avraham Adler <avraham.adler at gmail.com>
> wrote:
>> Hello, Simon.
>>
>> I ran into a similar problem before (with uniroot
>> <https://stat.ethz.ch/R-manual/R-devel/library/stats/html/uniroot.html>),
>> and what I found was that the function in R is actually written in C. I
>> found the original code, but not being a real programmer, I stopped working
>> on what I was doing and put it off until the point where I have enough time
>> to develop the skill to port that into RCpp (maybe 2047?). In your case, optimize()
>> <https://stat.ethz.ch/R-manual/R-devel/library/stats/html/optimize.html>itself
>> is a translation of an algorithm by Brent in FORTRAN which can be found here
>> (textfile) <http://www.netlib.org/fmm/fmin.f>. Perhaps you can translate
>> it into C++ and call it directly?
>>
>> Avi
>>
>> On Tue, 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
>>>
>>
>>
> -- 
> Simon Alexander Riddell
> London School of Economics
> linkedin.com/in/simonriddell <http://uk.linkedin.com/in/simonriddell>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141222/0550f93b/attachment.html>


More information about the Rcpp-devel mailing list