ex-5: implement the code for plotting σ_MC vs calls

move all the functions to lib.c and lib.h; add this two librearies to makefile;
make plot.py more easy to use for passing to one mode to the other (show or
save figure); create fit_plot and create the figure 5-fit.pdf.

ex-5/lib.c

@@ -0,0 +1,86 @@
+#include <gsl/gsl_fit.h>
+#include <stdio.h>
+#include <math.h>
+#include "lib.h"
+
+
+// Wrapper for the gsl_function structure
+//
+double function (double * x, size_t dim, void * params)
+  {
+  return exp(x[0]); 
+  }
+
+///////////////////////////////////////////////////////////////////////////////  
+
+// Results printer.
+//
+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);
+    }
+  } 
+
+///////////////////////////////////////////////////////////////////////////////  
+
+// Perform a fit in order to compare the data with the expected function:
+//
+//   y = a⋅x^b  →  ln(y) = ln(a) + b⋅ln(x)  →  Y = A + B⋅X
+//
+// with:
+//
+//   A = ln(a)  →  a = e^A
+//   B = b      →  b = B
+//
+// using a linear regression. For b, the results is hence compared with the
+// expected one:
+//
+//   b_exp = - 0.5
+//
+void fit (struct bag full_bag, double* p,
+          double* a, double* a_err,
+          double* b, double* b_err)
+  {
+  // Expected value.
+  //
+  double b_exp = -0.5;
+
+  // Parse arguments.
+  //
+  double size = full_bag.size;
+  double* x = full_bag.pokets.x;
+  double* y = full_bag.pokets.y;
+  for (size_t i = 0; i < size; i++)
+    {
+    x[i] = log(x[i]);
+    y[i] = log(y[i]);
+    }
+
+  // Do fit.
+  //
+  double A, B, A_err, B_err, AB_cov, sum2;
+  gsl_fit_linear(x, 1, y, 1, size, &A, &B, &A_err, &AB_cov, &B_err, &sum2);
+  
+  // Parse results.
+  //
+  A_err = sqrt(A_err);
+  B_err = sqrt(B_err);
+  *a = exp(A);
+  *b = B;
+  *a_err = *a * A_err;
+  *b_err = B_err;
+
+  // Check compatibility with expected values.
+  //
+  double t = fabs(*b - b_exp)/ *b_err;
+  *p = 1 - erf(t/sqrt(2));
+  }
+

ex-5/lib.h

@@ -0,0 +1,35 @@
+#pragma once
+
+
+// Subtruct for the results storage.
+//
+struct poket
+  {
+  double * x;  // Vector with the numbers of function calls.
+  double * y;  // Vector with Plain MC error estimates.
+  };
+
+// Struct for the results storage.
+//
+struct bag
+  {
+  struct poket pokets;  // Values.
+  size_t size;          // Quantity of values.
+  };
+
+///////////////////////////////////////////////////////////////////////////////
+
+// Wrapper for the gsl_function structure.
+//
+double function (double * x, size_t dim, void * params);
+
+// Results printer.
+//
+void results (size_t calls, double result, double error, double chi);
+
+// Do a linear minimization fitting the model y = a⋅x^b. The p-value of the
+// compatibility test of b with the expected value -0.5 is returned in p.
+//
+void fit (struct bag full_bag, double* p,
+          double* a, double* a_err,
+          double* b, double* b_err);

ex-5/main.c

@@ -2,43 +2,34 @@
 #include <gsl/gsl_monte_miser.h>
 #include <gsl/gsl_monte_vegas.h>
 #include <gsl/gsl_monte.h>
+#include <stdio.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);
-    }
-  } 
+#include "lib.h"
 
 
 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;
+  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.
   //
@@ -52,23 +43,18 @@ int main(int argc, char** argv)
   expo.dim = 1;
   expo.params = NULL;
 
-  // Print the results table header
+  // 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.
+  // Compute the integral with the three MC methods: plain MC, MISER and VEGAS.
   //
-  while (calls <= 5000000)
+  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);
@@ -77,6 +63,11 @@ int main(int argc, char** argv)
     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,
@@ -106,13 +97,28 @@ int main(int argc, char** argv)
     double chi = sVE->chisq;
     results(0, integral, error, chi);
 
-
     // Update function calls
+    //
     printf ("\n");
-    calls = calls*10;
+    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;

ex-5/plot.py

@@ -4,26 +4,29 @@ import matplotlib.pyplot as plt
 import numpy             as np
 import sys
 
+
 table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|')
 calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table
 exact = 1.7182818285
 
-plt.rcParams['font.size'] = 17
-
 plt.figure()
+# plt.figure(figsize=(5, 3))
+# plt.rcParams['font.size'] = 8
+
 plt.title('Plain MC', loc='right')
 plt.ylabel('$I^{oss}$')
 plt.xlabel('calls')
 plt.xscale('log')
 
-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.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.tight_layout()
+# plt.savefig('notes/images/5-test.pdf')
 plt.show()

ex-5/plots/fit_plot.py

@@ -0,0 +1,31 @@
+#!/usr/bin/env python
+
+import matplotlib.pyplot as plt
+import numpy             as np
+import sys
+
+
+table = np.loadtxt(sys.stdin, unpack=True, skiprows=2, delimiter='|')
+calls, I_MC, σ_MC, I_MI, σ_MI, I_VE, σ_VE, chi = table
+exact = 1.7182818285
+
+plt.figure()
+# plt.figure(figsize=(5, 3))
+# plt.rcParams['font.size'] = 8
+
+plt.title('Integral uncertainty', loc='right')
+plt.ylabel('$\sigma_I$')
+plt.xlabel('calls')
+plt.xscale('log')
+plt.yscale('log')
+
+plt.plot(calls, σ_MC, linestyle='', marker='o',
+         color='#92182b', label='$\sigma_I$')
+x = np.logspace(1, 8, 5)
+a = 0.48734
+b = -0.49938
+plt.plot(x, a*(x**b), color='gray')
+
+plt.tight_layout()
+# plt.savefig('notes/images/5-fit.pdf')
+plt.show()

makefile

@@ -27,7 +27,7 @@ ex-4/bin/main: ex-4/main.c ex-4/lib.c
 	$(CCOMPILE)
 
 ex-5: ex-5/bin/main
-ex-5/bin/%: ex-5/%.c
+ex-5/bin/%: ex-5/main.c ex-5/lib.c
 	$(CCOMPILE)
 
 ex-6: ex-6/bin/main ex-6/bin/test