Giù Marcer
f60f1f65d3
A for-cycle has been implemented in order to produce a bunch of results with different numbers of function-calls which can now be plotted with plot.py, which has been added to the folder.
99 lines
2.5 KiB
C
99 lines
2.5 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>
|
|
|
|
|
|
// 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;
|
|
}
|