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