Hi Ben,

I follow your instruction 1D numerical integration (again) and https://www.boost.org/doc/libs/1_67_0/libs/math/doc/html/quadrature.html to use boost library to write a function in .cpp with Rcpp to compute a integral of exp(x) function.

Here is the code.

#include <Rcpp.h>

#include <boost/math/quadrature/gauss_kronrod.hpp>

using namespace Rcpp;

// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]

double integral_GK(double ll, double ul,double nu) {

boost::math::quadrature::gauss_kronrod<double, 15> integrator;

auto f1 = [&](double x) {

return std::exp(x);

};

double error;

double ll_new=ll * nu;

double ul_new=ul * nu;

double Q = integrator.integrate(f1, ll_new, ul_new, 0, 0, &error)/nu;

return Q;

}

The code works well. What I’m trying is to translate this code to .hpp code as your instruction in https://cran.r-project.org/web/packages/rstan/vignettes/external.html so that I can inlcude in stan model, but I’m not successful to make the code run.

Here is the .hpp code that I tried:

template <typename T0__, typename T1__,typename T2__>

typename boost::math::tools::promote_args<T0__, T1__,T2__>::type

integral_GK(const T0__& ll, const T1__& ul,const T2__& nu, std::ostream* pstream__) {

boost::math::quadrature::gauss_kronrod<double, 15> integrator;

auto f1 = [&](double x) {

return std::exp(x);

};

double error;

double ll_new=ll * nu;

double ul_new=ul * nu;

double Q = integrator.integrate(f1, ll_new, ul_new, 0, 0, &error)/nu;

return Q;

}

and I go to the model_header.hpp and I included the header

#include <boost/math/quadrature/gauss_kronrod.hpp>

Here is a part of the error:

In file included from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/config.hpp:39:0,

from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/tools/config.hpp:13,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/var.hpp:7,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/core.hpp:12,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math/rev/mat.hpp:4,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/stan/math.hpp:4,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:4,

from file29e417e718b4.cpp:8:

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: “BOOST_NO_CXX11_RVALUE_REFERENCES” redefined

# define BOOST_NO_CXX11_RVALUE_REFERENCES

^

:0:0: note: this is the location of the previous definition

In file included from C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/quadrature/gauss_kronrod.hpp:18:0,

from C:/Users/nhatlth/Documents/R/win-library/3.3/StanHeaders/include/src/stan/model/model_header.hpp:20,

from file29e417e718b4.cpp:8:

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp: In member function ‘std::vector boost::math::legendre_stieltjes::zeros() const’:

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:205:18: error: ‘g’ does not name a type

auto g = [&](Real t) { return this->operator()(t); };

^

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:206:18: error: ‘p’ does not name a type

auto p = boost::math::tools::bisect(g, lower_bound, upper_bound, tol);

^

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:208:31: error: ‘p’ was not declared in this scope

Real x_nk_guess = p.first + (p.second - p.first)*half();

^

C:/Users/nhatlth/Documents/R/win-library/3.3/BH/include/boost/math/special_functions/legendre_stieltjes.hpp:211:18: error: ‘f’ does not name a type

auto f = [&] (Real x) { Real Pn = this->operator()(x);

^…

Can you please help me to translate this code into a code.hpp so that I can run in stan model. I know that I ask too much from you. But it seems that I’m stuck here since my C++ knowledge is very limited.

Thank you very much in advance!

Nhat