nnpp.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_LANG_NNPP_H
#define LINE_LANG_NNPP_H
 
/**
 * @file
 * @ingroup line_lang
 * NNPP: a nonstationary non-Poisson arrival process, generated by thinning a
 * majorizing stream (Lewis and Shedler 1979) whose epochs come from a Weibull
 * equilibrium renewal process inverted through the majorizing mean-value
 * function (SCIATA, Liu, Kuhl, Liu and Wilson, WSC 2015).
 *
 * The rate function lambda(t) lives on [0, S] and is either an EPTMP
 * (exponential of a polynomial plus sines) or a piecewise-linear interpolation
 * of knots. With dispersion C = 1 the renewal gaps are Exp(1) and the method
 * IS Lewis-Shedler thinning of an NHPP; with C != 1 the Weibull shape is chosen
 * so the renewal gaps have mean 1 and SCV C, which is what gives the arrival
 * counts their dispersion ratio.
 *
 * EVERY DETERMINISTIC QUANTITY HERE IS PINNED ACROSS THE FOUR CODEBASES at 1e-12
 * relative: the majorizer (golden section with a safety max, or the exact knot
 * maximum), the target mean value mu(t) (exact trapezoid, or composite Simpson
 * with a fixed panel rule), and the Weibull shape/scale (exactly 200 bisection
 * steps in log-shape). The operation order follows the shared specification,
 * so do not "simplify" an expression here without changing the other three.
 *
 * The sampler is stochastic and only statistically comparable across
 * languages. One `Stream` per arrival path carries the last renewal epoch in
 * majorizer space, the first-gap flag and the last accepted arrival.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iostream>
#include <limits>
#include <random>
#include <string>
#include <vector>
 
#include "line/solvers/ldes/ldes_ssj_variates.h"
#include "line/util/error.h"
#include "line/util/line_console.h"
 
namespace line {
namespace lang {
 
class Nnpp {
public:
    /** The rate-function family, with the codes `sn.proc` carries. */
    enum Kind { EPTMP = 1, PWLINEAR = 2 };
 
    /** One arrival path: renewal epoch in majorizer space, first-gap flag, last arrival. */
    struct Stream {
        double u = 0.0;
        bool first = true;
        double t_last = 0.0;
    };
 
    Nnpp() {}
 
    /**
     * EPTMP: lambda(t) = exp(sum_i alpha_i t^i + sum_k gamma_k sin(omega_k t + phi_k))
     * on [0, horizon].
     */
    static Nnpp eptmp(const std::vector<double>& alpha, const std::vector<double>& gamma,
                      const std::vector<double>& omega, const std::vector<double>& phi,
                      double horizon, double dispersion = 1.0, bool cyclic = true,
                      int Q = 200) {
        if (alpha.empty()) throw InputError("NNPP.eptmp: alpha must hold at least one coefficient");
        if (gamma.size() != omega.size() || gamma.size() != phi.size())
            throw InputError("NNPP.eptmp: gamma, omega and phi must have the same length");
        if (!(horizon > 0.0) || !std::isfinite(horizon))
            throw InputError("NNPP.eptmp: the horizon must be finite and positive");
        for (double v : alpha)
            if (!std::isfinite(v)) throw InputError("NNPP.eptmp: alpha must be finite");
        for (std::size_t k = 0; k < gamma.size(); ++k)
            if (!std::isfinite(gamma[k]) || !std::isfinite(omega[k]) || !std::isfinite(phi[k]))
                throw InputError("NNPP.eptmp: gamma, omega and phi must be finite");
        Nnpp n;
        n.kind_ = EPTMP;
        n.alpha_ = alpha;
        n.gamma_ = gamma;
        n.omega_ = omega;
        n.phi_ = phi;
        n.horizon_ = horizon;
        n.init(dispersion, cyclic, Q);
        return n;
    }
 
    /** PWLINEAR: lambda linearly interpolated between (times[j], rates[j]); horizon = times.back(). */
    static Nnpp piecewise_linear(const std::vector<double>& times, const std::vector<double>& rates,
                                 double dispersion = 1.0, bool cyclic = true, int Q = 200) {
        if (times.size() < 2 || rates.size() != times.size())
            throw InputError("NNPP.piecewiseLinear: times and rates must have the same length, at "
                             "least 2");
        if (times[0] != 0.0) throw InputError("NNPP.piecewiseLinear: times[0] must be 0");
        for (std::size_t j = 0; j + 1 < times.size(); ++j)
            if (!(times[j + 1] > times[j]))
                throw InputError("NNPP.piecewiseLinear: times must be strictly increasing");
        if (!std::isfinite(times.back()))
            throw InputError("NNPP.piecewiseLinear: times must be finite");
        for (double r : rates)
            if (!(r >= 0.0) || !std::isfinite(r))
                throw InputError("NNPP.piecewiseLinear: rates must be finite and non-negative");
        Nnpp n;
        n.kind_ = PWLINEAR;
        n.times_ = times;
        n.rates_ = rates;
        n.horizon_ = times.back();
        n.init(dispersion, cyclic, Q);
        return n;
    }
 
    /**
     * Rebuild from the `params` row of `sn.proc` ([m+1, p, alpha, gamma, omega, phi]
     * for EPTMP, [n, times, rates] for PWLINEAR).
     */
    static Nnpp from_params(int kind, const std::vector<double>& params, double horizon,
                            double dispersion, bool cyclic, int Q) {
        if (kind == EPTMP) {
            if (params.size() < 2) throw InputError("NNPP: malformed EPTMP params row");
            const std::size_t m1 = static_cast<std::size_t>(params[0]);
            const std::size_t p = static_cast<std::size_t>(params[1]);
            if (params.size() != 2 + m1 + 3 * p) throw InputError("NNPP: malformed EPTMP params row");
            std::vector<double> a(params.begin() + 2, params.begin() + 2 + m1);
            std::vector<double> g(params.begin() + 2 + m1, params.begin() + 2 + m1 + p);
            std::vector<double> w(params.begin() + 2 + m1 + p, params.begin() + 2 + m1 + 2 * p);
            std::vector<double> f(params.begin() + 2 + m1 + 2 * p, params.end());
            return eptmp(a, g, w, f, horizon, dispersion, cyclic, Q);
        }
        if (kind == PWLINEAR) {
            if (params.empty()) throw InputError("NNPP: malformed PWLINEAR params row");
            const std::size_t n = static_cast<std::size_t>(params[0]);
            if (params.size() != 1 + 2 * n) throw InputError("NNPP: malformed PWLINEAR params row");
            std::vector<double> t(params.begin() + 1, params.begin() + 1 + n);
            std::vector<double> r(params.begin() + 1 + n, params.end());
            return piecewise_linear(t, r, dispersion, cyclic, Q);
        }
        throw InputError("NNPP: unknown rate-function kind " + std::to_string(kind));
    }
 
    // ---------------------------------------------------------------------
    // Rate function and mean value
    // ---------------------------------------------------------------------
 
    /** lambda(t), cyclic aware; zero past the horizon when not cyclic. */
    double rate_at(double t) const {
        if (t < 0.0) return 0.0;
        if (cyclic_) return lam(reduce(t));
        if (t > horizon_) return 0.0;
        return lam(t);
    }
 
    /** mu(t) = int_0^t lambda, cyclic aware; constant past the horizon when not cyclic. */
    double mean_value(double t) const {
        if (!(t > 0.0)) return 0.0;
        if (cyclic_) {
            const double k = std::floor(t / horizon_);
            const double r = t - k * horizon_;
            return k * mu_s_ + mu_within(r);
        }
        if (t >= horizon_) return mu_s_;
        return mu_within(t);
    }
 
    double time_average_rate() const { return mu_s_ / horizon_; }
    double mean() const {
        return mu_s_ > 0.0 ? horizon_ / mu_s_ : std::numeric_limits<double>::infinity();
    }
    double rate() const { return mu_s_ / horizon_; }
    double scv() const { return std::numeric_limits<double>::quiet_NaN(); }
    double skewness() const { return std::numeric_limits<double>::quiet_NaN(); }
 
    /** The Q+1 majorizer breakpoints z_i = i*S/Q. */
    const std::vector<double>& majorizer_breakpoints() const { return z_; }
    /** The Q piecewise-constant majorizer rates. */
    const std::vector<double>& majorizer_rates() const { return lt_; }
    /** muTilde(S), the majorizer mass over one cycle. */
    double mean_value_majorizer() const { return M_.back(); }
    /** muTilde(t), cyclic aware. */
    double mean_value_majorizer(double t) const {
        if (!(t > 0.0)) return 0.0;
        double k = 0.0;
        if (cyclic_) {
            k = std::floor(t / horizon_);
            t = t - k * horizon_;
        } else if (t >= horizon_) {
            return M_.back();
        }
        std::size_t i = static_cast<std::size_t>(std::floor(t / h_));
        if (i >= lt_.size()) i = lt_.size() - 1;
        return k * M_.back() + M_[i] + lt_[i] * (t - z_[i]);
    }
    double acceptance_rate() const { return M_.back() > 0.0 ? mu_s_ / M_.back() : 0.0; }
    double dispersion() const { return disp_; }
    bool is_cyclic() const { return cyclic_; }
    double horizon() const { return horizon_; }
    double weibull_shape() const { return wshape_; }
    double weibull_scale() const { return wscale_; }
    int q() const { return static_cast<int>(lt_.size()); }
    int kind() const { return kind_; }
    std::string kind_name() const { return kind_ == EPTMP ? "eptmp" : "pwlinear"; }
    std::string method() const { return disp_ == 1.0 ? "thinning" : "sciata"; }
    std::string name() const { return "NNPP"; }
 
    const std::vector<double>& alpha() const { return alpha_; }
    const std::vector<double>& gamma() const { return gamma_; }
    const std::vector<double>& omega() const { return omega_; }
    const std::vector<double>& phi() const { return phi_; }
    const std::vector<double>& times() const { return times_; }
    const std::vector<double>& rates() const { return rates_; }
 
    /** The `params` row of `sn.proc`, see `from_params`. */
    std::vector<double> params() const {
        std::vector<double> p;
        if (kind_ == EPTMP) {
            p.push_back(static_cast<double>(alpha_.size()));
            p.push_back(static_cast<double>(gamma_.size()));
            p.insert(p.end(), alpha_.begin(), alpha_.end());
            p.insert(p.end(), gamma_.begin(), gamma_.end());
            p.insert(p.end(), omega_.begin(), omega_.end());
            p.insert(p.end(), phi_.begin(), phi_.end());
        } else {
            p.push_back(static_cast<double>(times_.size()));
            p.insert(p.end(), times_.begin(), times_.end());
            p.insert(p.end(), rates_.begin(), rates_.end());
        }
        return p;
    }
 
    /**
     * Weibull (shape, scale) whose mean is 1 and SCV is C: 200 bisection steps on
     * x = log(shape) in [log 0.02, log 200]. C == 1 is exactly (1, 1).
     */
    static void weibull_for_dispersion(double C, double& shape, double& scale) {
        if (!(C > 0.0) || !std::isfinite(C))
            throw InputError("NNPP: the dispersion ratio must be finite and positive");
        if (C == 1.0) {
            shape = 1.0;
            scale = 1.0;
            return;
        }
        auto g = [C](double x) {
            const double a = std::exp(x);
            return std::exp(std::lgamma(1.0 + 2.0 / a) - 2.0 * std::lgamma(1.0 + 1.0 / a)) - 1.0 - C;
        };
        double lo = std::log(0.02), hi = std::log(200.0);
        const double glo = g(lo), ghi = g(hi);
        if (!(glo > 0.0 && ghi < 0.0))
            throw InputError("NNPP: no Weibull shape in [0.02, 200] gives dispersion " +
                             std::to_string(C));
        for (int it = 0; it < 200; ++it) {
            const double mid = 0.5 * (lo + hi);
            if (g(mid) > 0.0)
                lo = mid;
            else
                hi = mid;
        }
        shape = std::exp(0.5 * (lo + hi));
        scale = std::exp(-std::lgamma(1.0 + 1.0 / shape));
    }
 
    // ---------------------------------------------------------------------
    // Majorizer inverse and sampling
    // ---------------------------------------------------------------------
 
    /** Majorizer rate in force at an in-cycle time tc. */
    double majorizer_at(double tc) const {
        std::size_t i = static_cast<std::size_t>(std::floor(tc / h_));
        if (i >= lt_.size()) i = lt_.size() - 1;
        return lt_[i];
    }
 
    /** muTilde^{-1}(u), absolute time, cyclic aware; +Inf when the majorizer has no mass. */
    double majorizer_inverse(double u) const {
        const double MQ = M_.back();
        if (!(MQ > 0.0)) return std::numeric_limits<double>::infinity();
        double k = 0.0;
        if (cyclic_) {
            k = std::floor(u / MQ);
            u = u - k * MQ;
        } else if (u > MQ) {
            return std::numeric_limits<double>::infinity();
        }
        const std::size_t Q = lt_.size();
        // smallest i in 1..Q with M[i] >= u
        std::size_t i = static_cast<std::size_t>(std::lower_bound(M_.begin() + 1, M_.end(), u) -
                                                 M_.begin());
        if (i > Q) i = Q;
        while (i < Q && !(lt_[i - 1] > 0.0)) ++i;
        while (i > 1 && !(lt_[i - 1] > 0.0)) --i;  // u == MQ landed past the last positive rate
        const double t = z_[i - 1] + (u - M_[i - 1]) / lt_[i - 1];
        return k * horizon_ + t;
    }
 
    /** The equilibrium (first) renewal gap: scale * Y^(1/shape), Y ~ Gamma(1/shape, 1). */
    template <class Unif>
    double equilibrium_gap(Unif& unif) const {
        const double u = unif();
        if (wshape_ == 1.0) return -std::log1p(-u);
        const double y = ldes::ssj::gamma_inverse(1.0 / wshape_, 1.0, u);
        return wscale_ * std::pow(y, 1.0 / wshape_);
    }
 
    /** A noninitial renewal gap: scale * (-log(1-U))^(1/shape). */
    template <class Unif>
    double renewal_gap(Unif& unif) const {
        const double u = unif();
        const double e = -std::log1p(-u);
        if (wshape_ == 1.0) return e;
        return wscale_ * std::pow(e, 1.0 / wshape_);
    }
 
    /**
     * The next accepted arrival epoch (absolute time) of the path `s`, or +Inf when
     * a non-cyclic process has no more arrivals. `unif()` returns U(0,1).
     */
    template <class Unif>
    double next_arrival(Stream& s, Unif& unif) const {
        if (!(M_.back() > 0.0)) return std::numeric_limits<double>::infinity();
        for (;;) {
            const double X = s.first ? equilibrium_gap(unif) : renewal_gap(unif);
            s.first = false;
            s.u = s.u + X;
            if (!cyclic_ && s.u > M_.back()) return std::numeric_limits<double>::infinity();
            const double t = majorizer_inverse(s.u);
            const double tc = cyclic_ ? reduce(t) : std::min(t, horizon_);
            const double maj = majorizer_at(tc);
            double p = (maj > 0.0) ? lam(tc) / maj : 0.0;
            if (p < 0.0) p = 0.0;
            if (p > 1.0) p = 1.0;
            if (unif() <= p) {
                s.t_last = t;
                return t;
            }
        }
    }
 
    /**
     * n successive interarrival times along one path from time 0. A non-cyclic
     * process that runs out of arrivals returns fewer than n.
     */
    template <class Unif>
    std::vector<double> sample(std::size_t n, Unif&& unif) const {
        std::vector<double> out;
        out.reserve(n);
        Stream s;
        double prev = 0.0;
        for (std::size_t j = 0; j < n; ++j) {
            const double t = next_arrival(s, unif);
            if (!std::isfinite(t)) break;
            out.push_back(t - prev);
            prev = t;
        }
        return out;
    }
 
    /** sample(n) off a std::mt19937_64, U drawn in (0,1). */
    std::vector<double> sample(std::size_t n, std::mt19937_64& rng) const {
        std::uniform_real_distribution<double> U(0.0, 1.0);
        auto unif = [&]() {
            double u = U(rng);
            while (!(u > 0.0)) u = U(rng);
            return u;
        };
        return sample(n, unif);
    }
 
private:
    void init(double dispersion, bool cyclic, int Q) {
        if (Q < 1) throw InputError("NNPP: the majorizer needs Q >= 1 subintervals");
        weibull_for_dispersion(dispersion, wshape_, wscale_);
        if (dispersion > 2.0 && util::LineConsole::get_verbose() != util::VerboseLevel::SILENT)
            std::cerr << "[LINE] Warning: NNPP: dispersion " << dispersion
                      << " exceeds 2; SCIATA is designed for C in (0, 2]." << std::endl;
        disp_ = dispersion;
        cyclic_ = cyclic;
        const double S = horizon_;
        h_ = S / Q;
        z_.assign(static_cast<std::size_t>(Q) + 1, 0.0);
        for (int i = 0; i <= Q; ++i) z_[static_cast<std::size_t>(i)] = i * S / Q;
        lt_.assign(static_cast<std::size_t>(Q), 0.0);
        for (int i = 1; i <= Q; ++i) {
            const double a = z_[static_cast<std::size_t>(i - 1)], b = z_[static_cast<std::size_t>(i)];
            lt_[static_cast<std::size_t>(i - 1)] = (kind_ == EPTMP) ? golden_max(a, b) : pwl_max(a, b);
        }
        M_.assign(static_cast<std::size_t>(Q) + 1, 0.0);
        for (int i = 1; i <= Q; ++i) {
            const std::size_t ii = static_cast<std::size_t>(i);
            M_[ii] = M_[ii - 1] + lt_[ii - 1] * (z_[ii] - z_[ii - 1]);
        }
        mu_s_ = (kind_ == EPTMP) ? simpson(S, 64 * static_cast<long>(Q)) : trapezoid(S);
    }
 
    /** Reduce an absolute time into [0, S). */
    double reduce(double t) const {
        const double k = std::floor(t / horizon_);
        double r = t - k * horizon_;
        if (r < 0.0) r = 0.0;
        if (r >= horizon_) r = 0.0;
        return r;
    }
 
    /** lambda at an in-cycle time t in [0, S]. */
    double lam(double t) const {
        if (kind_ == EPTMP) {
            const std::size_t m = alpha_.size() - 1;
            double s = alpha_[m];
            for (std::size_t i = m; i-- > 0;) s = s * t + alpha_[i];
            for (std::size_t k = 0; k < gamma_.size(); ++k)
                s += gamma_[k] * std::sin(omega_[k] * t + phi_[k]);
            return std::exp(s);
        }
        const std::size_t n = times_.size();
        // largest j <= n-2 with times[j] <= t
        std::size_t j = static_cast<std::size_t>(
            std::upper_bound(times_.begin(), times_.begin() + static_cast<long>(n - 1), t) -
            times_.begin());
        j = (j == 0) ? 0 : j - 1;
        return rates_[j] + (rates_[j + 1] - rates_[j]) * (t - times_[j]) / (times_[j + 1] - times_[j]);
    }
 
    double golden_max(double a, double b) const {
        const double phi = (1.0 + std::sqrt(5.0)) / 2.0;
        const double delta = 1e-4;
        double best = std::max(lam(a), lam(b));
        double x1 = a + (2.0 - phi) * (b - a);
        double x2 = a + (phi - 1.0) * (b - a);
        double f1 = lam(x1), f2 = lam(x2);
        best = std::max(best, std::max(f1, f2));
        while ((b - a) >= delta) {
            if (f1 < f2) {
                a = x1;
                x1 = x2;
                f1 = f2;
                x2 = a + (phi - 1.0) * (b - a);
                f2 = lam(x2);
                best = std::max(best, f2);
            } else {
                b = x2;
                x2 = x1;
                f2 = f1;
                x1 = a + (2.0 - phi) * (b - a);
                f1 = lam(x1);
                best = std::max(best, f1);
            }
        }
        return best;
    }
 
    double pwl_max(double a, double b) const {
        double best = std::max(lam(a), lam(b));
        for (std::size_t j = 0; j < times_.size(); ++j)
            if (a < times_[j] && times_[j] < b) best = std::max(best, rates_[j]);
        return best;
    }
 
    /** Composite Simpson of lambda on [0, t] with N (even) panels. */
    double simpson(double t, long N) const {
        const double h = t / static_cast<double>(N);
        double s = lam(0.0) + lam(t);
        for (long i = 1; i < N; ++i) s += ((i % 2 == 1) ? 4.0 : 2.0) * lam(static_cast<double>(i) * h);
        return s * h / 3.0;
    }
 
    /** Exact trapezoid of the piecewise-linear lambda on [0, t], t in [0, S]. */
    double trapezoid(double t) const {
        double s = 0.0;
        for (std::size_t j = 0; j + 1 < times_.size(); ++j) {
            const double a = times_[j];
            if (!(t > a)) break;
            const double b = std::min(t, times_[j + 1]);
            const double fb = (b == times_[j + 1]) ? rates_[j + 1] : lam(b);
            s += 0.5 * (rates_[j] + fb) * (b - a);
        }
        return s;
    }
 
    /** mu(t) for t in [0, S]. */
    double mu_within(double t) const {
        if (!(t > 0.0)) return 0.0;
        if (kind_ == PWLINEAR) return trapezoid(t);
        if (t == horizon_) return mu_s_;
        const long Q = static_cast<long>(lt_.size());
        long N = 2 * static_cast<long>(std::ceil(32.0 * Q * t / horizon_));
        if (N < 2) N = 2;
        return simpson(t, N);
    }
 
    int kind_ = EPTMP;
    std::vector<double> alpha_, gamma_, omega_, phi_;
    std::vector<double> times_, rates_;
    double horizon_ = 1.0;
    double disp_ = 1.0;
    bool cyclic_ = true;
    double wshape_ = 1.0, wscale_ = 1.0;
    double h_ = 1.0;
    std::vector<double> z_, lt_, M_;
    double mu_s_ = 0.0;
};
 
}  // namespace lang
}  // namespace line
 
#endif  // LINE_LANG_NNPP_H