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


// Wrapper for the gsl_function.
//
double function (double * x, size_t dim, void * params)
  {
  return exp(x[0]); 
  }


// Function which prints out the results when called.
//
void results (size_t calls, double result, double error)
  {
  if (calls != 0)
    {
    printf ("%ld\t", calls);
    }
  printf ("%.10f\t", result);
  printf ("%.10f\t", error);
  } 


int main(int argc, char** argv)
  {
  // Some useful variables.
  //
  size_t dims = 1;          // Integral dimension
  double lower[1] = {0};    // Integration lower limit
  double upper[1] = {1};    // Integration upper limit
  size_t calls = 5000;        // Initial number of function calls
  double integral, error;

  // Initialize an RNG.
  //
  gsl_rng_env_setup();
  gsl_rng *r = gsl_rng_alloc(gsl_rng_default);

  // Define GSL function.
  //
  gsl_monte_function expo = {&function, dims, NULL};
  expo.f = &function;
  expo.dim = 1;
  expo.params = NULL;

  // Compute the integral for the following number of function calls:
  //          50
  //         500
  //       5'000
  //      50'000
  //     500'000
  //   5'000'000
  //  50'000'000
  // with the three MC methods: plain MC, MISER and VEGAS.
  //
  while (calls < 50000001)
    {
    // Plain MC
    gsl_monte_plain_state *sMC = gsl_monte_plain_alloc (dims);
    gsl_monte_plain_integrate (&expo, lower, upper, dims, calls,
                               r, sMC, &integral, &error);
    gsl_monte_plain_free(sMC);
    results(calls, integral, error);

    // MISER
    gsl_monte_miser_state *sMI = gsl_monte_miser_alloc (dims);
    gsl_monte_miser_integrate (&expo, lower, upper, dims, calls,
                               r, sMI, &integral, &error);
    gsl_monte_miser_free(sMI);
    results(0, integral, error);

    // VEGAS
    gsl_monte_vegas_state *sVE = gsl_monte_vegas_alloc (dims);
    gsl_monte_vegas_integrate (&expo, lower, upper, dims, calls,
                               r, sVE, &integral, &error);
    do
      {
      gsl_monte_vegas_integrate (&expo, lower, upper, dims,
                                 calls/5, r, sVE, &integral, &error);
      } while (fabs (gsl_monte_vegas_chisq (sVE) - 1.0) > 0.5);
    gsl_monte_vegas_free(sVE);
    results(0, integral, error);

    // Update function calls
    printf ("\n");
    calls = calls*10;
    }

  // Free memory.
  gsl_rng_free(r);

  return EXIT_SUCCESS;
}