Implementing multidimensional integral for a custom function in C++

Question

I am not an expert with C++, but I am trying to implement a 4-dimensional integral using GSL numerical integration approach.

The code below shows the whole algorithm. Although it seems correct what I wrote, I have some difficulties to get results. It actually takes very long time.

Do you have any suggestion on what could be the issue? Also, any suggestion on how I could improve my code?

Here below, I describe in detail the algorithm.

Headings necessary for the implementation of the algorithm. GSL headers are used for implementing the 4-dimensional integral:

#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include<algorithm>
#include <limits>
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_vegas.h>

using namespace std;
using namespace Numer;

PDFfunction represents the probability density function that is used within the integrand function of the 4-dimensional integral:

double PDFfunction(double L, int t, double tau, double x0, float x, int n_lim) {

    const double c = (M_PI/2) * (M_PI/2) * ((2 * t) / tau);

    double res = 0;

    for(int n = 1; n <= n_lim; ++n){

        res += exp(-1 * (n * n) * c) * cos((n * M_PI * x) / L) * cos((n * M_PI * x0) / L);

    }

    return (1 / L) + ((2 / L) * res);
}

Composite_at_t is a function that makes use of the PDFfunction to compute pbA and pbB:

double Composite_at_t(double t, double B, double x0, double xt_pos, double y0, double yt_pos, int Ltot, double tau, int n_lim) {

double pbA = PDFfunction(Ltot, t, tau, x0, xt_pos, n_lim);
double pbB = PDFfunction(Ltot, t, tau, y0, yt_pos, n_lim);
double b1 = 2 * (2 * t) * exp(-2 * t * B);
return pbA * pbB * b1;
}

Composite_at_tCLASS is a Func class which computes a first integral over variable t.

class Composite_at_tCLASS: public Func{
private:
    double B;
    double x0;
    double xt_pos;
    double y0;
    double yt_pos;
    int Ltot;
    double tau;
    int n_lim;
public:
    Composite_at_tCLASS(double B_, double x0_, double xt_pos_, double y0_, double yt_pos_, int Ltot_, double tau_, int n_lim_) : B(B_), x0(x0_), xt_pos(xt_pos_), y0(y0_), yt_pos(yt_pos_), Ltot(Ltot_), tau(tau_), n_lim(n_lim_) {}
    double operator()(const double& t) const{
        return Composite_at_t(t, B, x0, xt_pos, y0, yt_pos, Ltot, tau, n_lim);
    }
};

integrate_CompositeCLASS is the actual function that uses the class Composite_at_tCLASS and perform the integral over t, between 0 and Infinity.

double integrate_CompositeCLASS(double B, double x0, double xt_pos, double y0, double yt_pos, int Ltot, double tau, int n_lim){

    Composite_at_tCLASS f(B, x0, xt_pos, y0, yt_pos, Ltot, tau, n_lim);
    double err_est;
    int err_code;
    double inf = std::numeric_limits<double>::infinity();
    double res = integrate(f, 0, inf, err_est, err_code);
    return res;
}

Since I want to compute an integral of the integrate_CompositeCLASS function over another variable (this time variable B), I create another Func class that would allow me to do so:

class CompositeCLASS_overB: public Func{
private:
    double x0;
    double xt_pos;
    double y0;
    double yt_pos;
    int Ltot;
    double tau;
    int n_lim;
public:
    CompositeCLASS_overB(double x0_, double xt_pos_, double y0_, double yt_pos_, int Ltot_, double tau_, int n_lim_) : x0(x0_), xt_pos(xt_pos_), y0(y0_), yt_pos(yt_pos_), Ltot(Ltot_), tau(tau_), n_lim(n_lim_) {}
    double operator()(const double& B) const{
        return integrate_CompositeCLASS(B, x0, xt_pos, y0, yt_pos, Ltot, tau, n_lim);
    }
};

integrate_CompositeCLASS_overB is the actual integral function, that uses the CompositeCLASS_overB class to perform the integral over B (in my test case, the range of the integral is between 0.01 and 0.02).

  double integrate_CompositeCLASS_overB(double lowB, double highB, double x0, double xt_pos, double y0, double yt_pos, int Ltot, double tau, int n_lim){

    CompositeCLASS_overB f(x0, xt_pos, y0, yt_pos, Ltot, tau, n_lim);
    double err_est;
    int err_code;
    double res = integrate(f, lowB, highB, err_est, err_code);
    return res;
}

Now I define the list of fixed parameters (my_f_params) that will be used in the integrand function of the 4-dimensional integral. g function defines the 4-dimensional integral, where the variables to integrate are given by k[0], k[1], k[2] and k[3].

struct my_f_params { double lowB; double highB; double tau; int n_lim; int Ltot; double ref_xA; double ref_xB; double ref_yA; double ref_yB;};

double g(double *k, size_t dim, void *p){

struct my_f_params * fp = (struct my_f_params *)p;

double* xAt{ new double };
double* xBt{ new double };
double* yAt{ new double };
double* yBt{ new double };

*xAt = k[0] + fp->ref_xA;
*xBt = k[1] + fp->ref_xB;
*yAt = k[2] + fp->ref_yA;
*yBt = k[3] + fp->ref_yB;

double res_time = integrate_CompositeCLASS_overB(fp->lowB, fp->highB, *xAt, *xBt, *yAt, *yBt, fp->Ltot, fp->tau, fp->n_lim);

delete xAt;
delete xBt;
delete yAt;
delete yBt;

return res_time;
}

Ultimately, integrate_integral is the actual 4-dimensional integral. In this case, and as mentioned above, the 4-dimensional integral is implemented using GSL (VEGAS algorithm): https://www.gnu.org/software/gsl/doc/html/montecarlo.html. In this integral, the range of integration is given by the arrays xl={0, 0, 0, 0} and xu={4, 4, 10, 10}, which are specified inside the function.

double integrate_integral(const double& lowB, const double& highB, const double& tau, const int& n_lim, const int& Ltot,
const double& ref_xA, const double& ref_xB, const double& ref_yA, const double& ref_yB){


double res, err;

double xl[4] = {0, 0, 0, 0};
double xu[4] = {4, 4, 10, 10};
const gsl_rng_type *T;
gsl_rng *r;

gsl_monte_function G;
struct my_f_params params = { lowB, highB, tau, n_lim, Ltot, ref_xA, ref_xB, ref_yA, ref_yB};

G.f = &g;
G.dim = 4;
G.params = &params;

size_t calls = 10000;

gsl_rng_env_setup ();

T = gsl_rng_default;
r = gsl_rng_alloc (T);


  gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (4);
  gsl_monte_vegas_integrate (&G, xl, xu, 4, 10000, r, s,
                             &res, &err);

  do
        {
          gsl_monte_vegas_integrate (&G, xl, xu, 4, calls/5, r, s,
                                     &res, &err);
        }
      while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);

      gsl_monte_vegas_free (s);
      gsl_rng_free (r);


return res;
}

Answer 1

Unnecessary memory allocations

I don't understand how you came to use new to allocate four doubles inside g(). There is absolutely no reason for that as far as I can tell, instead just can just declare them on the stack:

double xAt = k[0] + fp->ref_xA;
double XBt = ...
...
return integrate_CompositeCLASS_overB(fp->lowB, fp->highB, xAt, xBt, yAt, yBt, fp->Ltot, fp->tau, fp->n_lim);

Optimize the PDF

From what I can see, the most expensive part is likely the evaluation of the probability density function, as it has lots of divisions and transcedental functions to evaluate. Assuming you don't have to worry about infinities and denormals, consider compiling your code with -ffast-math, this allows the compiler to use a lot more optimization on maths.

My math is rusty, but you could consider expanding cos() into terms of cexp() and see if you can simplify the expression inside the loop that way. The next thing to consider is to get rid of divisions as much as possible. By calculating 1 / L up front and then multiplying by that instead of dividing by L already saves 3 div instructions, even when using -ffast-math. Even better would be if you would feed PDFfunction() with the reciprocals of L and tau, so it doesn't have to do the divisions itself.

I see the same PDF is used for both x and y. If you were integrating at regular intervals with the same domain for x and y, I would suggest using memoization, but if you are doing Monte Carlo integration that probably won't help. Another option, if the PDF is smooth enough, is to precalculate it for lots of points in the domain used in the integration, and then to use a simple interpolation to return the estimated PDF at arbitrary positions.

Answer 2

Following the suggestion of @G.Sliepen, I have edited the PDFfunction to reduce divisions and replace cos(x) with an expansion of two complex exponentials, such that:

cos(x) = 1/2 * ( cexp(I * x) + cexp(-I * x) )

Since in PDFfunction two cos(x) are multiplied, we can replace the product of two exponentials as one exponential with argument equal to the sum of the arguments of the two exponentials:

cos(x) * cos(y) = 1/2 * ( cexp(I * x) + cexp(-I * x) ) * 1/2 * ( cexp(I * y) + cexp(-I * y) );

which then becomes:

cos(x) * cos(y) = 1/4 * cexp( (I * x) + (I * y) ) + cexp( (I * x) + (-I * y) ) + cexp( (-I * x) + (I * y) ) + cexp( (-I * x) + (-I * y) )

Using these expressions within PDFfunction and using the C++ implementation of complex exponential, we obtain:

const complex<double> I{ 0.0, 1.0 };

double PDFfunction(double invL, int t, double invtau, double x0, double x, double n_lim) {

    const double c =  M_PI * (M_PI/4) * ((2 * t) * invtau);
    complex<double> res = 0;

    for(double n = 1; n <= n_lim; ++n){

      res += 0.25 * (exp((-1 * n * n * c)  + (I * n * M_PI * x * invL) + (I * n * M_PI * x0 * invL)) +  exp((-1 * n * n * c)  + (I * n * M_PI * x * invL) + (-I * n * M_PI * x0 * invL)) +
      exp((-1 * n * n * c)  + (-I * n * M_PI * x * invL) + (I * n * M_PI * x0 * invL)) + exp((-1 * n * n * c)  + (-I * n * M_PI * x * invL) + (-I * n * M_PI * x0 * invL)));
    }



    return invL + ((2 * invL) * std::real(res));
}

However, the replacement of cos(x) with complex exponential did not really speed up my code. You may want to give it a try and see if I am wrong.