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 <stdio.h>
#include <math.h>
#include "lib.h"


int main(int argc, char** argv)
  {
  // Some useful variables.
  //
  size_t dims = 1;               // Integral dimension
  double lower[1] = {0};         // Integration range lower limit
  double upper[1] = {1};         // Integration range upper limit
  double integral, error;        // Result and error

  // An integral is estimated with a number c_i of function calls.
  // Different number of function calls are tested.
  //
  size_t c_0 = 50;                                   // Initial number of function calls
  size_t c_f = 50000000;                             // Final number of function calls
  size_t factor = 10;                                // c_(i+1) = c_i * factor
  size_t size = round(1 + log(c_f/c_0)/log(factor)); // Number of integrations

  // A fit will be performed. This struct is needed to accomplish it.
  //
  struct bag full_bag;
  full_bag.pokets.x = calloc(size, sizeof(double));
  full_bag.pokets.y = calloc(size, sizeof(double));
  full_bag.size = size;

  // 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;

  // Print the results table header.
  //
  printf(" Calls |  Plain   |  Error   |  MISER   |"
         "  Error   |  VEGAS   |  Error   | χ²\n");
  printf(" ------|----------|----------|----------|"
         "----------|----------|----------|---------\n");

  // Compute the integral with the three MC methods: plain MC, MISER and VEGAS.
  //
  size_t calls = c_0;
  size_t i = 0;
  while (calls <= c_f)
    {
    // 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);

    // Update the struct for the fit.
    //
    full_bag.pokets.x[i] = calls;
    full_bag.pokets.y[i++] = 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, 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*factor;
    }

  // Do a fit of the Plain MC errors.
  //
  double p, a, a_err, b, b_err;
  fit(full_bag, &p, &a, &a_err, &b, &b_err);

  // Print the fit results.
  //
  fprintf (stderr, "\n## Fit results:\n\n");
  fprintf (stderr, "a = %.5f\t",  a);
  fprintf (stderr, "δa = %.5f\n", a_err);
  fprintf (stderr, "b = %.5f\t",  b);
  fprintf (stderr, "δb = %.5f\n", b_err);
  fprintf (stderr, "p-value = %.5f\n", p);

  // Free memory.
  //
  gsl_rng_free(r);

  return EXIT_SUCCESS;
}