Using GSL Monte Carlo Integration in C++ for a multidimensional function

Question

I'm trying to use GSL Monte Carlo Integration in a C++ code that I'm generating. The idea would be to have a brownian motion function (brownian), which is used in another function (g) for performing 4-dim numerical integration.

double brownian(const double &x,const double &x0,const double & sigma,const double & t) {
   double a1 = (1/sqrt(2.0 * M_PI * sigma * t));
   double b1 = exp(-((x - x0) * (x - x0))/(2.0 * sigma * t));
   double res = a1 * b1;
  return res;
}

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

 const double& xA0 = *(double *)p[0];
 const double& xB0 = *(double *)p[1];
 const double& yA0 = *(double *)p[2];
 const double& yB0 = *(double *)p[3];
 const double& sigma = *(double *)p[4];
 const double& t1 = *(double *)p[5];


  double temp_pbxA = brownian(k[0], xA0, sigma,t1);
  double temp_pbxB = brownian(k[1], xB0, sigma, t1);
  double temp_pbyA = brownian(k[2], yA0, sigma,t1);
  double temp_pbyB = brownian(k[3], yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);
}

double integrate_test(const double xl, const double xu, const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){

  double res, err;

 const gsl_rng_type *T;
  gsl_rng *r;

  gsl_monte_function G;
  G.f = &g;
  G.dim = 4;
  G.params = { xA0, xB0, yA0, yB0, sigma, t1};

  size_t calls = 500000;

  gsl_rng_env_setup ();

  T = gsl_rng_default;
  r = gsl_rng_alloc (T);


    gsl_monte_plain_state *s = gsl_monte_plain_alloc (4);
    gsl_monte_plain_integrate (&G, xl, xu, 4, calls, r, s,
                               &res, &err);
    gsl_monte_plain_free (s);

return res;
    }

However, when I try to compile the code I get the following error:

error: assigning to 'double (*)(double *, size_t, void *)' (aka 'double (*)(double *, unsigned long, void *)') from incompatible type 'double (*)(double *, size_t, void **)' (aka 'double (*)(double *, unsigned long, void **)'): type mismatch at 3rd parameter ('void *' vs 'void **')
  G.f = &g;

I do not understand why it takes the void *p[] as it was a void**. Any suggestion?

Answer 1

This is a diff between your code and the one that compiles. This does not mean that my code is correct: you only get what you asked for - since you did not provide a complete minimum working example, anyone trying to help you is forced to guess what your problem really is.

13c13
< double g(double *k, size_t dim, void *p){
---
> double g(double *k, size_t dim, void *p[]){
15,20c15,20
<  const double& xA0 = ((double *)p)[0];
<  const double& xB0 = ((double *)p)[1];
<  const double& yA0 = ((double *)p)[2];
<  const double& yB0 = ((double *)p)[3];
<  const double& sigma = ((double *)p)[4];
<  const double& t1 = ((double *)p)[5];
---
>  const double& xA0 = *(double *)p[0];
>  const double& xB0 = *(double *)p[1];
>  const double& yA0 = *(double *)p[2];
>  const double& yB0 = *(double *)p[3];
>  const double& sigma = *(double *)p[4];
>  const double& t1 = *(double *)p[5];
31c31
< double integrate_test(const double xl[], const double xu[], const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){
---
> double integrate_test(const double xl, const double xu, const double& xA0, const double& xB0, const double& yA0, const double& yB0, const double& t1, const double& sigma){
41,42c41
<   double params[] = { xA0, xB0, yA0, yB0, sigma, t1};
<   G.params = params;
---
>   G.params = { xA0, xB0, yA0, yB0, sigma, t1};

As you can see, there are at least 2 problems that your compiler complains about.

1. gsl_monte_function expect the component f to have the signature double (* f) (double * x, size_t dim, void * params), see https://www.gnu.org/software/gsl/doc/html/montecarlo.html . I must stress that normally parameters in GSL are passed via a void pointer to struct that must be defined somewhere beforehand, but here you're lucky, all parameters are just doubles, so to save extra code you can use an old, plain C array

double params[] = { xA0, xB0, yA0, yB0, sigma, t1};

and pass its address as the G.params. Then you'll retrieve them inside g using C-syntax like this

const double& xA0 = ((double *)p)[0];

(convert the void pointer to double* and then treat it as a name of a C-style array).

2. gsl_monte_plain_integrate expects its 2nd and 3rd parameters to be C-style arrays. See https://www.gnu.org/software/gsl/doc/html/montecarlo.html . That's why you need to change the types of xl and xu in integrate_test.

GSL is hard to use. Don't hesitate to ask if in doubt.