me_sample.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MAM_ME_SAMPLE_H
#define LINE_API_MAM_ME_SAMPLE_H
 
/**
 * @file
 * @ingroup api_mam
 * Sample a matrix exponential by numerical inversion of its exact CDF.
 *
 * Templated port of matlab/lib/kpctoolbox/map/me_sample.m, cross-checked against
 * jar/src/main/java/jline/api/mam/Me_sample.java.
 *
 * AN ME IS A RENEWAL PROCESS, which is what separates this from `rap_sample`.
 * The inter-arrival times are i.i.d. draws from F(t) = 1 - pie exp(D0 t) e, so
 * the entry law is the SAME at every sample and the inversion table built from
 * it can be reused. `rap_sample` instead advances the entry law after each
 * arrival, which is correct for a RAP and wrong for an ME: it delivers a
 * correlated stream where the reference delivers an independent one.
 *
 * WHY A TABLE AND NOT A ROOT SOLVE. Routing ME through `rap_sample` costs a
 * full matrix exponential per bisection step, some fifty of them per variate;
 * on a 15-phase CME that is 140x the Java engine on the same model, enough for
 * the LDES wrapper to hit its wall-clock budget and return an empty result. The
 * reference tabulates the CDF once on a grid whose horizon is doubled until the
 * survival is below TAILTOL, locates the variate by binary search and polishes
 * it with Newton steps on the exact CDF and density, so the answer is not
 * limited by the linear interpolation of the table.
 *
 * ONE UNIFORM PER VARIATE, as in `rap_sample`, so a caller switching between the
 * two consumes the generator at the same rate.
 *
 * The grid itself is generated by a single expm(D0 h) and repeated
 * vector-matrix products, and the Newton step propagates the row vector by a
 * scaled Taylor series rather than re-exponentiating: both are O(K^2) where a
 * matrix exponential is O(K^3).
 *
 * ARITHMETIC: transcendental.
 */
 
#include <cmath>
#include <cstddef>
#include <vector>
 
#include "line/api/mam/map_moment.h"
#include "line/api/pfqn/pfqn_mc_common.h"
#include "line/num/number.h"
#include "line/util/eig.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mam {
 
namespace mesampledetail {
 
/** Number of tabulated grid points, including t = 0 (me_sample.m GRIDPTS). */
const std::size_t grid_points = 1000;
/** Cap on horizon doublings (me_sample.m MAXDBL). */
const int max_doublings = 40;
/** Survival mass left beyond the horizon (me_sample.m TAILTOL). */
const double tail_mass_tol = 1e-12;
/** Newton refinements applied to each inverted quantile (me_sample.m NEWTON). */
const int newton_steps = 3;
/** Largest ||A||*d a single Taylor substep handles. */
const double taylor_theta = 0.5;
/** Substep cap before falling back to a matrix exponential. */
const int taylor_max_steps = 64;
/** Taylor terms per substep. */
const int taylor_max_terms = 30;
 
/** w <- w M, the row-vector product the grid and the Taylor walk both need. */
template <class T>
std::vector<T> vec_times_mat(const std::vector<T>& w, const Matrix<T>& M) {
    const std::size_t n = w.size();
    std::vector<T> out(n, num_traits<T>::from_int(0));
    for (std::size_t j = 0; j < n; ++j)
        for (std::size_t i = 0; i < n; ++i) out[j] += w[i] * M(i, j);
    return out;
}
 
/**
 * w exp(A d) by a scaled Taylor series on the VECTOR.
 *
 * Falls back to a matrix exponential when the substep count would exceed the
 * cap, which is the only regime where the series is the more expensive of the
 * two.
 */
template <class T>
std::vector<T> expm_propagate(const std::vector<T>& w, const Matrix<T>& A, double norm_a, double d) {
    const std::size_t n = w.size();
    if (!(d > 0.0)) return w;
    int steps = 1;
    const double theta = norm_a * d;
    if (theta > taylor_theta) {
        steps = static_cast<int>(std::ceil(theta / taylor_theta));
        if (steps > taylor_max_steps)
            return vec_times_mat(w, expm(A, num_traits<T>::from_double(d)));
    }
    const double ds = d / steps;
    std::vector<T> acc = w;
    for (int s = 0; s < steps; ++s) {
        std::vector<T> term = acc, next = acc;
        for (int k = 1; k <= taylor_max_terms; ++k) {
            term = vec_times_mat(term, A);
            const T c = num_traits<T>::from_double(ds / k);
            double max_term = 0.0, max_acc = 0.0;
            for (std::size_t j = 0; j < n; ++j) {
                term[j] = T(term[j] * c);
                next[j] += term[j];
                const double at = std::fabs(num_traits<T>::to_double(term[j]));
                if (at > max_term) max_term = at;
                const double an = std::fabs(num_traits<T>::to_double(next[j]));
                if (an > max_acc) max_acc = an;
            }
            if (max_term <= 1e-18 * (max_acc > 1e-300 ? max_acc : 1e-300)) break;
        }
        acc = next;
    }
    return acc;
}
 
/**
 * The dominant (least negative) eigenvalue of D0, which governs the tail decay.
 *
 * The spectrum is taken in double precision whatever T is: the extrapolation it
 * feeds is a tail approximation, so its own last digits do not matter, and the
 * fallback -1/mean is what the reference uses when the decomposition fails.
 */
template <class T>
double dominant_rate(const Matrix<T>& A, double mean) {
    double eta = -std::numeric_limits<double>::infinity();
    const std::size_t n = A.rows();
    Matrix<double> Ad(n, n);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) Ad(i, j) = num_traits<T>::to_double(A(i, j));
    try {
        const std::vector<std::complex<double> > ev = eig_values(Ad);
        for (std::size_t i = 0; i < ev.size(); ++i)
            if (ev[i].real() > eta) eta = ev[i].real();
    } catch (const std::exception&) {
        eta = -std::numeric_limits<double>::infinity();
    }
    if (!(eta < 0.0) || !std::isfinite(eta)) eta = (mean > 0.0) ? -1.0 / mean : -1.0;
    return eta;
}
 
}  // namespace mesampledetail
 
/**
 * Stateful ME sampler holding the inversion table.
 *
 * The table is what makes the LDES usage pattern -- one variate per service
 * event, from a generator that lives as long as the station does -- cost a
 * binary search rather than a matrix exponential.
 */
template <class T>
class MeSampler {
public:
    MeSampler() : n_(0), exponential_(false), exp_rate_(0.0), h_(0.0), t_end_(0.0), s_end_(0.0),
                  eta_(-1.0), norm_a_(0.0) {}
 
    /** Build the table of `m`, whose D0 is the ME matrix and whose entry law is map_pie. */
    explicit MeSampler(const Map<T>& m, const std::vector<T>& a0 = std::vector<T>()) {
        static_assert(num_traits<T>::has_transcendental,
                      "MeSampler inverts a matrix-exponential distribution function");
        n_ = m.order();
        if (n_ == 0) throw InputError("me_sample: the representation is empty");
        A_ = m.D0;
        std::vector<T> alpha = a0.empty() ? map_pie(m) : a0;
        if (alpha.size() != n_) throw InputError("me_sample: the entry law has the wrong length");
 
        row_sum_.assign(n_, num_traits<T>::from_int(0));
        norm_a_ = 0.0;
        for (std::size_t i = 0; i < n_; ++i) {
            double abs_row = 0.0;
            for (std::size_t j = 0; j < n_; ++j) {
                row_sum_[i] += A_(i, j);
                abs_row += std::fabs(num_traits<T>::to_double(A_(i, j)));
            }
            if (abs_row > norm_a_) norm_a_ = abs_row;
        }
 
        // A one-phase ME with alpha = 1 IS an exponential; inverting it in
        // closed form skips a thousand-point table no caller would gain from.
        exponential_ = (n_ == 1) && std::fabs(num_traits<T>::to_double(alpha[0]) - 1.0) < 1e-12 &&
                       num_traits<T>::to_double(A_(0, 0)) < 0.0;
        exp_rate_ = exponential_ ? -num_traits<T>::to_double(A_(0, 0)) : 0.0;
        if (exponential_) {
            h_ = t_end_ = s_end_ = 0.0;
            eta_ = -exp_rate_;
            return;
        }
 
        const double mean = num_traits<T>::to_double(map_mean(m));
        const double var = num_traits<T>::to_double(map_var(m));
        const double sigma = (var > 0.0) ? std::sqrt(var) : 0.0;
        double horizon = mean + 10.0 * sigma;
        if (!(horizon > 0.0) || !std::isfinite(horizon)) horizon = 1.0;
        for (int k = 0; k < mesampledetail::max_doublings; ++k) {
            if (survival_by_expm(alpha, horizon) < mesampledetail::tail_mass_tol) break;
            horizon *= 2.0;
        }
 
        const std::size_t g = mesampledetail::grid_points;
        h_ = horizon / static_cast<double>(g - 1);
        t_.resize(g);
        F_.resize(g);
        w_.resize(g);
        // ONE matrix exponential for the whole grid: exp(D0 (k+1)h) applied to
        // the entry law is exp(D0 kh) applied to it and then multiplied by
        // exp(D0 h), so the table is a walk of row-vector products.
        const Matrix<T> Eh = expm(A_, num_traits<T>::from_double(h_));
        w_[0] = alpha;
        t_[0] = 0.0;
        F_[0] = clamp_unit(1.0 - sum_of(alpha));
        for (std::size_t i = 1; i < g; ++i) {
            w_[i] = mesampledetail::vec_times_mat(w_[i - 1], Eh);
            t_[i] = static_cast<double>(i) * h_;
            // Forced nondecreasing: roundoff breaks monotonicity in the tail,
            // where consecutive values differ by less than eps.
            double c = clamp_unit(1.0 - sum_of(w_[i]));
            if (c < F_[i - 1]) c = F_[i - 1];
            F_[i] = c;
        }
        t_end_ = t_[g - 1];
        const double surv = 1.0 - F_[g - 1];
        s_end_ = (surv > 0.0) ? surv : 0.0;
        eta_ = mesampledetail::dominant_rate(A_, mean);
    }
 
    /** One variate, consuming exactly one uniform. */
    double next(pfqn::McRng& rng) const {
        const double u = pfqn::mc_uniform01(rng);
        return quantile(u);
    }
 
    /** The inverse CDF at `u`, exposed so a caller can supply its own uniform. */
    double quantile(double u) const {
        if (exponential_) return -std::log(1.0 - u) / exp_rate_;
        const std::size_t g = F_.size();
        if (u <= F_[0]) return 0.0;
        if (u >= F_[g - 1]) {
            // Exponential tail S(x) ~ S(tEnd) exp(eta (x - tEnd)). Clamping to
            // the grid endpoint instead truncates the tail and biases the mean.
            const double tail = 1.0 - u;
            if (s_end_ <= 0.0 || tail <= 0.0 || !(eta_ < 0.0)) return t_end_;
            const double x = t_end_ + std::log(s_end_ / tail) / (-eta_);
            return (x > t_end_) ? x : t_end_;
        }
        std::size_t lo = 0, hi = g - 1;
        while (hi - lo > 1) {
            const std::size_t mid = (lo + hi) / 2;
            if (F_[mid] <= u)
                lo = mid;
            else
                hi = mid;
        }
        const double den = F_[lo + 1] - F_[lo];
        double x = (den > 0.0) ? t_[lo] + (u - F_[lo]) / den * h_ : t_[lo];
        const double left = t_[lo], right = t_[lo] + h_;
 
        for (int k = 0; k < mesampledetail::newton_steps; ++k) {
            const std::vector<T> w =
                mesampledetail::expm_propagate(w_[lo], A_, norm_a_, x - left);
            const double surv = sum_of(w);
            double f = 0.0;
            for (std::size_t j = 0; j < n_; ++j)
                f -= num_traits<T>::to_double(w[j]) * num_traits<T>::to_double(row_sum_[j]);
            if (!(f > 0.0)) break;
            const double err = (1.0 - surv) - u;
            if (std::fabs(err) < 1e-14) break;
            const double xn = x - err / f;
            if (!(xn > left) || !(xn < right)) break;
            const bool converged =
                std::fabs(xn - x) <= 1e-15 * (std::fabs(x) > 1.0 ? std::fabs(x) : 1.0);
            x = xn;
            if (converged) break;
        }
        return x;
    }
 
private:
    double survival_by_expm(const std::vector<T>& alpha, double t) const {
        const std::vector<T> w =
            mesampledetail::vec_times_mat(alpha, expm(A_, num_traits<T>::from_double(t)));
        const double s = sum_of(w);
        return (s > 0.0) ? s : 0.0;
    }
 
    static double sum_of(const std::vector<T>& v) {
        double s = 0.0;
        for (std::size_t i = 0; i < v.size(); ++i) s += num_traits<T>::to_double(v[i]);
        return s;
    }
 
    static double clamp_unit(double x) { return x < 0.0 ? 0.0 : (x > 1.0 ? 1.0 : x); }
 
    std::size_t n_;
    Matrix<T> A_;
    std::vector<T> row_sum_;
    bool exponential_;
    double exp_rate_;
    std::vector<double> t_;
    std::vector<double> F_;
    std::vector<std::vector<T> > w_;
    double h_;
    double t_end_;
    double s_end_;
    double eta_;
    double norm_a_;
};
 
/**
 * `me_sample`: n INDEPENDENT variates of the matrix exponential (D0, D1).
 *
 * Passing a MAP or a RAP here samples its stationary marginal independently,
 * which is a deliberate renewal approximation: the autocorrelation is dropped.
 * Use `rap_sample` to retain it.
 *
 * @param m   the process, read through its entry law map_pie and its D0
 * @param n   number of inter-arrival times to draw
 * @param rng generator, advanced by one uniform per variate
 * @param a0  entry law; empty selects map_pie
 */
template <class T>
std::vector<T> me_sample(const Map<T>& m, std::size_t n, pfqn::McRng& rng,
                         const std::vector<T>& a0 = std::vector<T>()) {
    const MeSampler<T> s(m, a0);
    std::vector<T> out;
    out.reserve(n);
    for (std::size_t i = 0; i < n; ++i)
        out.push_back(num_traits<T>::from_double(s.next(rng)));
    return out;
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_API_MAM_ME_SAMPLE_H