I am working on a large C++ framework to do some particle physics computation and it would be nice to have a method to do N-dimensional integration of a generic function with N variables. The integration happens on the N-cube [0:1]^N.
The integrand is unfortunately a member function and so I am passing an std::bind object to a templated function (lambdas work well too). When N > 1, the integration method binds again the function fixing one variable and passes the new function to the method for N-1.
I have the feeling that this can be done recursively by having N as a template parameter, but this is where I get stuck: is there a way to "templetize" the binding part such that it adds the right amount of placeholders with either the arity of the function or N? Then pass the new bound method to integral<N-1>, etc. I believe the problem lies in that an std::bind object has no signature. Maybe with lambdas one could exploit variadic pack expansions, but again the template has to resolve the signature somehow.
Is there a way to make this work? If possible in c++11 and without boost libraries.
This is a very simplified version of what I have now without recursions, but instead different definitions of the integration method for each dimension (up to 3 now)
#include <iostream>
#include <functional>
#include <cmath>
// integrate between [0:1]
template <typename Functor>
double integral1D(Functor fn, size_t steps = 100) {
double sum = 0.;
double step = 1./steps;
for (size_t s = 0; s < steps; ++s)
sum += fn(s * step);
return sum * step;
}
template <typename Functor>
double integral2D(Functor fn, size_t steps = 100) {
auto subfn = [&](double v) ->double {
auto subint = std::bind(fn, v, std::placeholders::_1);
return integral1D(subint, steps);
};
return integral1D(subfn, steps);
}
template <typename Functor>
double integral3D(Functor fn, size_t steps = 100) {
auto subfn = [&](double v) ->double {
auto subint = std::bind(fn, v, std::placeholders::_1, std::placeholders::_2);
return integral2D(subint, steps);
};
return integral1D(subfn, steps);
}
struct A
{
double gaus1(double x, double range) { // computes jacobian on [-range:range]
x = range * (2. * x - 1.);
return 2 * range * std::exp(- x*x);
}
double gaus2(double x, double y, double range) {
return gaus1(x, range) * gaus1(y, range);
}
double gaus3(double x, double y, double z, double range) {
return gaus2(x, y, range) * gaus1(z, range);
}
double gaus1_integrate(double range) {
auto func = std::bind(&A::gaus1, this, std::placeholders::_1, range);
return integral1D(func);
}
double gaus2_integrate(double range) {
auto func = std::bind(&A::gaus2, this, std::placeholders::_1,
std::placeholders::_2, range);
return integral2D(func);
}
double gaus3_integrate(double range) {
auto func = std::bind(&A::gaus3, this, std::placeholders::_1,
std::placeholders::_2,
std::placeholders::_3, range);
return integral3D(func);
}
};
int main() {
A a;
std::cout << "1D integral is " << a.gaus1_integrate(50) << "
";
std::cout << "2D integral is " << a.gaus2_integrate(50) << "
";
std::cout << "3D integral is " << a.gaus3_integrate(50) << "
";
return 0;
}
The above example works and gives expected results. I know it is not recommended to do Riemann sums (or equivalent) for more complicated functions with more dimensions, but it would be nice to see if something like described above can work.
ANSWER
Based on Guillaume's helpful suggestion, here is the c++11 working example. There were minor edits needed.
#include <iostream>
#include <functional>
#include <cmath>
template<int n, typename Functor>
struct integralNDsubint {
double v;
const Functor &fn;
template<typename... Args>
double operator()(Args... rest) const {
return fn(v, rest...);
}
};
template <int n, typename Functor, typename std::enable_if<(n > 1), bool>::type = true>
double integralND(Functor fn, size_t steps = 100) {
auto subfn = [&](double v) -> double {
auto subint = integralNDsubint<n, Functor> {v, fn};
// following block can be used in c++14 without
// need of helper struct integralNDsubint
//auto subint = [v, &fn](auto... rest) {
//return fn(v, rest...);
//};
return integralND<n - 1>(subint, steps);
};
return integralND<1>(subfn, steps);
}
template <int n, typename Functor, typename std::enable_if<(n == 1), bool>::type = true>
double integralND(Functor fn, size_t steps = 100) {
double sum = 0.;
double step = 1./steps;
for (size_t s = 0; s < steps; ++s) {
sum += fn(s * step);
}
return sum * step;
}
struct A
{
double gaus1(double x, double range) {
x = range * (2. * x - 1.);
return 2 * range * std::exp(- x*x);
}
double gaus2(double x, double y, double range) {
return gaus1(x, range) * gaus1(y, range);
}
double gaus3(double x, double y, double z, double range) {
return gaus2(x, y, range) * gaus1(z, range);
}
double gaus1_integrate(double range) {
auto func = [=](double x) { return gaus1(x, range); };
return integralND<1>(func);
}
double gaus2_integrate(double range) {
auto func = [=](double x, double y) { return gaus2(x, y, range); } ;
return integralND<2>(func);
}
double gaus3_integrate(double range) {
auto func = [=](double x, double y, double z) { return gaus3(x, y, z, range); } ;
return integralND<3>(func);
}
};
int main() {
A a;
std::cout << "1D integral is " << a.gaus1_integrate(50) << "
";
std::cout << "2D integral is " << a.gaus2_integrate(50) << "
";
std::cout << "3D integral is " << a.gaus3_integrate(50) << "
";
return 0;
}
question from:
https://stackoverflow.com/questions/65907638/write-templated-recursive-integration-method-that-accepts-generic-n-dimensional 
