120 lines
3.1 KiB
C
120 lines
3.1 KiB
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>
|
|
|
|
|
|
// exp() wrapper for the gsl_function structure
|
|
//
|
|
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 ("%6.0e | ", (double)calls);
|
|
}
|
|
printf ("%5f | ", result);
|
|
printf ("%5f | ", error);
|
|
if (chi != 0)
|
|
{
|
|
printf ("%5f", 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;
|
|
|
|
// Print the results table header
|
|
printf(" Calls | Plain | Error | MISER |"
|
|
" Error | VEGAS | Error | χ²\n");
|
|
printf(" ------|----------|----------|----------|"
|
|
"----------|----------|----------|---------\n");
|
|
|
|
// 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 <= 5000000)
|
|
{
|
|
// 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, ¶ms);
|
|
// params.verbose = 1;
|
|
// params.ostream = stderr;
|
|
// gsl_monte_vegas_params_set(sVE, ¶ms);
|
|
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;
|
|
}
|