ex-5: modified

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.

ex-5/main.c

@@ -5,87 +5,90 @@
 #include <math.h>
 
 
-// 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)
+//
+void results (size_t calls, double result, double error)
   {
-  printf ("---------------------------\n%s\n", title);
-  printf ("result:  % .10f\n", result);
-  printf ("sigma:   % .10f\n", error);
+  if (calls != 0)
+    {
+    printf ("%ld\t", calls);
+    }
+  printf ("%.10f\t", result);
+  printf ("%.10f\t", 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)
+int main(int argc, char** argv)
   {
-    fprintf(stderr, "usage: %s D\n", argv[0]);
-    fprintf(stderr, "Computes the integral generating D points.\n");
-    return EXIT_FAILURE;
-    }
-
   // 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;
 
+  // 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)
     {
-    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);
+    // 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.

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()

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*}

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:
 $$

notes/todo

@@ -2,3 +2,4 @@
 - aggiungere citazioni e referenze
 - rifare grafici senza bordino
 - leggere l'articolo di Lucy
+- riscrivere il readme del 5