<div dir="ltr">Nice, thanks for sharing!<div><br></div><div>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.</div>
<div><br></div><div><a href="https://github.com/baptiste/cubature/blob/master/minimal.c">https://github.com/baptiste/cubature/blob/master/minimal.c</a><br></div><div><br></div><div>(probably very poorly implemented, I was very surprised when I got it to work!)</div>
<div><br></div><div>I use it in my planar package to do a 2D integration of a 3D complex vector field (thus returning 6 real components),</div><div><br></div><div><a href="https://github.com/baptiste/planar/blob/master/src/gaussian_beam.cpp#L478">https://github.com/baptiste/planar/blob/master/src/gaussian_beam.cpp#L478</a><br>
</div><div><br></div><div>Cheers,</div><div><br></div><div>b.</div><div><br></div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On 30 July 2014 22:25, Yixuan Qiu <span dir="ltr"><<a href="mailto:yixuan.qiu@cos.name" target="_blank">yixuan.qiu@cos.name</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><div><div><div>Hi Dave,<br>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().<br>

<br></div>You may try the code below:<br><br><span style="font-family:courier new,monospace"><br>incl ='<br><br>#include <R_ext/Applic.h><br><br>List Rcpp_integrate(integr_fn f, void *ex, double lower, double upper,<br>

                    int subdiv = 100, double eps_abs = 1e-7, double eps_rel = 1e-7)<br>{<br>    int lenw = 4 * subdiv;<br>    int *iwork = new int[subdiv];<br>    double *work = new double[lenw];<br>    <br>    double value;<br>

    double abs_err;<br>    int subdiv_used;<br>    int neval;<br>    int error;<br>    <br>    Rdqags(f, ex, &lower, &upper, &eps_abs, &eps_rel,<br>           &value, &abs_err, &neval, &error, &subdiv, &lenw,<br>

           &subdiv_used, iwork, work);<br>    <br>    delete [] iwork;<br>    delete [] work;<br>    <br>    return List::create(Named("value") = value,<br>                        Named("abs.err") = abs_err,<br>

                        Named("subdivisions") = subdiv_used,<br>                        Named("neval") = neval);<br>}<br><br>typedef struct {<br>    double shape;<br>    double scale;<br>} Params;<br>
<br>
void integrand(double *t, int n, void *ex)<br>{<br>    Params *param = (Params *) ex;<br>    for(int i = 0; i < n; i++)<br>        t[i] = t[i] * R::dweibull(t[i], param->shape, param->scale, 0);<br>}<br><br>'<br>

<br>cppFunction('<br><br>List integrate_example()<br>{<br>    Params param = {2, 3};<br>    return Rcpp_integrate(integrand, &param, 0, 5);<br>}<br><br>', includes = incl)<br><br>integrate(function(x) x * dweibull(x, 2, 3), 0, 5)<br>

integrate_example()</span><br><br><br><br><br></div>Also, there may be easier ways as KK and Dirk suggested. Find one that's most comfortable to you.<br><br><br></div>Best,<br></div>Yixuan<br></div><div class="gmail_extra">

<br><br><div class="gmail_quote">2014-07-30 15:15 GMT-04:00 Silkworth,David J. <span dir="ltr"><<a href="mailto:SILKWODJ@airproducts.com" target="_blank">SILKWODJ@airproducts.com</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div><div class="h5">







<div>
<font face="Calibri"><span style="font-size:11pt">
<div>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.</div>
<div> </div>
<div>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.</div>
<div> </div>
<div>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.</div>
<div> </div>
<div>Any ideas that could help me?</div>
<div> </div>
<div>Dave Silkworth</div>
<div> </div>
</span></font>
</div>

<br></div></div><div class="">_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></div></blockquote></div><span class="HOEnZb"><font color="#888888"><br>
<br clear="all"><br>
-- <br>Yixuan Qiu <<a href="mailto:yixuan.qiu@cos.name" target="_blank">yixuan.qiu@cos.name</a>><br>Department of Statistics,<br>Purdue University<br>
</font></span></div>
<br>_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br></blockquote></div><br></div>