/* 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` and
 * other optimisation functions because the GSL
 * routines expect 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 or large values
  if (x < -200)   x = -200;
  if (x > 100000) x = 100000;

  // 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.001);
  }

  fprintf(stderr, "QDF error: %.3g\n", upp - low);
  double root = gsl_root_fsolver_root(s);

  // free memory
  gsl_root_fsolver_free(s);

  return root;
}