gsl-example.ipynb

Kernel: C++17

XEUS C++17 running GSL

• https://github.com/QuantStack/xeus

• https://www.gnu.org/software/gsl/doc/html/index.html

Tipp: Use the doube-arrow ("fast forward") to restart and run all cells after any changes!

#pragma cling add_include_path("/usr/include/")

#pragma cling load("gsl")
#pragma cling load("gslcblas")

#include <stdlib.h>
#include <iostream>
#include <gsl/gsl_math.h>
#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_miser.h>
#include <gsl/gsl_monte_plain.h>
#include <gsl/gsl_monte_vegas.h>

Computation of the integral,

$$I = \int_{(-pi,-pi,-pi)}^{(+pi, +pi, +pi)} \frac{\mathrm{d}x \mathrm{d}y \mathrm{d}z}{(2 \pi)^3} \frac{1}{1-cos(x) \cos(y) \cos(z)}$$

The exact answer is

$$\frac{\Gamma(1/4)^4}{4 pi^3}$$

This example is taken from C.Itzykson, J.M.Drouffe, "Statistical Field Theory - Volume 1", Section 1.1, p21, which cites the original paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74 1800 (1977)

For simplicity we compute the integral over the region (0,0,0) -> (pi,pi,pi) and multiply by 8

double exact = 1.3932039296856768591842462603255;

double g(double *k, size_t dim, void *params) {
  (void)(dim); /* avoid unused parameter warnings */
  (void)(params);
  double A = 1.0 / (M_PI * M_PI * M_PI);
  return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
}

void display_results(char *title, double result, double error) {
   std::cout << "================== " << title << " ==================\n";
   std::cout << "result = " << result << "\n";
   std::cout << "sigma  = " << error << "\n";
   std::cout << "exact  = " << exact << "\n";
   std::cout << "error  = " << result - exact << " = " << fabs(result - exact) / error << " sigma\n";
}

double res, err;

double xl[3] = {0, 0, 0};
double xu[3] = {M_PI, M_PI, M_PI};

const gsl_rng_type *T;
gsl_rng *r;

gsl_monte_function G = {&g, 3, 0};

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(3);
gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_plain_free(s);

display_results((char*)"plain", res, err);
}

{
gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
gsl_monte_miser_free(s);

display_results((char*)"miser", res, err);
}

{
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
display_results((char*)"vegas warm-up", res, err);

printf("converging...\n");

do {
  gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
  printf(
      "result = % .6f sigma = % .6f "
      "chisq/dof = %.1f\n",
      res, err, gsl_monte_vegas_chisq(s));
} while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);

display_results((char*)"vegas final", res, err);

gsl_monte_vegas_free(s);
}

gsl_rng_free(r);

================== plain ==================
result = 1.41221
sigma  = 0.0134359
exact  = 1.3932
error  = 0.0190048 = 1.41448 sigma
================== miser ==================
result = 1.39132
sigma  = 0.00346056
exact  = 1.3932
error  = -0.00188235 = 0.543944 sigma
================== vegas warm-up ==================
result = 1.39267
sigma  = 0.00341041
exact  = 1.3932
error  = -0.000531339 = 0.155799 sigma
================== vegas final ==================
result = 1.39328
sigma  = 0.00036248
exact  = 1.3932
error  = 7.74556e-05 = 0.213682 sigma