help.C

#include "random.h"
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <random>

/*
    GenRan: non-Gaussian random field generator.
    Copyright (C) 2006 Peter Grassl
    Distributed under the MIT licence.

    Mathematical helpers used throughout the program: summary moments
    (moment), uniform random numbers (ran1), log-gamma (gammaln), the
    standard-normal CDF / inverse-CDF (normal_01_cdf, normal_01_cdf_inv,
    plus shift-and-scale wrappers), polynomial evaluation (dpoly_value),
    and the grafted-Weibull-Gaussian rescaling fixed point
    (computeRescaling).

    Provenance:
      - moment, ran1, gammaln are independent or standard-library
        implementations.
      - normal_*_cdf{,_inv}, dpoly_value are adapted from John Burkardt's
        "prob" library (https://people.sc.fsu.edu/~jburkardt/), MIT.
        normal_01_cdf_inv is the Wichura AS 241 algorithm.
*/


// ---------------------------------------------------------------------------
//  moment(): mean, average deviation, standard deviation, variance, skewness
//  and kurtosis of data[1..n] (legacy 1-based indexing kept for callers).
//
//  Implemented with an online Welford / Pebay update, which is numerically
//  stable for the very long sequences this program produces (e.g. n = 2^19).
// ---------------------------------------------------------------------------
void moment(double data[], int n, double *ave, double *adev, double *sdev,
            double *var, double *skew, double *curt)
{
  if (n < 2) {
    std::fprintf(stderr, "moment: n must be at least 2\n");
    std::exit(1);
  }

  double mean = 0.0;
  double m2 = 0.0;   // sum of (x_i - mean)^2
  double m3 = 0.0;   // sum of (x_i - mean)^3
  double m4 = 0.0;   // sum of (x_i - mean)^4

  for (int i = 1; i <= n; ++i) {
    const double k = static_cast<double>(i);
    const double delta = data[i] - mean;
    const double delta_n = delta / k;
    const double delta_n2 = delta_n * delta_n;
    const double term = delta * delta_n * (k - 1.0);
    m4 += term * delta_n2 * (k * k - 3.0 * k + 3.0)
          + 6.0 * delta_n2 * m2 - 4.0 * delta_n * m3;
    m3 += term * delta_n * (k - 2.0) - 3.0 * delta_n * m2;
    m2 += term;
    mean += delta_n;
  }

  // Average absolute deviation has no online formula, so make a second pass.
  double mad = 0.0;
  for (int i = 1; i <= n; ++i) mad += std::fabs(data[i] - mean);

  *ave = mean;
  *adev = mad / n;
  const double variance = m2 / static_cast<double>(n - 1);
  *var = variance;
  *sdev = std::sqrt(variance);
  if (variance > 0.0) {
    *skew = (m3 / n) / std::pow(variance, 1.5);
    *curt = (m4 / n) / (variance * variance) - 3.0;
  } else {
    std::fprintf(stderr, "moment: cannot compute skew/kurtosis when variance is zero\n");
    std::exit(1);
  }
}


// ---------------------------------------------------------------------------
//  ran1(): drop-in replacement for the legacy uniform [0, 1) generator.
//
//  Backed by std::mt19937_64. The function preserves the call-site
//  convention used elsewhere in the program: a non-positive *idum reseeds
//  the engine with abs(*idum) (or 1 if *idum is zero) and *idum is then
//  set to 1 to mark the engine as seeded. Subsequent calls with positive
//  *idum just return the next value without touching the engine state.
// ---------------------------------------------------------------------------
double ran1(long *idum)
{
  static std::mt19937_64 engine;
  static std::uniform_real_distribution<double> dist(0.0, 1.0);
  static bool seeded = false;

  if (*idum <= 0 || !seeded) {
    long s = (*idum <= 0) ? -(*idum) : *idum;
    if (s == 0) s = 1;
    engine.seed(static_cast<std::uint64_t>(s));
    seeded = true;
    *idum = s;
  }
  return dist(engine);
}


// ---------------------------------------------------------------------------
//  gammaln(): natural log of the gamma function. Thin wrapper over the
//  C++11 standard-library implementation.
// ---------------------------------------------------------------------------
double gammaln(double xx)
{
  return std::lgamma(xx);
}


// ===========================================================================
//  The remainder of this file: normal-distribution CDF and inverse-CDF
//  functions adapted from John Burkardt's "prob" library
//  (https://people.sc.fsu.edu/~jburkardt/cpp_src/prob/prob.html) which is
//  distributed under the MIT licence. The inverse CDF is the Wichura
//  AS 241 algorithm.
// ===========================================================================

double normal_01_cdf ( double x )
{
  double a1 = 0.398942280444;
  double a2 = 0.399903438504;
  double a3 = 5.75885480458;
  double a4 = 29.8213557808;
  double a5 = 2.62433121679;
  double a6 = 48.6959930692;
  double a7 = 5.92885724438;
  double b0 = 0.398942280385;
  double b1 = 3.8052E-08;
  double b2 = 1.00000615302;
  double b3 = 3.98064794E-04;
  double b4 = 1.98615381364;
  double b5 = 0.151679116635;
  double b6 = 5.29330324926;
  double b7 = 4.8385912808;
  double b8 = 15.1508972451;
  double b9 = 0.742380924027;
  double b10 = 30.789933034;
  double b11 = 3.99019417011;
  double cdf;
  double q;
  double y;

  if (fabs(x) <= 1.28){
    y = 0.5*x*x;
    q = 0.5-fabs(x)*(a1-a2*y/(y+a3-a4/(y+a5
                                       +a6/(y+a7))));
  }
  else if (fabs(x)<=12.7){
    y=0.5*x*x;
    q=exp(-y)*b0/(fabs(x)-b1+b2/(fabs(x)+b3+b4/(fabs(x)-b5+b6/(fabs(x)+b7-b8/(fabs(x)+b9+b10/(fabs(x)+b11))))));
  }
  else{
    q = 0.0;
  }
  if(x < 0.0){
    cdf = q;
  }
  else{
    cdf = 1.-q;
  }

  return cdf;
}



double normal_cdf ( double x, double a, double b )
{
  double cdf;
  double y;

  y = ( x - a ) / b;

  cdf = normal_01_cdf ( y );

  return cdf;
}



double dpoly_value ( int n, double a[], double x )
{
  int i;
  double value;

  value = 0.0;

  for (i=n-1; 0<=i; i--){
    value = value*x + a[i];
  }

  return value;
}




double normal_cdf_inv ( double cdf, double a, double b )
{
  double x;
  double x2;

  if ( cdf < 0.0 || 1.0 < cdf )
  {
    printf( "\n");
    printf("NORMAL_CDF_INV - Fatal error!\n");
    printf("  CDF < 0 or 1 < CDF.\n");
    exit ( 1 );
  }

  x2 = normal_01_cdf_inv ( cdf );

  x = a + b * x2;

  return x;
}



double normal_01_cdf_inv ( double p )
{
  double a[8] = {
    3.3871328727963666080,     1.3314166789178437745e+2,
    1.9715909503065514427e+3,  1.3731693765509461125e+4,
    4.5921953931549871457e+4,  6.7265770927008700853e+4,
    3.3430575583588128105e+4,  2.5090809287301226727e+3 };
  double b[8] = {
    1.0,                       4.2313330701600911252e+1,
    6.8718700749205790830e+2,  5.3941960214247511077e+3,
    2.1213794301586595867e+4,  3.9307895800092710610e+4,
    2.8729085735721942674e+4,  5.2264952788528545610e+3 };
  double c[8] = {
    1.42343711074968357734,     4.63033784615654529590,
    5.76949722146069140550,     3.64784832476320460504,
    1.27045825245236838258,     2.41780725177450611770e-1,
    2.27238449892691845833e-2,  7.74545014278341407640e-4 };
  double const1 = 0.180625;
  double const2 = 1.6;
  double d[8] = {
    1.0,                        2.05319162663775882187,
    1.67638483018380384940,     6.89767334985100004550e-1,
    1.48103976427480074590e-1,  1.51986665636164571966e-2,
    5.47593808499534494600e-4,  1.05075007164441684324e-9 };
  double e[8] = {
    6.65790464350110377720,     5.46378491116411436990,
    1.78482653991729133580,     2.96560571828504891230e-1,
    2.65321895265761230930e-2,  1.24266094738807843860e-3,
    2.71155556874348757815e-5,  2.01033439929228813265e-7 };
  double f[8] = {
    1.0,                        5.99832206555887937690e-1,
    1.36929880922735805310e-1,  1.48753612908506148525e-2,
    7.86869131145613259100e-4,  1.84631831751005468180e-5,
    1.42151175831644588870e-7,  2.04426310338993978564e-15 };
  double q;
  double r;
  double split1 = 0.425;
  double split2 = 5.0;
  double value;

  if ( p <= 0.0 )
  {
    value = -HUGE_VAL;
    return value;
  }

  if ( 1.0 <= p )
  {
    value = HUGE_VAL;
    return value;
  }

  q = p - 0.5;

  if ( fabs ( q ) <= split1 )
  {
    r = const1 - q * q;
    value = q * dpoly_value ( 8, a, r ) / dpoly_value ( 8, b, r );
  }
  else
  {
    if ( q < 0.0 )
    {
      r = p;
    }
    else
    {
      r = 1.0 - p;
    }

    if ( r <= 0.0 )
    {
      value = -1.0;
      exit ( 1 );
    }

    r = sqrt ( -log ( r ) );

     if ( r <= split2 )
    {
      r = r - const2;
      value = dpoly_value ( 8, c, r ) / dpoly_value ( 8, d, r );
     }
     else
     {
       r = r - split2;
       value = dpoly_value ( 8, e, r ) / dpoly_value ( 8, f, r );
    }

    if ( q < 0.0 )
    {
      value = -value;
    }

  }

  return value;
}


// Fixed-point iteration that determines the rescaling factor, Gaussian-core
// mean, and Gaussian-core std for the grafted Weibull-Gaussian CDF
// (typeOfCDF == 3). The empirical formula linking the core std to the
// coefficient of variation comes from Bažant & Pang (2006); the
// "+ standardDeviation*squarePart" sign corrects a minus that appears in
// Pang's thesis.
#define MYPI 3.14159265

double computeRescaling(double graftedParameter[], double &mean, double &standardDeviation)
{
  const double modulus               = graftedParameter[1];
  const double graftingProbability   = graftedParameter[2];
  const double coefficientOfVariation = graftedParameter[3];
  const double tolerance = 1.e-3;
  const int    maxIterations = 200;

  double rfNew = 1., rf = 0., length = 0., helpVal = 0., gaussianProbability = 0.;
  double squarePart = 0.;
  double error = 1.;
  int iteration = 0;
  while (error > tolerance) {
    if (++iteration > maxIterations) {
      std::fprintf(stderr,
                   "computeRescaling: did not converge after %d iterations "
                   "(last error %.3e, tolerance %.3e)\n",
                   maxIterations, error, tolerance);
      std::exit(1);
    }
    rf = rfNew;
    length = pow(-log(1. - graftingProbability / rf), 1. / modulus);
    helpVal = 1000. * graftingProbability / rf / (109. * graftingProbability / rf + 0.133);
    standardDeviation = exp(-3.25 + 11.6 * coefficientOfVariation
                            - helpVal * pow(coefficientOfVariation, 2.));
    squarePart = pow(-2. * log(sqrt(2. * MYPI) * modulus * standardDeviation
                               * pow(length, modulus - 1.)
                               * exp(-pow(length, modulus))),
                     0.5);
    mean = length + standardDeviation * squarePart;
    gaussianProbability = normal_01_cdf((length - mean) / standardDeviation);
    rfNew = 1. / (1. - gaussianProbability + graftingProbability / rf);
    error = fabs(rfNew - rf);
  }
  return rfNew;
}