CoCalc Public Filestest.c
Author: Harald Schilly
Views : 454
1#include <stdlib.h>
2#include <gsl/gsl_math.h>
3#include <gsl/gsl_monte.h>
4#include <gsl/gsl_monte_miser.h>
5#include <gsl/gsl_monte_plain.h>
6#include <gsl/gsl_monte_vegas.h>
7
8/* Computation of the integral,
9
10      I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))
11
12   over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
13   is Gamma(1/4)^4/(4 pi^3).  This example is taken from
14   C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
15   Volume 1", Section 1.1, p21, which cites the original
16   paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
17   1800 (1977) */
18
19/* For simplicity we compute the integral over the region
20   (0,0,0) -> (pi,pi,pi) and multiply by 8 */
21
22double exact = 1.3932039296856768591842462603255;
23
24double g(double *k, size_t dim, void *params) {
25  (void)(dim); /* avoid unused parameter warnings */
26  (void)(params);
27  double A = 1.0 / (M_PI * M_PI * M_PI);
28  return A / (1.0 - cos(k[0]) * cos(k[1]) * cos(k[2]));
29}
30
31void display_results(char *title, double result, double error) {
32  printf("%s ==================\n", title);
33  printf("result = % .6f\n", result);
34  printf("sigma  = % .6f\n", error);
35  printf("exact  = % .6f\n", exact);
36  printf("error  = % .6f = %.2g sigma\n", result - exact,
37         fabs(result - exact) / error);
38}
39
40int main(void) {
41  double res, err;
42
43  double xl[3] = {0, 0, 0};
44  double xu[3] = {M_PI, M_PI, M_PI};
45
46  const gsl_rng_type *T;
47  gsl_rng *r;
48
49  gsl_monte_function G = {&g, 3, 0};
50
51  size_t calls = 500000;
52
53  gsl_rng_env_setup();
54
55  T = gsl_rng_default;
56  r = gsl_rng_alloc(T);
57
58  {
59    gsl_monte_plain_state *s = gsl_monte_plain_alloc(3);
60    gsl_monte_plain_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
61    gsl_monte_plain_free(s);
62
63    display_results("plain", res, err);
64  }
65
66  {
67    gsl_monte_miser_state *s = gsl_monte_miser_alloc(3);
68    gsl_monte_miser_integrate(&G, xl, xu, 3, calls, r, s, &res, &err);
69    gsl_monte_miser_free(s);
70
71    display_results("miser", res, err);
72  }
73
74  {
75    gsl_monte_vegas_state *s = gsl_monte_vegas_alloc(3);
76    gsl_monte_vegas_integrate(&G, xl, xu, 3, 10000, r, s, &res, &err);
77    display_results("vegas warm-up", res, err);
78
79    printf("converging...\n");
80
81    do {
82      gsl_monte_vegas_integrate(&G, xl, xu, 3, calls / 5, r, s, &res, &err);
83      printf(
84          "result = % .6f sigma = % .6f "
85          "chisq/dof = %.1f\n",
86          res, err, gsl_monte_vegas_chisq(s));
87    } while (fabs(gsl_monte_vegas_chisq(s) - 1.0) > 0.5);
88
89    display_results("vegas final", res, err);
90
91    gsl_monte_vegas_free(s);
92  }
93
94  gsl_rng_free(r);
95
96  return 0;
97}
98