analistica/ex-5/main.c

#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, double chi)
  {
  if (calls != 0)
    {
    printf ("%.0e\t", (double)calls);
    }
  printf ("%.10f\t", result);
  printf ("%.10f\t", error);
  if (chi != 0)
    {
    printf ("%.3f\t", chi);
    }
  } 


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 = 50;        // 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, 0);

    // 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, 0);

    // VEGAS
    gsl_monte_vegas_state *sVE = gsl_monte_vegas_alloc (dims);
    gsl_monte_vegas_integrate (&expo, lower, upper, dims, calls,
                               r, sVE, &integral, &error);
      
    // In order to see the parameters of the VEGAS integration, decomment this
    // paragraph:
    //
    // gsl_monte_vegas_params params; 
    // gsl_monte_vegas_params_get(sVE, &params);
    // params.verbose = 1;
    // params.ostream = stderr;
    // gsl_monte_vegas_params_set(sVE, &params);
    do
      {
      gsl_monte_vegas_integrate (&expo, lower, upper, dims,
                                 calls, r, sVE, &integral, &error);
      } while (fabs (gsl_monte_vegas_chisq (sVE) - 1.0) > 0.5);
    gsl_monte_vegas_free(sVE);
    double chi = sVE->chisq;
    results(0, integral, error, chi);


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

  // Free memory.
  gsl_rng_free(r);

  return EXIT_SUCCESS;
}