#include <iostream>
#include <iomanip>
#include <string>

#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>

using namespace std;

double func (double *x, size_t dim, void *params);
void display_results (string title, double result, double error);

int main ()
{
  const size_t dim=3;
  double result,error;
  string method;
  gsl_rng *r;
  gsl_monte_function G = { &func, dim, 0 };
  double a[dim] = { 0, 0, 0 };
  double b[dim] = { M_PI, M_PI, M_PI };
  size_t calls = 500000;


  // random number generator  
  gsl_rng_env_setup ();                
  r = gsl_rng_alloc (gsl_rng_default);

  {
  method="plain";
  gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim);
  gsl_monte_plain_integrate (&G, a, b, dim, calls, r, s, 
                                    &result, &error);
  gsl_monte_plain_free (s);
  display_results (method , result, error);
  }


  {
  method="MISER";
  gsl_monte_miser_state *s = gsl_monte_miser_alloc (dim);
  gsl_monte_miser_integrate (&G, a, b, dim, calls, r, s,
			     &result, &error);
  gsl_monte_miser_free (s);
  display_results (method , result, error);
  }


  {
  method="VEGAS warmup";
  gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (dim);
  // warmup
  gsl_monte_vegas_integrate (&G, a, b, dim, 10000, r, s,
			     &result, &error);
  display_results (method, result, error);
  method="VEGAS";


  do
    {
      gsl_monte_vegas_integrate (&G, a, b, dim, calls/5, r, s,
				 &result, &error);
      display_results (method, result, error);
      method = "VEGAS continue";
    }
  while (fabs(gsl_monte_vegas_chisq(s)-1.0) > 0.5);
  
  gsl_monte_vegas_free (s);
  }
  return 0;
}
double func (double *x, size_t dim, void *params)
{  
  double A = 1.0 / (M_PI * M_PI * M_PI);
  return A/(1.0 - cos (x[0]) * cos (x[1]) * cos (x[2]));
}
void display_results (string title, double result, double error)
{
  const double exact = 1.3932039296856768591842462603255;

  cout<<setiosflags(ios::fixed);
  cout<<setfill('=')<<setw(50)<<"="<<endl;
  cout<<"Method : "<<title<<endl;
  cout<<setfill('=')<<setw(50)<<"="<<endl;
  cout.precision(8);
  cout<<"result = "<<result<<" +/- "<<error<<endl;
  cout<<" exact = "<<exact<<endl;
  cout<<"|diff| = "<<fabs(result-exact)<<endl;
}