[Rcpp-devel] Rf_integrate?

[Baptiste Auguie]

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 <yixuan.qiu at cos.name> 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, &param, 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. <SILKWODJ at airproducts.com>:
>
>>  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
>> Rcpp-devel at lists.r-forge.r-project.org
>> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>>
>
>
>
> --
> Yixuan Qiu <yixuan.qiu at cos.name>
> Department of Statistics,
> Purdue University
>
> _______________________________________________
> 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/20140731/436118d7/attachment.html>