[Rcpp-devel] Rf_integrate?

[Yixuan Qiu]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140730/7d59dae7/attachment-0001.html>