I'm trying to implement a 4-dimensional numerical integration in a C++ code to be, then, sourced in R using the Rcpp package.
I have followed the the following tutorial Rcpp Numerical Integration, and I managed to get some results. However, when I compare these results with those I obtain when using cubature R package (cubintegrate function, method='cuhre'), I do get very different results. Do you know why? And also, would you suggest an alternative way to perform fast 4-dimensional integration in C++/Rcpp?
Below there is my code:
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <RcppNumerical.h>
#include <Rcpp.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>
using namespace std;
using namespace Rcpp;
using namespace Numer;
// [[Rcpp::depends(BH)]]
// [[Rcpp::export]]
double gaussian_free_diffusion(double x,double x0, double sigma, double t) {
double pi = 2 * acos(0.0);
double a1 = (1/sqrt(2.0 * pi * sigma * t));
double b1 = exp(-pow((x - x0), 2.0)/(2.0 * sigma * t));
double res = a1 * b1;
return res;
}
// [[Rcpp::export]]
double integral_CppFunctionCUBATURE(NumericVector Zt_pos, const double& xA0, const double& xB0, const double& yA0,
const double& yB0, const double& t1, const double& sigma){
double xAt_pos = Zt_pos[0];
double xBt_pos = Zt_pos[1];
double yAt_pos = Zt_pos[2];
double yBt_pos = Zt_pos[3];
double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);
return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);
};
class integral_CppFunction: public MFunc{
private:
const double xA0;
const double xB0;
const double yA0;
const double yB0;
const double t1;
const double sigma;
public:
integral_CppFunction(const double& xA0_, const double& xB0_,
const double& yA0_, const double& yB0_,const double& t1_, const double& sigma_) : xA0(xA0_), xB0(xB0_),
yA0(yA0_), yB0(yB0_), t1(t1_), sigma(sigma_) {}
double operator()(Constvec& Zt_pos)
{
double xAt_pos = Zt_pos[0];
double xBt_pos = Zt_pos[1];
double yAt_pos = Zt_pos[2];
double yBt_pos = Zt_pos[3];
double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);
return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);
}
};
// [[Rcpp::export]]
Rcpp::List integrate_CppFunction(Eigen::VectorXd lower, Eigen::VectorXd upper, const double& xA0, const double& xB0, const double& yA0,
const double& yB0, const double& t1, const double& sigma)
{
integral_CppFunction f(xA0, xB0, yA0, yB0, t1, sigma);
double err_est;
int err_code;
const double res = integrate(f, lower, upper, err_est, err_code);
return Rcpp::List::create(
Rcpp::Named("approximate") = res,
Rcpp::Named("error_estimate") = err_est,
Rcpp::Named("error_code") = err_code
);
};
integrate_CppFunction is the C++ function for integration using RcppEigen as suggested in Rcpp Numerical Integration, and integral_CppFunctionCUBATURE is the function I use with cubintegrate. In R, I would then do:
sourceCpp('./myCppcode.cpp')
xA0<-3
xB0<-2
yA0<-5
yB0<-10
t<-500
sigma<-1
integrate_CppFunction(lower=rep(-1000,4), upper=rep(1000,4), xA0=xA0, xB0=xB0, yA0=yA0,
yB0=yB0, t1=t, sigma=sigma)$approximate
Result: -106814.2
cubature::cubintegrate(integral_CppFunctionCUBATURE,lower=rep(-1000,4), upper=rep(1000,4),
xA0=xA0, xB0=xB0, yA0=yA0,yB0=yB0, t1=t, sigma=sigma,method='cuhre')$integral
Result: 0.9999978
I have explored the details of the integrate function implemented im RcppNumerical, and it seems to me that one possible explanation for the difference in the results between RcppNnumerical and cubature is the maxeval argument, corresponding to the maximum number of function evaluations needed. In the integrate function of RcppNnumerical, maxeval is set to 1000, whereas in cubintegrate maxeval is set to 10^6.
I have tried to add the maxeval argument in the integrate_CppFunction function:
// [[Rcpp::export]]
Rcpp::List integrate_CppFunction(Eigen::VectorXd lower, Eigen::VectorXd upper, const double& xA0, const double& xB0, const double& yA0,
const double& yB0, const double& t1, const double& sigma, const int maxeval)
{
integral_CppFunction f(xA0, xB0, yA0, yB0, t1, sigma);
double err_est;
int err_code;
const double res = integrate(f, lower, upper, err_est, err_code, maxeval);
return Rcpp::List::create(
Rcpp::Named("approximate") = res,
Rcpp::Named("error_estimate") = err_est,
Rcpp::Named("error_code") = err_code
);
};
but with no success. Any thought on this?
