ctmc_accumulated.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MC_CTMC_ACCUMULATED_H
#define LINE_API_MC_CTMC_ACCUMULATED_H
 
/**
 * @file
 * @ingroup api_mc
 * Accumulated state probabilities L(t) = int_0^t pi(s) ds of a CTMC by
 * uniformisation (Reibman and Trivedi 1989; Haverkort 1998, section 15.2.4).
 *
 * Templated port of matlab/src/api/mc/ctmc_accumulated.m, twin of
 * jline.api.mc.Ctmc_accumulated and line_solver.api.mc.ctmc.ctmc_accumulated.
 * With P = I + Q/q, q = 1.1 max|q_ii| and Poisson weights w_j at q h,
 *
 *   int_0^h pi0 exp(Q s) ds = (1/q) sum_j (1 - W_j) pi0 P^j,  W_j = sum_{i<=j} w_i
 *
 * so L(t) r is the expected reward accumulated over [0, t] for a rate reward
 * r, and L(t) r / t its time average. Successive time points are marched
 * interval by interval, so a grid costs one series per gap. The series is
 * truncated at the first K whose Poisson tail 1 - W_K is below tol and whose
 * neglected integral mass q h - sum_{j<=K}(1 - W_j) is below tol q h. An
 * interval with q h > 500 is split into equal segments, as in
 * ctmc_uniformization.h, and each segment end is rescaled to the mass of its
 * start so the truncation loss does not compound across segments.
 *
 * The weights involve exp(-q h): no exact instantiation.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
 
#include "line/api/mc/ctmc_uniformization.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
template <class T>
struct AccumulatedResult {
    Matrix<T> L;       ///< row k is int_0^{t_k} pi(s) ds
    Matrix<T> piT;     ///< row k is pi(t_k)
    std::size_t kmax;  ///< largest number of Poisson terms used by any interval
};
 
namespace detail {
 
template <class T>
void accumulate_interval(std::vector<T>& v, const Matrix<T>& P, const T& q, const T& h, double tol,
                         long maxiter, std::vector<T>& Lh, std::size_t& kmax) {
    using std::exp;
    using std::expm1;
    const std::size_t n = v.size();
    const T zero = num_traits<T>::from_int(0);
    Lh.assign(n, zero);
    kmax = 0;
    if (h == zero) return;
    if (q == zero) {
        for (std::size_t i = 0; i < n; ++i) Lh[i] = h * v[i];
        return;
    }
    const T qt = q * h;
    const double qtd = num_traits<T>::to_double(qt);
    const double MAXQT = 500.0;
    if (qtd > MAXQT) {
        const long nSeg = static_cast<long>(std::ceil(qtd / MAXQT));
        const T hs = h / num_traits<T>::from_int(nSeg);
        std::vector<T> dL;
        for (long s = 0; s < nSeg; ++s) {
            std::size_t K = 0;
            accumulate_interval(v, P, q, hs, tol, maxiter, dL, K);
            for (std::size_t i = 0; i < n; ++i) Lh[i] += dL[i];
            kmax = std::max(kmax, K);
        }
        return;
    }
    if (maxiter <= 0)
        maxiter = static_cast<long>(std::max(100.0, std::ceil(qtd + 10.0 * std::sqrt(qtd) + 20.0)));
    const T tolT = num_traits<T>::from_double(tol);
 
    // Truncation point.
    T w = exp(-qt);
    T tail = -expm1(-qt);  // 1 - W_0
    T S = tail;            // sum_{j<=K} (1 - W_j)
    long K = 0;
    while (K < maxiter && !(tail <= tolT && T(qt - S) <= T(tolT * qt))) {
        ++K;
        w = w * qt / num_traits<T>::from_int(K);
        tail = tail - w;
        if (tail < zero) tail = zero;
        S += tail;
    }
    kmax = static_cast<std::size_t>(K);
 
    // Series.
    w = exp(-qt);
    tail = -expm1(-qt);
    std::vector<T> Pk = v, vEnd(n), Lsum(n);
    for (std::size_t i = 0; i < n; ++i) {
        vEnd[i] = w * Pk[i];
        Lsum[i] = tail * Pk[i];
    }
    for (long j = 1; j <= K; ++j) {
        Pk = vecmat(Pk, P);
        w = w * qt / num_traits<T>::from_int(j);
        tail = tail - w;
        if (tail < zero) tail = zero;
        for (std::size_t i = 0; i < n; ++i) {
            vEnd[i] += w * Pk[i];
            Lsum[i] += tail * Pk[i];
        }
    }
    T sv = zero, se = zero;
    for (std::size_t i = 0; i < n; ++i) {
        sv += v[i];
        se += vEnd[i];
    }
    if (se != zero && sv != zero)
        for (T& e : vEnd) e = e * (sv / se);
    for (std::size_t i = 0; i < n; ++i) Lh[i] = Lsum[i] / q;
    v = vEnd;
}
 
}  // namespace detail
 
/**
 * @brief Accumulated state probabilities int_0^t pi(s) ds at each time point.
 *
 * @param pi0 initial distribution (row vector of length n)
 * @param Q   generator
 * @param t   nonempty, finite, nonnegative, nondecreasing time points
 * @param tol truncation tolerance (MATLAB default 1e-12)
 * @param maxiter Poisson term cap per interval, <= 0 for the MATLAB heuristic
 */
template <class T>
AccumulatedResult<T> ctmc_accumulated(const std::vector<T>& pi0, const Matrix<T>& Q, const std::vector<T>& t,
                                      double tol = 1e-12, long maxiter = -1) {
    static_assert(num_traits<T>::has_transcendental,
                  "ctmc_accumulated requires transcendental arithmetic: the Poisson weights involve "
                  "exp(-q t)");
    const std::size_t n = Q.rows();
    if (Q.cols() != n) throw InputError("ctmc_accumulated: generator is not square");
    if (pi0.size() != n) throw InputError("ctmc_accumulated: pi0 has the wrong length");
    const T zero = num_traits<T>::from_int(0);
    if (t.empty()) throw InputError("ctmc_accumulated: t must be nonempty");
    for (std::size_t k = 0; k < t.size(); ++k) {
        const double tk = num_traits<T>::to_double(t[k]);
        if (!std::isfinite(tk) || t[k] < zero || (k > 0 && t[k] < t[k - 1]))
            throw InputError("ctmc_accumulated: t must hold finite, nonnegative, nondecreasing time points");
    }
    const T q = detail::uniformization_rate(Q);
    const Matrix<T> P = (q == zero) ? Matrix<T>(n, n) : detail::uniformized_matrix(Q, q);
 
    AccumulatedResult<T> res{Matrix<T>(t.size(), n), Matrix<T>(t.size(), n), 0};
    std::vector<T> v = pi0, acc(n, zero), dL;
    T tprev = zero;
    for (std::size_t k = 0; k < t.size(); ++k) {
        std::size_t K = 0;
        detail::accumulate_interval(v, P, q, T(t[k] - tprev), tol, maxiter, dL, K);
        res.kmax = std::max(res.kmax, K);
        for (std::size_t i = 0; i < n; ++i) {
            acc[i] += dL[i];
            res.L(k, i) = acc[i];
            res.piT(k, i) = v[i];
        }
        tprev = t[k];
    }
    return res;
}
 
/** Single time point overload. */
template <class T>
AccumulatedResult<T> ctmc_accumulated(const std::vector<T>& pi0, const Matrix<T>& Q, const T& t,
                                      double tol = 1e-12, long maxiter = -1) {
    return ctmc_accumulated(pi0, Q, std::vector<T>{t}, tol, maxiter);
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_ACCUMULATED_H