htswanalysis/MACS/lib/gsl/gsl-1.11/doc/examples/monte.c

   1 #include <stdlib.h>
   2 #include <gsl/gsl_math.h>
   3 #include <gsl/gsl_monte.h>
   4 #include <gsl/gsl_monte_plain.h>
   5 #include <gsl/gsl_monte_miser.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
  22 double exact = 1.3932039296856768591842462603255;
  23
  24 double
  25 g (double *k, size_t dim, void *params)
  26 {
  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
  31 void
  32 display_results (char *title, double result, double error)
  33 {
  34   printf ("%s ==================\n", title);
  35   printf ("result = % .6f\n", result);
  36   printf ("sigma  = % .6f\n", error);
  37   printf ("exact  = % .6f\n", exact);
  38   printf ("error  = % .6f = %.1g sigma\n", result - exact,
  39           fabs (result - exact) / error);
  40 }
  41
  42 int
  43 main (void)
  44 {
  45   double res, err;
  46
  47   double xl[3] = { 0, 0, 0 };
  48   double xu[3] = { M_PI, M_PI, M_PI };
  49
  50   const gsl_rng_type *T;
  51   gsl_rng *r;
  52
  53   gsl_monte_function G = { &g, 3, 0 };
  54
  55   size_t calls = 500000;
  56
  57   gsl_rng_env_setup ();
  58
  59   T = gsl_rng_default;
  60   r = gsl_rng_alloc (T);
  61
  62   {
  63     gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
  64     gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s,
  65                                &res, &err);
  66     gsl_monte_plain_free (s);
  67
  68     display_results ("plain", res, err);
  69   }
  70
  71   {
  72     gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
  73     gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
  74                                &res, &err);
  75     gsl_monte_miser_free (s);
  76
  77     display_results ("miser", res, err);
  78   }
  79
  80   {
  81     gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
  82
  83     gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
  84                                &res, &err);
  85     display_results ("vegas warm-up", res, err);
  86
  87     printf ("converging...\n");
  88
  89     do
  90       {
  91         gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
  92                                    &res, &err);
  93         printf ("result = % .6f sigma = % .6f "
  94                 "chisq/dof = %.1f\n", res, err, s->chisq);
  95       }
  96     while (fabs (s->chisq - 1.0) > 0.5);
  97
  98     display_results ("vegas final", res, err);
  99
 100     gsl_monte_vegas_free (s);
 101   }
 102
 103   gsl_rng_free (r);
 104
 105   return 0;
 106 }