Nice, thanks for sharing! FWIW, I wrote a wrapper for cubature to work at the C++ level with armadillo objects, but cubature is not particularly efficient in 1D. It has some nice features though, such as vector-valued integrals.
https://github.com/baptiste/cubature/blob/master/minimal.c (probably very poorly implemented, I was very surprised when I got it to work!) I use it in my planar package to do a 2D integration of a 3D complex vector field (thus returning 6 real components), https://github.com/baptiste/planar/blob/master/src/gaussian_beam.cpp#L478 Cheers, b. On 30 July 2014 22:25, Yixuan Qiu <[email protected]> wrote: > Hi Dave, > I wrote a small wrapper Rcpp_integrate(), which simplifies the interface > to the R internal function Rdqags. Basically you need to define your own > integrand and parameters passed to that function. Thanks to the default > arguments introduced in C++, in most cases you only need to specify four > arguments for Rcpp_integrate(), and memory allocation stuffs are done > inside Rcpp_integrate(). > > You may try the code below: > > > incl =' > > #include <R_ext/Applic.h> > > List Rcpp_integrate(integr_fn f, void *ex, double lower, double upper, > int subdiv = 100, double eps_abs = 1e-7, double > eps_rel = 1e-7) > { > int lenw = 4 * subdiv; > int *iwork = new int[subdiv]; > double *work = new double[lenw]; > > double value; > double abs_err; > int subdiv_used; > int neval; > int error; > > Rdqags(f, ex, &lower, &upper, &eps_abs, &eps_rel, > &value, &abs_err, &neval, &error, &subdiv, &lenw, > &subdiv_used, iwork, work); > > delete [] iwork; > delete [] work; > > return List::create(Named("value") = value, > Named("abs.err") = abs_err, > Named("subdivisions") = subdiv_used, > Named("neval") = neval); > } > > typedef struct { > double shape; > double scale; > } Params; > > void integrand(double *t, int n, void *ex) > { > Params *param = (Params *) ex; > for(int i = 0; i < n; i++) > t[i] = t[i] * R::dweibull(t[i], param->shape, param->scale, 0); > } > > ' > > cppFunction(' > > List integrate_example() > { > Params param = {2, 3}; > return Rcpp_integrate(integrand, ¶m, 0, 5); > } > > ', includes = incl) > > integrate(function(x) x * dweibull(x, 2, 3), 0, 5) > integrate_example() > > > > > Also, there may be easier ways as KK and Dirk suggested. Find one that's > most comfortable to you. > > > Best, > Yixuan > > > 2014-07-30 15:15 GMT-04:00 Silkworth,David J. <[email protected]>: > >> I have a current project desire to move something like R’s integrate >> function inside a loop in Rcpp code. The performance hit to call back to R >> seems to kill the advantage of Rcpp in the first place. Actually my >> integrand is t*pdf(t), very similar indeed to pweibull which integrates >> pdf(t). It has been hard to find previous discussion on this topic since >> the title of articles and book include the word Integration in another >> context. >> >> I realize that eventually R’s integrate function calls rdqags (over a >> definite interval), but there is a lot of memory management that is taken >> care of before the call. This is over my head. >> >> I could try to incorporate the GSL, but this too seems daunting (even >> with RcppGSL). I think there may be integration support in the boost >> headers, but my head is too small for this yet. >> >> Any ideas that could help me? >> >> Dave Silkworth >> >> >> _______________________________________________ >> Rcpp-devel mailing list >> [email protected] >> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel >> > > > > -- > Yixuan Qiu <[email protected]> > Department of Statistics, > Purdue University > > _______________________________________________ > Rcpp-devel mailing list > [email protected] > https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel >
_______________________________________________ Rcpp-devel mailing list [email protected] https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
