diff --git a/ex-5/main.c b/ex-5/main.c index cf05f46..bfb6014 100644 --- a/ex-5/main.c +++ b/ex-5/main.c @@ -5,88 +5,91 @@ #include -// Function which must be passed to the gsl_function. -double function (double * x, size_t dim, void * params){ - // Avoid warning when no parameters are passed to the function. - (void)(params); +// 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 (char *title, double result, double error) -{ - printf ("---------------------------\n%s\n", title); - printf ("result: % .10f\n", result); - printf ("sigma: % .10f\n", error); -} - - -int main(int argc, char** argv){ - - // If no argument is given is input, an error signal is displayed - // and the program quits. - - if (argc != 2) +// +void results (size_t calls, double result, double error) + { + if (calls != 0) { - fprintf(stderr, "usage: %s D\n", argv[0]); - fprintf(stderr, "Computes the integral generating D points.\n"); - return EXIT_FAILURE; + printf ("%ld\t", calls); } + printf ("%.10f\t", result); + printf ("%.10f\t", error); + } + +int main(int argc, char** argv) + { // Some useful variables. - double integral, error; - size_t calls = 5000; - size_t dims = 1; + // + 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; - { - gsl_monte_vegas_state *s = gsl_monte_vegas_alloc (dims); - gsl_monte_vegas_integrate (&expo, lower, upper, dims, calls, - r, s, &integral, &error); - do { gsl_monte_vegas_integrate (&expo, lower, upper, dims, - calls/5, r, s, &integral, &error); - } while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5); - - // Free memory. - gsl_monte_vegas_free(s); - - // Display results - results("Vegas", integral, error); - } - - { - gsl_monte_miser_state *s = gsl_monte_miser_alloc (dims); - gsl_monte_miser_integrate (&expo, lower, upper, dims, calls, - r, s, &integral, &error); - // Free memory. - gsl_monte_miser_free (s); - - // Display results - results ("Miser", integral, error); - } - - { - gsl_monte_plain_state *s = gsl_monte_plain_alloc (dims); + // 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, s, &integral, &error); - // Display results - results ("Montecarlo", integral, error); + r, sMC, &integral, &error); + gsl_monte_plain_free(sMC); + results(calls, integral, error); - // Free memory. - gsl_monte_plain_free (s); - } + // 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); diff --git a/ex-5/plot.py b/ex-5/plot.py new file mode 100644 index 0000000..4c74c51 --- /dev/null +++ b/ex-5/plot.py @@ -0,0 +1,20 @@ +import matplotlib.pyplot as plt +import numpy as np +import sys + +calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE = np.loadtxt(sys.stdin, unpack=True) +exact = 1.7182818285 + +plt.rcParams['font.size'] = 17 +plt.figure() +plt.title('Plain MC', loc='right') +plt.ylabel('$I^{oss}$') +plt.xlabel('calls') +plt.axhline(y=exact, color='#c556ea', linestyle='-', label='Exact value') + +plt.errorbar(calls, I_MC, linestyle='', marker='o', yerr=σ_MC, color='#92182b', label='Plain MC') +plt.errorbar(calls, I_MI, linestyle='', marker='o', yerr=σ_MI, color='black', label='MISER') +plt.errorbar(calls, I_VE, linestyle='', marker='o', yerr=σ_VE, color='gray', label='VEGAS') + +plt.legend() +plt.show() diff --git a/notes/images/MC_all.pdf b/notes/images/MC_all.pdf new file mode 100644 index 0000000..6363da3 Binary files /dev/null and b/notes/images/MC_all.pdf differ diff --git a/notes/images/MC_few.pdf b/notes/images/MC_few.pdf new file mode 100644 index 0000000..72f95ad Binary files /dev/null and b/notes/images/MC_few.pdf differ diff --git a/notes/images/MC_some.pdf b/notes/images/MC_some.pdf new file mode 100644 index 0000000..f995a9a Binary files /dev/null and b/notes/images/MC_some.pdf differ diff --git a/notes/sections/3.md b/notes/sections/3.md index 3535d47..c752526 100644 --- a/notes/sections/3.md +++ b/notes/sections/3.md @@ -2,7 +2,7 @@ ## Meson decay events generation -A number of $N = 50'000$ points on the unit sphere, each representing a +A number of $N = 50000$ points on the unit sphere, each representing a particle detection event, must be generated according to then angular probability distribution function $F$: \begin{align*} diff --git a/notes/sections/4.md b/notes/sections/4.md index ed8f274..f1972dc 100644 --- a/notes/sections/4.md +++ b/notes/sections/4.md @@ -156,7 +156,7 @@ Namely: ## Monte Carlo simulation The same distribution should be found by generating and binning points in a -proper way. A number of $N = 50'000$ points were generated as a couple of values +proper way. A number of $N = 50000$ points were generated as a couple of values ($P$, $\theta$), with $P$ evenly distributed between 0 and $P_{\text{max}}$ and $\theta$ given by the same procedure described in @sec:3, namely: $$ diff --git a/notes/todo b/notes/todo index 9d88081..2742125 100644 --- a/notes/todo +++ b/notes/todo @@ -2,3 +2,4 @@ - aggiungere citazioni e referenze - rifare grafici senza bordino - leggere l'articolo di Lucy +- riscrivere il readme del 5