#include <stdlib.h>

#include <gsl/gsl_math.h>

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

#include <gsl/gsl_sf_gamma.h>

/*
 * Computation of the integral,
 *
 *     I = int (dx dy dz)/(2pi)^3  1/(1-cos(x)cos(y)cos(z))
 *
 * over (-pi,-pi,-pi) to (+pi, +pi, +pi).  The exact answer
 * is 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)
{
  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, double exact)
{
  printf ("==== %s ==============\n\n", title);
  printf ("  result    = % .6f\n", result);
  printf ("  err.est.  = % .6f\n", error);
  printf ("  exact     = % .6f\n", exact);
  printf ("  error     = % .6f = %.2g err.est.\n", result - exact,
          fabs (result - exact) / error);
  printf ("\n");
}

int main (void)
{
  double res, err;

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

  double exact = 0.25*pow(gsl_sf_gamma (.25), 4)/pow(M_PI, 3);    /* Gamma(1/4)^4/(4 pi^3) */

  gsl_rng *r;
  unsigned long seed = 1UL;

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

  size_t calls = 500000;

  /* allocate random number generator */
  r = gsl_rng_alloc(gsl_rng_taus2);

  /* seed the random number generator */
  gsl_rng_set(r, seed);

  /* monte plain */
  gsl_monte_plain_state *sp = gsl_monte_plain_alloc (3);
  gsl_monte_plain_integrate (&G, xl, xu, 3, calls, r, sp, &res, &err);
  gsl_monte_plain_free (sp);

  display_results ("plain", res, err, exact);

  /* monte miser */
  gsl_monte_miser_state *sm = gsl_monte_miser_alloc (3);
  gsl_monte_miser_integrate (&G, xl, xu, 3, calls, r, sm, &res, &err);
  gsl_monte_miser_free (sm);

  display_results ("miser", res, err, exact);

  /* monte vegas */
  gsl_monte_vegas_state *sv = gsl_monte_vegas_alloc (3);

  gsl_monte_vegas_integrate (&G, xl, xu, 3, 10000, r, sv, &res, &err);
  display_results ("vegas warm-up", res, err, exact);

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

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

  display_results ("vegas final", res, err, exact);

  gsl_monte_vegas_free (sv);

  gsl_rng_free (r);

  return 0;
}