I am looking for a deterministic threadsafe Rcpp algorithm for 2-D numerical integration. RcppNumerical provides a partial interface to Cuba for multidimensional integration, but from my trials that appears not to be threadsafe in RcppParallel, and it probably uses a Monte Carlo method. That throws me back on repeated 1-dimensional integration. I have used this successfully with the (not threadsafe) R function Rdqags, but my (possibly naive) coding for RcppNumerical fails to compile because the nested class is abstract. Perhaps due to the operator() virtual function.
Can anyone suggest a way around this in RcppNumerical, or some alternative?
My test code emulating the 2-D example from https://github.com/yixuan/RcppNumerical is below. It gives errors like
cannot declare variable 'f2' to be of abstract type 'Normal2'
cannot declare variable 'f1' to be of abstract type 'Normal1'
Murray
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;
// P(a1 < X1 < b1, a2 < X2 < b2), (X1, X2) ~ N([0], [1 rho])
// ([0], [rho 1])
class Normal2: public Func
{
private:
const double rho;
const double x;
double const1; // 2 * (1 - rho^2)
double const2; // 1 / (2 * PI) / sqrt(1 - rho^2)
public:
Normal2(const double& rho_, const double& x_) : rho(rho_), x(x_)
{
const1 = 2.0 * (1.0 - rho * rho);
const2 = 1.0 / (2 * M_PI) / std::sqrt(1.0 - rho * rho);
}
// PDF of bivariate normal
double operator()(const double& y)
{
double z = x * x - 2 * rho * x * y + y * y;
return const2 * std::exp(-z / const1);
}
};
class Normal1: public Func
{
private:
const double rho;
double a2, b2;
public:
Normal1(const double& rho_, const double& a2_, const double& b2_) : rho(rho_), a2(a2_), b2(b2_) {}
// integral in y dimension for given x
double operator()(const double& x)
{
Normal2 f2(rho, x);
double err_est;
int err_code;
const double res = integrate(f2, a2, b2, err_est, err_code);
return res;
}
};
// [[Rcpp::export]]
Rcpp::List integrate_test3()
{
double a1 = -1.0;
double b1 = 1.0;
double a2 = -1.0;
double b2 = 1.0;
Normal1 f1(0.5, a2, b2); // rho = 0.5
double err_est;
int err_code;
const double res = integrate(f1, a1, b1, err_est, err_code);
return Rcpp::List::create(
Rcpp::Named("approximate") = res,
Rcpp::Named("error_estimate") = err_est,
Rcpp::Named("error_code") = err_code
);
}
The Numer::Func class is an abstract class because of one undefined method:
virtual double operator()(const double& x) const = 0;
Now you are providing an implementation for
double operator()(const double& x)
which leaves the above method undefined and hence the class abstract. You should change this to
double operator()(const double& x) const
for both Normal1 and Normal2 to have your code compile.
BTW, my compiler (gcc 9.2) is even quite explicit about this problem:
59094915.cpp: In member function ‘double Normal1::operator()(const double&)’:
59094915.cpp:43:17: error: cannot declare variable ‘f2’ to be of abstract type ‘Normal2’
43 | Normal2 f2(rho, x);
| ^~
59094915.cpp:9:7: note: because the following virtual functions are pure within ‘Normal2’:
9 | class Normal2: public Func
| ^~~~~~~
In file included from /usr/local/lib/R/site-library/RcppNumerical/include/integration/wrapper.h:13,
from /usr/local/lib/R/site-library/RcppNumerical/include/RcppNumerical.h:16,
from 59094915.cpp:3:
/usr/local/lib/R/site-library/RcppNumerical/include/integration/../Func.h:26:20: note: ‘virtual double Numer::Func::operator()(const double&) const’
26 | virtual double operator()(const double& x) const = 0;
| ^~~~~~~~
59094915.cpp: In function ‘Rcpp::List integrate_test3()’:
59094915.cpp:58:13: error: cannot declare variable ‘f1’ to be of abstract type ‘Normal1’
58 | Normal1 f1(0.5, a2, b2); // rho = 0.5
| ^~
59094915.cpp:32:7: note: because the following virtual functions are pure within ‘Normal1’:
32 | class Normal1: public Func
| ^~~~~~~
In file included from /usr/local/lib/R/site-library/RcppNumerical/include/integration/wrapper.h:13,
from /usr/local/lib/R/site-library/RcppNumerical/include/RcppNumerical.h:16,
from 59094915.cpp:3:
/usr/local/lib/R/site-library/RcppNumerical/include/integration/../Func.h:26:20: note: ‘virtual double Numer::Func::operator()(const double&) const’
26 | virtual double operator()(const double& x) const = 0;
| ^~~~~~~~