ldes_ssj_variates.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_LDES_LDES_SSJ_VARIATES_H
#define LINE_SOLVERS_LDES_LDES_SSJ_VARIATES_H
 
/**
 * @file
 * @ingroup line_solvers
 * The variate layer of the Java LDES engine, reproduced: SSJ's `randvar`
 * generators as inverse-CDF functions of one uniform.
 *
 * THE MEASUREMENT THAT SHAPED THIS FILE. Every `*Gen` instance `Solver_ssj`
 * builds -- Exponential, Erlang, Uniform, Weibull, Pareto, Lognormal, Gamma,
 * Poisson, Binomial, Bernoulli -- consumes EXACTLY ONE uniform per draw. That
 * was measured, not assumed: a probe generated five variates from a seeded
 * MRG32k3a and then searched the raw stream for the next value, and every
 * family had advanced the stream by exactly five. SSJ's instance generators
 * invert; none of them rejects. The consequence is the whole reason a
 * bit-compatible C++ engine is feasible: the stream can never desynchronize, so
 * a family whose quantile this file computes to 1e-13 rather than to the last
 * bit still leaves every LATER draw identical, and only perturbs its own sample
 * below the tolerance of any event comparison.
 *
 * WHAT IS EXACT AND WHAT IS TO 1e-13. Exponential, Uniform, Pareto, Weibull and
 * Bernoulli have closed-form quantiles and agree with SSJ to the last bit.
 * Normal (hence Lognormal) uses Wichura's AS 241, Gamma and Erlang inverting
 * the regularized incomplete gamma by Newton with a Wilson-Hilferty start, and
 * Poisson and Binomial by discrete inversion -- the two discrete families return
 * INTEGERS and so are exact wherever the cdf comparison lands on the same side,
 * which it does everywhere the test looks.
 *
 * WHY NOT CALL SSJ'S OWN FORMULAE BLINDLY. `ExponentialDist.inverseF` is
 * `-log1p(-u)/lambda` and not `-log(1-u)/lambda`; the two differ in the last
 * bits for small u, which is exactly where an interarrival time matters most.
 * Where a spelling like that decides agreement, this file uses SSJ's.
 *
 * AND THE SPELLING IS NOT ENOUGH: IT MUST BE THE SAME log1p. `std::log1p` is
 * allowed a 1-ulp error, and glibc changed which representative it returns
 * between 2.35 and 2.39, so the SAME `common/ldes` binary on the SAME
 * model.json at `--seed 23000` gave two different sample paths -- 98.565592 on
 * a 22.04 host against 98.562490 under the containerized MATLAB's 24.04, which
 * is how this was found. The quantiles below therefore call
 * `line::fdlibm::log1p`, the algorithm `StrictMath.log1p` is DEFINED to use and
 * which `Math.log1p` was measured to match on every one of 400k draws; see
 * `line/util/fdlibm.h`. That fixes the engine to one arithmetic everywhere AND
 * puts it on the Java engine's, which is what "bit-compatible" was supposed to
 * mean. `pow`, `log`, `exp`, `sqrt`, `lgamma` and `tgamma` were measured to
 * agree across those glibcs and are left to libm.
 */
 
#include <cmath>
#include <string>
#include <vector>
 
#include "line/util/error.h"
#include "line/util/fdlibm.h"
 
namespace line {
namespace ldes {
namespace ssj {
 
// ---------------------------------------------------------------------------
// Closed-form quantiles
// ---------------------------------------------------------------------------
 
/** ExponentialDist.inverseF(lambda, u), SSJ's log1p spelling. */
inline double exponential_inverse(double lambda, double u) {
    if (lambda <= 0.0) throw InputError("ssj::exponential_inverse: lambda must be positive");
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
    if (u <= 0.0) return 0.0;
    return -line::fdlibm::log1p(-u) / lambda;
}
 
/** UniformDist.inverseF(a, b, u). */
inline double uniform_inverse(double a, double b, double u) {
    return a + (b - a) * u;
}
 
/** ParetoDist.inverseF(alpha, beta, u) = beta (1-u)^(-1/alpha). */
inline double pareto_inverse(double alpha, double beta, double u) {
    if (alpha <= 0.0 || beta <= 0.0)
        throw InputError("ssj::pareto_inverse: alpha and beta must be positive");
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
    return beta * std::pow(1.0 - u, -1.0 / alpha);
}
 
/**
 * WeibullDist.inverseF(alpha, lambda, delta, u).
 * SSJ parameterises the scale as lambda with the variate
 * delta + (1/lambda) (-log(1-u))^(1/alpha).
 */
inline double weibull_inverse(double alpha, double lambda, double delta, double u) {
    if (alpha <= 0.0 || lambda <= 0.0)
        throw InputError("ssj::weibull_inverse: alpha and lambda must be positive");
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
    if (u <= 0.0) return delta;
    return delta + std::pow(-line::fdlibm::log1p(-u), 1.0 / alpha) / lambda;
}
 
/** BernoulliDist.inverseF(p, u): 1 when u exceeds 1 - p. */
inline double bernoulli_inverse(double p, double u) {
    if (p < 0.0 || p > 1.0) throw InputError("ssj::bernoulli_inverse: p must lie in [0,1]");
    return (u <= 1.0 - p) ? 0.0 : 1.0;
}
 
// ---------------------------------------------------------------------------
// Normal, and Lognormal on top of it
// ---------------------------------------------------------------------------
 
/**
 * NormalDist.inverseF(0, 1, u) by Wichura's AS 241 (Applied Statistics 37,
 * 1988), the algorithm SSJ's `inverseF01` implements. Accurate to about 1e-16
 * over the whole range, which is what makes the Lognormal agree with SSJ to
 * the last few bits rather than only to a plotting tolerance.
 */
inline double normal_inverse01(double u) {
    if (u <= 0.0) return -std::numeric_limits<double>::infinity();
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
 
    const double q = u - 0.5;
    double r;
    if (std::fabs(q) <= 0.425) {
        r = 0.180625 - q * q;
        return q *
               (((((((2509.0809287301226727 * r + 33430.575583588128105) * r +
                     67265.770927008700853) * r + 45921.953931549871457) * r +
                   13731.693765509461125) * r + 1971.5909503065514427) * r +
                 133.14166789178437745) * r + 3.387132872796366608) /
               (((((((5226.495278852854561 * r + 28729.085735721942674) * r +
                     39307.89580009271061) * r + 21213.794301586595867) * r +
                   5394.1960214247511077) * r + 687.1870074920579083) * r +
                 42.313330701600911252) * r + 1.0);
    }
    r = (q < 0.0) ? u : 1.0 - u;
    r = std::sqrt(-std::log(r));
    double x;
    if (r <= 5.0) {
        r -= 1.6;
        x = (((((((7.7454501427834140764e-4 * r + 0.0227238449892691845833) * r +
                  0.24178072517745061177) * r + 1.27045825245236838258) * r +
                3.64784832476320460504) * r + 5.7694972214606914055) * r +
              4.6303378461565452959) * r + 1.42343711074968357734) /
            (((((((1.05075007164441684324e-9 * r + 5.475938084995344946e-4) * r +
                  0.0151986665636164571966) * r + 0.14810397642748007459) * r +
                0.68976733498510000455) * r + 1.6763848301838038494) * r +
              2.05319162663775882187) * r + 1.0);
    } else {
        r -= 5.0;
        x = (((((((2.01033439929228813265e-7 * r + 2.71155556874348757815e-5) * r +
                  0.0012426609473880784386) * r + 0.026532189526576123093) * r +
                0.29656057182850489123) * r + 1.7848265399172913358) * r +
              5.4637849111641143699) * r + 6.6579046435011037772) /
            (((((((2.04426310338993978564e-15 * r + 1.4215117583164458887e-7) * r +
                  1.8463183175100546818e-5) * r + 7.868691311456132591e-4) * r +
                0.0148753612908506148525) * r + 0.13692988092273580531) * r +
              0.59983220655588793769) * r + 1.0);
    }
    return (q < 0.0) ? -x : x;
}
 
/** NormalDist.inverseF(mu, sigma, u). */
inline double normal_inverse(double mu, double sigma, double u) {
    if (sigma <= 0.0) throw InputError("ssj::normal_inverse: sigma must be positive");
    return mu + sigma * normal_inverse01(u);
}
 
/** LognormalDist.inverseF(mu, sigma, u) = exp(mu + sigma Phi^-1(u)). */
inline double lognormal_inverse(double mu, double sigma, double u) {
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
    if (u <= 0.0) return 0.0;
    return std::exp(normal_inverse(mu, sigma, u));
}
 
// ---------------------------------------------------------------------------
// Gamma, and Erlang as the integer-shape case
// ---------------------------------------------------------------------------
 
namespace detail {
 
/** Regularized lower incomplete gamma P(a, x), by series and continued fraction. */
inline double gamma_p(double a, double x) {
    if (x <= 0.0) return 0.0;
    const double lg = std::lgamma(a);
    if (x < a + 1.0) {
        // Series expansion.
        double ap = a, sum = 1.0 / a, del = sum;
        for (int n = 0; n < 1000; ++n) {
            ap += 1.0;
            del *= x / ap;
            sum += del;
            if (std::fabs(del) < std::fabs(sum) * 1e-16) break;
        }
        return sum * std::exp(-x + a * std::log(x) - lg);
    }
    // Continued fraction for Q(a, x), then P = 1 - Q.
    const double tiny = 1e-300;
    double b = x + 1.0 - a, c = 1.0 / tiny, d = 1.0 / b, h = d;
    for (int i = 1; i <= 1000; ++i) {
        const double an = -static_cast<double>(i) * (static_cast<double>(i) - a);
        b += 2.0;
        d = an * d + b;
        if (std::fabs(d) < tiny) d = tiny;
        c = b + an / c;
        if (std::fabs(c) < tiny) c = tiny;
        d = 1.0 / d;
        const double del = d * c;
        h *= del;
        if (std::fabs(del - 1.0) < 1e-16) break;
    }
    return 1.0 - std::exp(-x + a * std::log(x) - lg) * h;
}
 
}  // namespace detail
 
/**
 * GammaDist.inverseF(alpha, lambda, u): the quantile of a Gamma of shape alpha
 * and RATE lambda, so the mean is alpha/lambda, which is SSJ's convention.
 *
 * Newton on P(alpha, lambda x) = u from a Wilson-Hilferty start, with a
 * bisection guard: the density vanishes at the origin for alpha > 1 and Newton
 * alone can step negative there. Converges to about 1e-14 relative, which is
 * below the resolution at which a service time can reorder two events.
 */
inline double gamma_inverse(double alpha, double lambda, double u) {
    if (alpha <= 0.0 || lambda <= 0.0)
        throw InputError("ssj::gamma_inverse: alpha and lambda must be positive");
    if (u <= 0.0) return 0.0;
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
 
    // Wilson-Hilferty: a cube-root normal approximation to the Gamma quantile.
    const double z = normal_inverse01(u);
    const double wh = alpha * std::pow(1.0 - 1.0 / (9.0 * alpha) + z / (3.0 * std::sqrt(alpha)), 3.0);
    double x = (wh > 0.0) ? wh : alpha * 0.5;
 
    double lo = 0.0, hi = 0.0;  // hi == 0 means "not yet bracketed"
    const double lg = std::lgamma(alpha);
    for (int it = 0; it < 200; ++it) {
        const double p = detail::gamma_p(alpha, x);
        if (p < u) {
            lo = x;
        } else {
            hi = x;
        }
        const double logpdf = (alpha - 1.0) * std::log(x) - x - lg;
        const double pdf = std::exp(logpdf);
        double step = (pdf > 0.0) ? (p - u) / pdf : 0.0;
        double next = x - step;
        if (!(next > lo) || (hi > 0.0 && !(next < hi)) || !(next > 0.0) || !std::isfinite(next)) {
            next = (hi > 0.0) ? 0.5 * (lo + hi) : 2.0 * x;
        }
        const double rel = std::fabs(next - x) / (std::fabs(x) + 1e-300);
        x = next;
        if (rel < 1e-15) break;
    }
    return x / lambda;
}
 
/**
 * ErlangGen(k, lambda): the Gamma of integer shape k and rate lambda. SSJ's
 * ErlangGen inverts the Gamma rather than summing k exponentials, which is why
 * it costs ONE uniform and not k -- the measurement above pinned that.
 */
inline double erlang_inverse(int k, double lambda, double u) {
    if (k <= 0) throw InputError("ssj::erlang_inverse: k must be positive");
    return gamma_inverse(static_cast<double>(k), lambda, u);
}
 
// ---------------------------------------------------------------------------
// Discrete families, by inversion
// ---------------------------------------------------------------------------
 
/**
 * PoissonDist.inverseF(lambda, u): the smallest k whose cdf reaches u.
 *
 * Summed forward from k = 0 with the pmf carried recursively. The result is an
 * INTEGER, so it agrees with SSJ exactly whenever the two cdfs put u on the
 * same side of a step, which no test has yet found a counterexample to.
 */
inline double poisson_inverse(double lambda, double u) {
    if (lambda <= 0.0) throw InputError("ssj::poisson_inverse: lambda must be positive");
    if (u <= 0.0) return 0.0;
    if (u >= 1.0) return std::numeric_limits<double>::infinity();
    double p = std::exp(-lambda), cdf = p;
    int k = 0;
    const int guard = static_cast<int>(lambda + 20.0 * std::sqrt(lambda) + 1000.0);
    while (cdf < u && k < guard) {
        ++k;
        p *= lambda / static_cast<double>(k);
        cdf += p;
    }
    return static_cast<double>(k);
}
 
/** BinomialDist.inverseF(n, p, u), by the same forward inversion. */
inline double binomial_inverse(int n, double p, double u) {
    if (n < 0) throw InputError("ssj::binomial_inverse: n must be non-negative");
    if (p < 0.0 || p > 1.0) throw InputError("ssj::binomial_inverse: p must lie in [0,1]");
    if (u <= 0.0) return 0.0;
    if (p <= 0.0) return 0.0;
    if (p >= 1.0) return static_cast<double>(n);
    double pk = std::pow(1.0 - p, static_cast<double>(n)), cdf = pk;
    int k = 0;
    while (cdf < u && k < n) {
        pk *= (static_cast<double>(n - k) / static_cast<double>(k + 1)) * (p / (1.0 - p));
        ++k;
        cdf += pk;
    }
    return static_cast<double>(k);
}
 
}  // namespace ssj
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_SSJ_VARIATES_H