From 0ee7f9adb31ad6128c6400462e4f16eb70d540bb Mon Sep 17 00:00:00 2001
From: rnhmjoj <rnhmjoj@inventati.org>
Date: Tue, 24 Mar 2020 01:17:39 +0000
Subject: [PATCH] ex-1: split into several files

---
 ex-1/landau.c |  95 +++++++++++++++++
 ex-1/landau.h |  24 +++++
 ex-1/main.c   | 285 +++++++-------------------------------------------
 ex-1/tests.c  | 158 ++++++++++++++++++++++++++++
 ex-1/tests.h  |  28 +++++
 makefile      |   4 +-
 6 files changed, 346 insertions(+), 248 deletions(-)
 create mode 100644 ex-1/landau.c
 create mode 100644 ex-1/landau.h
 create mode 100644 ex-1/tests.c
 create mode 100644 ex-1/tests.h

diff --git a/ex-1/landau.c b/ex-1/landau.c
new file mode 100644
index 0000000..9f7b3d2
--- /dev/null
+++ b/ex-1/landau.c
@@ -0,0 +1,95 @@
+/* This file contains functions to
+ * compute the Landau distribution
+ * PDF, CDF and QDF functions.
+ */
+
+#include <gsl/gsl_roots.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_integration.h>
+
+
+/* This is a wrapper needed by `landau_cdf` because
+ * the numerical integration expects a function
+ * with parameters.
+ */
+double landau_pdf(double x, void* params) {
+  return gsl_ran_landau_pdf(x);
+}
+
+
+/* The cumulative function of the Landau distribution
+ * calculated by numerical integration.
+ */
+double landau_cdf(double x, void* params) {
+  // create the integrand
+  gsl_function pdf;
+  pdf.function = &landau_pdf;
+  pdf.params   = NULL;
+
+  // set up the integration
+  double res, err;
+  size_t iter = 500;
+  gsl_integration_workspace* w =
+    gsl_integration_workspace_alloc(iter);
+
+  // clip values too small
+  if (x < -200) x = -200;
+
+  // We integrate the pdf in [x, +∞) instead
+  // of (-∞, x] to avoid a singularity and
+  // then return 1 - P.
+  gsl_integration_qagiu(
+    &pdf,        // integrand
+    x,           // lower bound
+    1e-10, 1e-6, // abs and rel error
+    iter,        // max iteration
+    w,           // "workspace"
+    &res, &err); // result, abs error
+
+  // free the memory
+  gsl_integration_workspace_free(w);
+
+  return 1 - res;
+}
+
+
+/* Another wrapper: this time it is needed by
+ * `landau_qdf` to solve the equation:
+ * `landau_cdf(x) = p0` for x
+ */
+double landau_cdf_root(double x, void* params) {
+  double p0 = *(double*)params;
+  return p0 - landau_cdf(x, NULL);
+}
+
+
+/* The quantile function (inverse CDF) of the Landau
+ * distribution calculated by numerical root method.
+ */
+double landau_qdf(double p0) {
+  // create function
+  gsl_function cdf;
+  cdf.function = &landau_cdf_root;
+  cdf.params   = &p0;
+
+  // use the Brent method
+  gsl_root_fsolver* s =
+    gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
+
+  // search interval
+  double low = -1000,   // lower value
+         upp =  100000; // upper value
+  gsl_root_fsolver_set(s, &cdf, low, upp);
+
+  // iterative search
+  size_t iter = 1000; // max iteration
+  int stat = GSL_CONTINUE;
+  for (size_t i=0; stat==GSL_CONTINUE && i<iter; i++) {
+    stat = gsl_root_fsolver_iterate(s);
+    low  = gsl_root_fsolver_x_lower(s);
+    upp  = gsl_root_fsolver_x_upper(s);
+    stat = gsl_root_test_interval(low, upp, 0, 0.01);
+  }
+
+  return gsl_root_fsolver_root(s);
+}
diff --git a/ex-1/landau.h b/ex-1/landau.h
new file mode 100644
index 0000000..19b5a5e
--- /dev/null
+++ b/ex-1/landau.h
@@ -0,0 +1,24 @@
+/* This file contains functions to
+ * compute the Landau distribution
+ * PDF, CDF and QDF functions.
+ */
+
+#pragma once
+
+/* This is a wrapper needed by `landau_cdf` because
+ * the numerical integration expects a function
+ * with parameters.
+ */
+double landau_pdf(double x, void* params);
+
+
+/* The cumulative function of the Landau distribution
+ * calculated by numerical integration.
+ */
+double landau_cdf(double x, void* params);
+
+
+/* The quantile function (inverse CDF) of the Landau
+ * distribution calculated by numerical root method.
+ */
+double landau_qdf(double p0);
diff --git a/ex-1/main.c b/ex-1/main.c
index 4341cf2..c45396e 100644
--- a/ex-1/main.c
+++ b/ex-1/main.c
@@ -1,11 +1,11 @@
 #include <stdio.h>
 #include <stdlib.h>
-#include <gsl/gsl_min.h>
-#include <gsl/gsl_sum.h>
-#include <gsl/gsl_roots.h>
+#include <math.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_histogram.h>
-#include <gsl/gsl_integration.h>
+
+#include "landau.h"
+#include "tests.h"
 
 
 /* Function that compare doubles for sorting:
@@ -20,233 +20,6 @@ int cmp_double (const void *xp, const void *yp) {
 }
 
 
-/* This is a wrapper needed by `landau_cdf` because
- * the numerical integration expects a function
- * with parameters.
- */
-double landau_pdf(double x, void* params) {
-  return gsl_ran_landau_pdf(x);
-}
-
-
-/* The cumulative function of the Landau distribution
- * calculated by numerical integration.
- */
-double landau_cdf(double x, void* params) {
-  // create the integrand
-  gsl_function pdf;
-  pdf.function = &landau_pdf;
-  pdf.params   = NULL;
-
-  // set up the integration
-  double res, err;
-  size_t iter = 500;
-  gsl_integration_workspace* w =
-    gsl_integration_workspace_alloc(iter);
-
-  // clip values too small
-  if (x < -200) x = -200;
-
-  // We integrate the pdf in [x, +∞) instead
-  // of (-∞, x] to avoid a singularity and
-  // then return 1 - P.
-  gsl_integration_qagiu(
-    &pdf,        // integrand
-    x,           // lower bound
-    1e-10, 1e-6, // abs and rel error
-    iter,        // max iteration
-    w,           // "workspace"
-    &res, &err); // result, abs error
-
-  // free the memory
-  gsl_integration_workspace_free(w);
-
-  return 1 - res;
-}
-
-
-/* Another wrapper: this time it is needed by
- * `landau_qdf` to solve the equation:
- * `landau_cdf(x) = p0` for x
- */
-double landau_cdf_root(double x, void* params) {
-  double p0 = *(double*)params;
-  return p0 - landau_cdf(x, NULL);
-}
-
-
-/* The quantile function (inverse CDF) of the Landau
- * distribution calculated by numerical root method.
- */
-double landau_qdf(double p0) {
-  // create function
-  gsl_function cdf;
-  cdf.function = &landau_cdf_root;
-  cdf.params   = &p0;
-
-  // use the Brent method
-  gsl_root_fsolver* s =
-    gsl_root_fsolver_alloc(gsl_root_fsolver_brent);
-
-  // search interval
-  double low = -1000,   // lower value
-         upp =  100000; // upper value
-  gsl_root_fsolver_set(s, &cdf, low, upp);
-
-  // iterative search
-  size_t iter = 1000; // max iteration
-  int stat = GSL_CONTINUE;
-  for (size_t i=0; stat==GSL_CONTINUE && i<iter; i++) {
-    stat = gsl_root_fsolver_iterate(s);
-    low  = gsl_root_fsolver_x_lower(s);
-    upp  = gsl_root_fsolver_x_upper(s);
-    stat = gsl_root_test_interval(low, upp, 0, 0.01);
-  }
-
-  return gsl_root_fsolver_root(s);
-}
-
-
-/* Kolmogorov distribution CDF
- * for sample size n and statistic D
- */
-double kolmogorov_cdf(double D, int n) {
-  double x = sqrt(n) * D;
-
-  // trick to reduce estimate error
-  x += 1/(6 * sqrt(n)) + (x - 1)/(4 * n);
-
-  // calculate the first n_terms of the series
-  // Σ_k=1 exp(-(2k - 1)²π²/8x²)
-  size_t n_terms = 30;
-  double *terms = calloc(n_terms, sizeof(double));
-  for (size_t k=1; k<n_terms; k++) {
-    terms[k] = exp(-pow((2*(double)k - 1)*M_PI/x, 2) / 8);
-  }
-
-  // do a transform to accelerate the convergence
-  double sum, abserr;
-  gsl_sum_levin_utrunc_workspace* s = gsl_sum_levin_utrunc_alloc(n_terms);
-  gsl_sum_levin_utrunc_accel(terms, n_terms, s, &sum, &abserr);
-
-  fprintf(stderr, "\n# Kolmogorov CDF\n");
-  fprintf(stderr, "accel sum: %f\n", sum);
-  fprintf(stderr, "plain sum: %f\n", s->sum_plain);
-  fprintf(stderr, "err: %f\n",   abserr);
-
-  gsl_sum_levin_utrunc_free(s);
-  free(terms);
-
-  return sqrt(2*M_PI)/x * sum;
-}
-
-
-/* This is a wrapper needed by `numeric_mode` because
- * the minimization expects a function to be minimized and not
- * maximized.
- */
-double neg_landau_pdf(double x, void* params) {
-  return (-1) * gsl_ran_landau_pdf(x);
-}
-
-
-/* Compute numerically the mode of the Landau distribution.
- */
-double numeric_mode (double min, double max) {
-  // create funtion
-  gsl_function pdf;
-  pdf.function = &neg_landau_pdf;
-  pdf.params   = NULL;
-
-  // initialize minimization
-  double guess = 0;
-  int    iter = 0;
-  int    max_iter = 100;
-  double prec = 1e-7;
-  int    status;
-  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
-  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc (T);
-  gsl_min_fminimizer_set(s, &pdf, guess, min, max);
-
-  // minimization
-  do {
-      iter++;
-      status = gsl_min_fminimizer_iterate (s);
-      guess = gsl_min_fminimizer_x_minimum (s);
-      min = gsl_min_fminimizer_x_lower (s);
-      max = gsl_min_fminimizer_x_upper (s);
-      status = gsl_min_test_interval (min, max, prec, prec);
-
-  } while (status == GSL_CONTINUE && iter < max_iter);
-
-  // Free memory
-  gsl_min_fminimizer_free (s);
-  return guess;
-}
-
-
-/* This is the function to be minimized in `numeric_FWHM`.
- */
-double abs_landau_pdf(double x, void* params_) {
-  double* params = ((double *) params_);
-  return fabs(gsl_ran_landau_pdf(x) - params[0]);
-}
-
-
-/* Compute numerically the FWHM of the Landau distribution.
- */
-double numeric_fwhm (double min_lim, double mode, double max_lim) {
-  // create funtion
-  gsl_function pdf;
-  pdf.function = &abs_landau_pdf;
-  double params [1];
-  params[0]= gsl_ran_landau_pdf(mode)/2;
-  pdf.params = params;
-
-  // initialize minimization for x₋
-  double guess = mode - 1;
-  double min, max;
-  int    iter = 0;
-  int    max_iter = 100;
-  double prec = 1e-7;
-  int    status;
-  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
-  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc (T);
-  gsl_min_fminimizer_set(s, &pdf, guess, min_lim, mode);
-
-  // minimization
-  do {
-      iter++;
-      status = gsl_min_fminimizer_iterate (s);
-      guess = gsl_min_fminimizer_x_minimum (s);
-      min = gsl_min_fminimizer_x_lower (s);
-      max = gsl_min_fminimizer_x_upper (s);
-      status = gsl_min_test_interval (min, max, prec, prec);
-  } while (status == GSL_CONTINUE && iter < max_iter);
-  double x_low = guess;
-
-  // initialize minimization for x₊
-  guess = mode + 1;
-  gsl_min_fminimizer_set(s, &pdf, guess, mode, max_lim);
-
-  // minimization
-  do {
-      iter++;
-      status = gsl_min_fminimizer_iterate (s);
-      guess = gsl_min_fminimizer_x_minimum (s);
-      min = gsl_min_fminimizer_x_lower (s);
-      max = gsl_min_fminimizer_x_upper (s);
-      status = gsl_min_test_interval (min, max, prec, prec);
-
-  } while (status == GSL_CONTINUE && iter < max_iter);
-  double x_upp = guess;
-
-  // Free memory
-  gsl_min_fminimizer_free (s);
-  return x_upp - x_low;
-}
-
-
 /* Here we generate random numbers in a uniform 
  * range and by using the quantile we map them
  * to a Landau distribution. Then we generate an
@@ -266,20 +39,31 @@ int main(int argc, char** argv) {
   gsl_histogram* hist = gsl_histogram_alloc(bins);
   gsl_histogram_set_ranges_uniform(hist, min, max);
 
-  /* sample points from the Landau distribution
-   * and fill the histogram
+
+  /* Sample generation
+   *
+   * Sample points from the Landau
+   * distribution and fill the histogram.
    */
+  fprintf(stderr, "# Sampling\n");
+  fprintf(stderr, "generating %ld points... ", samples);
   double x;
   for(size_t i=0; i<samples; i++) {
     x = gsl_ran_landau(r);
     sample[i] = x;
     gsl_histogram_increment(hist, x);
   }
+  fprintf(stderr, "done\n");
 
   // sort the sample
   qsort(sample, samples, sizeof(double), &cmp_double);
 
-  // calculate Kolmogorov-Smirnov D
+
+  /* Kolmogorov-Smirnov test
+   *
+   * Compute the D statistic and its
+   * associated probability.
+   */
   double D = 0;
   double d;
   for(size_t i=0; i<samples; i++) {
@@ -287,16 +71,19 @@ int main(int argc, char** argv) {
     if (d > D)
       D = d;
   }
-
-  // Kolmogorov-Smirnov test
+  fprintf(stderr, "\n\n# Kolmogorov-Smirnov test\n");
   double beta = kolmogorov_cdf(D, samples);
-  fprintf(stderr, "\n# K-S test results\n");
+
+  // print the results
+  fprintf(stderr, "\n## Results\n");
   fprintf(stderr, "D: %f\n", D);
   fprintf(stderr, "α: %g\n", 1 - beta);
 
-  // Mode comparison
-  //
-  // find the bin with the maximum number of events
+
+  /* Mode comparison
+   *
+   * Find the bin with the maximum number of events
+   */
   double mode_o, maxbin = 0;
   double f_mode_o = 0;
   double m1, m2   = 0;
@@ -308,19 +95,22 @@ int main(int argc, char** argv) {
       f_mode_o = hist->bin[i];
     }
   }
-  fprintf(stderr, "\nstep:%.2f\n ", (max - min)/bins);
+  fprintf(stderr, "\n\n# Mode comparison\n");
+  fprintf(stderr, "\nstep: %.2f\n ", (max - min)/bins);
   f_mode_o = f_mode_o/samples;
   mode_o = min + (maxbin + 0.5)*(max - min)/bins;
 
   // print the results
   double mode_e = numeric_mode(min, max);
-  fprintf(stderr, "\n# Numeric mode\n");
+  fprintf(stderr, "\n## Results\n");
   fprintf(stderr, "expected mode: %.7f\n", mode_e);
   fprintf(stderr, "observed mode: %.3f\n", mode_o);
 
-  // FWHM comparison
-  //
-  // find the bins x₋ and x₊
+
+  /* FWHM comparison
+   *
+   * Find the bins x₋ and x₊.
+   */
   double half = f_mode_o*samples/2;
   m2 = samples;
   double x_low = 0;
@@ -348,8 +138,9 @@ int main(int argc, char** argv) {
   double fwhm_o = x_upp - x_low;
 
   // print the results
-  fprintf(stderr, "\n# Numeric FWHM\n");
-  double fwhm_e = numeric_fwhm(min, mode_e, max);
+  fprintf(stderr, "\n\n# FWHM comparison\n");
+  double fwhm_e = numeric_fwhm(min, max, mode_e);
+  fprintf(stderr, "\n# Results\n");
   fprintf(stderr, "expected FWHM: %.7f\n", fwhm_e);
   fprintf(stderr, "observed FWHM: %.3f\n", fwhm_o);
 
diff --git a/ex-1/tests.c b/ex-1/tests.c
new file mode 100644
index 0000000..09bb892
--- /dev/null
+++ b/ex-1/tests.c
@@ -0,0 +1,158 @@
+/* This file contains functions to perform
+ * statistical tests on the points sampled
+ * from a Landau distribution.
+ */
+
+#include <stdio.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_min.h>
+#include <gsl/gsl_sum.h>
+
+#include "landau.h"
+
+
+/* Kolmogorov distribution CDF
+ * for sample size n and statistic D
+ */
+double kolmogorov_cdf(double D, int n) {
+  double x = sqrt(n) * D;
+
+  // trick to reduce estimate error
+  x += 1/(6 * sqrt(n)) + (x - 1)/(4 * n);
+
+  // calculate the first n_terms of the series
+  // Σ_k=1 exp(-(2k - 1)²π²/8x²)
+  size_t n_terms = 30;
+  double *terms = calloc(n_terms, sizeof(double));
+  for (size_t k=0; k<n_terms; k++) {
+    terms[k] = exp(-pow((2*(double)(k + 1) - 1)*M_PI/x, 2) / 8);
+  }
+
+  // do a transform to accelerate the convergence
+  double sum, abserr;
+  gsl_sum_levin_utrunc_workspace* s = gsl_sum_levin_utrunc_alloc(n_terms);
+  gsl_sum_levin_utrunc_accel(terms, n_terms, s, &sum, &abserr);
+
+  fprintf(stderr, "\n## Kolmogorov CDF\n");
+  fprintf(stderr, "accel sum: %f\n", sum);
+  fprintf(stderr, "plain sum: %f\n", s->sum_plain);
+  fprintf(stderr, "err: %f\n",   abserr);
+
+  gsl_sum_levin_utrunc_free(s);
+  free(terms);
+
+  return sqrt(2*M_PI)/x * sum;
+}
+
+
+/* This is a wrapper needed by `numeric_mode` because
+ * the minimization expects a function to be minimized and not
+ * maximized.
+ */
+double neg_landau_pdf(double x, void* params) {
+  return (-1) * gsl_ran_landau_pdf(x);
+}
+
+
+/* Numerically computes the mode of a Landau
+ * distribution by maximising the derivative.
+ * The min,max parameters are the initial search
+ * interval for the optimisation.
+ */
+double numeric_mode(double min, double max) {
+  // create funtion
+  gsl_function pdf;
+  pdf.function = &neg_landau_pdf;
+  pdf.params   = NULL;
+
+  // initialize minimization
+  double guess = 0;
+  int    iter = 0;
+  int    max_iter = 100;
+  double prec = 1e-7;
+  int    status;
+  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
+  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc(T);
+  gsl_min_fminimizer_set(s, &pdf, guess, min, max);
+
+  // minimization
+  do {
+      iter++;
+      status = gsl_min_fminimizer_iterate(s);
+      guess = gsl_min_fminimizer_x_minimum(s);
+      min = gsl_min_fminimizer_x_lower(s);
+      max = gsl_min_fminimizer_x_upper(s);
+      status = gsl_min_test_interval(min, max, prec, prec);
+
+  } while (status == GSL_CONTINUE && iter < max_iter);
+
+  // Free memory
+  gsl_min_fminimizer_free(s);
+  return guess;
+}
+
+
+/* This is the function to be minimized in `numeric_FWHM`.
+ */
+double abs_landau_pdf(double x, void* params_) {
+  double* params = ((double *) params_);
+  return fabs(gsl_ran_landau_pdf(x) - params[0]);
+}
+
+
+/* Numerically computes the FWHM of Landau
+ * distribution by maximising the derivative.
+ * The `min,max` parameters are the initial search
+ * interval for the optimisation. `mode` can be
+ * computer with `numeric_mode(min, max)`.
+ */
+double numeric_fwhm(double min, double max, double mode) {
+  // create funtion
+  gsl_function pdf;
+  pdf.function = &abs_landau_pdf;
+  double params [1];
+  params[0]= gsl_ran_landau_pdf(mode)/2;
+  pdf.params = params;
+
+  // initialize minimization for x₋
+  double guess = mode - 1;
+  double fmin, fmax;
+  int    iter = 0;
+  int    max_iter = 100;
+  double prec = 1e-7;
+  int    status;
+  const gsl_min_fminimizer_type * T = gsl_min_fminimizer_goldensection;
+  gsl_min_fminimizer * s = gsl_min_fminimizer_alloc(T);
+  gsl_min_fminimizer_set(s, &pdf, guess, min, mode);
+
+  // minimization
+  do {
+      iter++;
+      status = gsl_min_fminimizer_iterate(s);
+      guess = gsl_min_fminimizer_x_minimum(s);
+      fmin = gsl_min_fminimizer_x_lower(s);
+      fmax = gsl_min_fminimizer_x_upper(s);
+      status = gsl_min_test_interval(fmin, fmax, prec, prec);
+  } while (status == GSL_CONTINUE && iter < max_iter);
+  double x_low = guess;
+
+  // initialize minimization for x₊
+  guess = mode + 1;
+  gsl_min_fminimizer_set(s, &pdf, guess, mode, max);
+
+  // minimization
+  do {
+      iter++;
+      status = gsl_min_fminimizer_iterate(s);
+      guess = gsl_min_fminimizer_x_minimum(s);
+      fmin = gsl_min_fminimizer_x_lower(s);
+      fmax = gsl_min_fminimizer_x_upper(s);
+      status = gsl_min_test_interval(fmin, fmax, prec, prec);
+
+  } while (status == GSL_CONTINUE && iter < max_iter);
+  double x_upp = guess;
+
+  // Free memory
+  gsl_min_fminimizer_free(s);
+  return x_upp - x_low;
+}
diff --git a/ex-1/tests.h b/ex-1/tests.h
new file mode 100644
index 0000000..17aed2f
--- /dev/null
+++ b/ex-1/tests.h
@@ -0,0 +1,28 @@
+/* This file contains functions to perform
+ * statistical tests on the points sampled
+ * from a Landau distribution.
+ */
+
+#pragma once
+
+/* Kolmogorov distribution CDF
+ * for sample size n and statistic D
+ */
+double kolmogorov_cdf(double D, int n);
+
+
+/* Numerically computes the mode of a Landau
+ * distribution by maximising the derivative.
+ * The min,max parameters are the initial search
+ * interval for the optimisation.
+ */
+double numeric_mode(double min, double max);
+
+
+/* Numerically computes the FWHM of Landau
+ * distribution by maximising the derivative.
+ * The `min,max` parameters are the initial search
+ * interval for the optimisation. `mode` can be
+ * computer with `numeric_mode(min, max)`.
+ */
+double numeric_fwhm(double min, double max, double mode);
diff --git a/makefile b/makefile
index 8dadabb..9ce2042 100644
--- a/makefile
+++ b/makefile
@@ -6,7 +6,9 @@ CCOMPILE = \
 	mkdir -p $(@D); \
 	$(CC) $(CFLAGS) $^ -o $@
 
-ex-1/bin/%: ex-1/%.c
+ex-1/bin/main: ex-1/main.c ex-1/landau.c ex-1/tests.c
+	$(CCOMPILE)
+ex-1/bin/pdf: ex-1/pdf.c
 	$(CCOMPILE)
 
 ex-2/bin/%: ex-2/%.c