// g++ -std=c++11 -Wall gsl_monte.carlo.cpp `gsl-config --libs`
#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>

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

int main ()
{
  const size_t dim=3;
  double result,error;
  std::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";
    auto 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";
    auto 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";
    auto 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);
  }
}
double func (double *x, size_t dim, void *params)
{  
  static const 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 (std::string title, double result, double error)
{
  const double exact = 1.3932039296856768591842462603255;
  std::cout<<std::setiosflags(std::ios::fixed);
  std::cout<<std::setfill('=')<<std::setw(50)<<"="<<"\n";
  std::cout<<"Method : "<<title<<"\n";
  std::cout<<std::setfill('=')<<std::setw(50)<<"=\n";
  std::cout.precision(8);
  std::cout<<"result = "<<result<<" +/- "<<error<<"\n";
  std::cout<<" exact = "<<exact<<"\n";
  std::cout<<"|diff| = "<<fabs(result-exact)<<"\n";
}