|
Packit |
67cb25 |
#include <stdlib.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_math.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte_plain.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte_miser.h>
|
|
Packit |
67cb25 |
#include <gsl/gsl_monte_vegas.h>
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* Computation of the integral,
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
I = int (dx dy dz)/(2pi)^3 1/(1-cos(x)cos(y)cos(z))
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
over (-pi,-pi,-pi) to (+pi, +pi, +pi). The exact answer
|
|
Packit |
67cb25 |
is Gamma(1/4)^4/(4 pi^3). This example is taken from
|
|
Packit |
67cb25 |
C.Itzykson, J.M.Drouffe, "Statistical Field Theory -
|
|
Packit |
67cb25 |
Volume 1", Section 1.1, p21, which cites the original
|
|
Packit |
67cb25 |
paper M.L.Glasser, I.J.Zucker, Proc.Natl.Acad.Sci.USA 74
|
|
Packit |
67cb25 |
1800 (1977) */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
/* For simplicity we compute the integral over the region
|
|
Packit |
67cb25 |
(0,0,0) -> (pi,pi,pi) and multiply by 8 */
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double exact = 1.3932039296856768591842462603255;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double
|
|
Packit |
67cb25 |
g (double *k, size_t dim, void *params)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
(void)(dim); /* avoid unused parameter warnings */
|
|
Packit |
67cb25 |
(void)(params);
|
|
Packit |
67cb25 |
double A = 1.0 / (M_PI * M_PI * M_PI);
|
|
Packit |
67cb25 |
return A / (1.0 - cos (k[0]) * cos (k[1]) * cos (k[2]));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
void
|
|
Packit |
67cb25 |
display_results (char *title, double result, double error)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
printf ("%s ==================\n", title);
|
|
Packit |
67cb25 |
printf ("result = % .6f\n", result);
|
|
Packit |
67cb25 |
printf ("sigma = % .6f\n", error);
|
|
Packit |
67cb25 |
printf ("exact = % .6f\n", exact);
|
|
Packit |
67cb25 |
printf ("error = % .6f = %.2g sigma\n", result - exact,
|
|
Packit |
67cb25 |
fabs (result - exact) / error);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
int
|
|
Packit |
67cb25 |
main (void)
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
double res, err;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
double xl[3] = { 0, 0, 0 };
|
|
Packit |
67cb25 |
double xu[3] = { M_PI, M_PI, M_PI };
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
const gsl_rng_type *T;
|
|
Packit |
67cb25 |
gsl_rng *r;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_monte_function G = { &g, 3, 0 };
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
size_t calls = 500000;
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_env_setup ();
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
T = gsl_rng_default;
|
|
Packit |
67cb25 |
r = gsl_rng_alloc (T);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_monte_plain_state *s = gsl_monte_plain_alloc (3);
|
|
Packit |
67cb25 |
gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, s,
|
|
Packit |
67cb25 |
&res, &err;;
|
|
Packit |
67cb25 |
gsl_monte_plain_free (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
display_results ("plain", res, err);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_monte_miser_state *s = gsl_monte_miser_alloc (3);
|
|
Packit |
67cb25 |
gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, s,
|
|
Packit |
67cb25 |
&res, &err;;
|
|
Packit |
67cb25 |
gsl_monte_miser_free (s);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
display_results ("miser", res, err);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (3);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, s,
|
|
Packit |
67cb25 |
&res, &err;;
|
|
Packit |
67cb25 |
display_results ("vegas warm-up", res, err);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
printf ("converging...\n");
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
do
|
|
Packit |
67cb25 |
{
|
|
Packit |
67cb25 |
gsl_monte_vegas_integrate (&G, xl, xu, 3, calls/5, r, s,
|
|
Packit |
67cb25 |
&res, &err;;
|
|
Packit |
67cb25 |
printf ("result = % .6f sigma = % .6f "
|
|
Packit |
67cb25 |
"chisq/dof = %.1f\n", res, err, gsl_monte_vegas_chisq (s));
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
display_results ("vegas final", res, err);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_monte_vegas_free (s);
|
|
Packit |
67cb25 |
}
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
gsl_rng_free (r);
|
|
Packit |
67cb25 |
|
|
Packit |
67cb25 |
return 0;
|
|
Packit |
67cb25 |
}
|