ctmc_gth.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_GTH_H
#define LINE_API_MC_CTMC_GTH_H
 
/**
 * @file
 * @ingroup api_mc
 * Grassmann-Taksar-Heyman (GTH) elimination for the stationary vector, and the
 * limit matrix J = lim expm(Q t) of any finite chain built on it.
 *
 * Templated port of matlab/src/api/mc/ctmc_gth.m, dtmc_gth.m, ctmc_cesaro.m,
 * dtmc_cesaro.m and private/gth_irreducible.m; JAR twins jline.api.mc.Ctmc_gth
 * and Ctmc_cesaro, Python twin line_solver.api.mc.gth.
 *
 * GTH censors the states from last to first, a_ij += a_ik a_kj / s_k with the
 * pivot s_k = sum_{j<k} a_kj taken from the off-diagonal rates, never from the
 * diagonal, and back-substitutes x_k = sum_{i<k} x_i a_ik / s_k. There is no
 * subtraction, so every probability keeps small RELATIVE error, also on stiff
 * and nearly completely decomposable chains. Only the off-diagonal entries are
 * read, so dtmc_gth(P) = ctmc_gth(P) and periodic chains are fine. GTH uses only
 * +, * and /, so the exact instantiation returns the exact rational vector.
 * (Grassmann, Taksar, Heyman, Oper. Res. 33(5):1107-1116, 1985.)
 *
 * ctmc_cesaro gives every closed class C_k its GTH vector pi_k and every
 * transient state i the row J(i, C_k) = H(i,k) pi_k, with H = (-Q_TT)^{-1}
 * Q_{T,C} 1 the absorption probabilities. J is the matrix of maximum out-forests
 * of the transition digraph (Agaev and Chebotarev, Autom. Remote Control
 * 61(9):1424-1450, 2000): stochastic, idempotent, Q J = J Q = 0, and for a DTMC
 * the Cesaro limit of P^k whatever the periodicity.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/api/mc/stronglyconncomp.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
template <class T>
struct GthResult {
    std::vector<T> p;  ///< stationary vector (empty when flag == 1)
    int flag;          ///< 0 success, 1 several closed classes (not unique)
};
 
namespace detail {
 
/** Off-diagonal rates of a square matrix, validated nonnegative (and finite). */
template <class T>
Matrix<T> gth_offdiag(const Matrix<T>& Q, const char* fn) {
    const std::size_t n = Q.rows();
    if (Q.cols() != n) throw InputError(std::string(fn) + ": the generator must be a square matrix");
    const T zero = num_traits<T>::from_int(0);
    Matrix<T> A(n, n, zero);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            if (i == j) continue;
            const T v = Q(i, j);
            bool bad = v < zero;
            if constexpr (!num_traits<T>::is_exact) bad = bad || !std::isfinite(num_traits<T>::to_double(v));
            if (bad) throw InputError(std::string(fn) + ": the off-diagonal rates must be finite and nonnegative");
            A(i, j) = v;
        }
    return A;
}
 
/**
 * GTH on the off-diagonal rates of an irreducible chain (the diagonal is never
 * read). ok is false when a pivot is not positive (or not finite); A is a copy.
 */
template <class T>
std::vector<T> gth_irreducible(Matrix<T> A, bool& ok) {
    const std::size_t n = A.rows();
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
    ok = true;
    if (n == 1) return std::vector<T>{one};
    for (std::size_t i = 0; i < n; ++i) A(i, i) = zero;
    for (std::size_t k = n - 1; k >= 1; --k) {
        T s = zero;
        for (std::size_t j = 0; j < k; ++j) s += A(k, j);
        bool bad = !(s > zero);
        if constexpr (!num_traits<T>::is_exact) bad = bad || !std::isfinite(num_traits<T>::to_double(s));
        if (bad) {
            ok = false;
            return std::vector<T>();
        }
        for (std::size_t i = 0; i < k; ++i) A(i, k) /= s;
        for (std::size_t i = 0; i < k; ++i) {
            const T aik = A(i, k);
            if (aik == zero) continue;
            for (std::size_t j = 0; j < k; ++j) A(i, j) += aik * A(k, j);
        }
    }
    std::vector<T> x(n, zero);
    x[0] = one;
    for (std::size_t k = 1; k < n; ++k) {
        T v = zero;
        for (std::size_t i = 0; i < k; ++i) v += x[i] * A(i, k);
        x[k] = v;
    }
    T tot = zero;
    for (const T& v : x) tot += v;
    for (T& v : x) v /= tot;
    return x;
}
 
template <class T>
Matrix<T> gth_sub(const Matrix<T>& A, const std::vector<std::size_t>& r, const std::vector<std::size_t>& c) {
    Matrix<T> B(r.size(), c.size());
    for (std::size_t a = 0; a < r.size(); ++a)
        for (std::size_t b = 0; b < c.size(); ++b) B(a, b) = A(r[a], c[b]);
    return B;
}
 
}  // namespace detail
 
/**
 * @brief Stationary vector of a CTMC by GTH elimination, with a flag.
 *
 * Unichain generators are accepted (transient states get zero mass). With
 * several closed classes the vector is not unique: flag 1 and an empty p.
 */
template <class T>
GthResult<T> ctmc_gth_flag(const Matrix<T>& Q) {
    const Matrix<T> A = detail::gth_offdiag(Q, "ctmc_gth");
    const std::size_t n = A.rows();
    const T zero = num_traits<T>::from_int(0);
    if (n == 0) throw InputError("ctmc_gth: empty generator");
    if (n == 1) return GthResult<T>{std::vector<T>{num_traits<T>::from_int(1)}, 0};
    const SccResult s = stronglyconncomp(A);
    std::size_t nclosed = 0, closed = 0;
    for (std::size_t c = 0; c < s.recurrent.size(); ++c)
        if (s.recurrent[c]) {
            ++nclosed;
            closed = c;
        }
    if (nclosed != 1) return GthResult<T>{std::vector<T>(), 1};
    std::vector<std::size_t> S = s.members[closed];
    std::sort(S.begin(), S.end());
    bool ok = true;
    const std::vector<T> x = detail::gth_irreducible(detail::gth_sub(A, S, S), ok);
    if (!ok)
        throw NumericError(ErrorKind{"Degenerate"},
                           "ctmc_gth: GTH elimination met a zero pivot on a closed class: the rates underflowed");
    std::vector<T> p(n, zero);
    for (std::size_t k = 0; k < S.size(); ++k) p[S[k]] = x[k];
    return GthResult<T>{p, 0};
}
 
/** @brief Stationary vector by GTH; throws InputError on several closed classes. */
template <class T>
std::vector<T> ctmc_gth(const Matrix<T>& Q) {
    GthResult<T> r = ctmc_gth_flag(Q);
    if (r.flag != 0)
        throw InputError("ctmc_gth: the chain has several closed classes, so its stationary distribution is "
                         "not unique. Use ctmc_cesaro(Q) for the limit from every initial state, or "
                         "ctmc_solve_reducible(Q, pi0).");
    return r.p;
}
 
/** @brief Stationary vector of a DTMC by GTH; self-loops are not read. */
template <class T>
std::vector<T> dtmc_gth(const Matrix<T>& P) {
    return ctmc_gth(P);
}
 
/** @brief Limit matrix J = lim expm(Q t) of any finite CTMC (Agaev-Chebotarev 2000). */
template <class T>
Matrix<T> ctmc_cesaro(const Matrix<T>& Q) {
    const Matrix<T> A = detail::gth_offdiag(Q, "ctmc_cesaro");
    const std::size_t n = A.rows();
    const T zero = num_traits<T>::from_int(0);
    if (n == 0) throw InputError("ctmc_cesaro: empty generator");
    Matrix<T> J(n, n, zero);
    if (n == 1) {
        J(0, 0) = num_traits<T>::from_int(1);
        return J;
    }
    const SccResult s = stronglyconncomp(A);
    std::vector<std::size_t> closed, T_;
    for (std::size_t c = 0; c < s.recurrent.size(); ++c)
        if (s.recurrent[c]) closed.push_back(c);
    for (std::size_t v = 0; v < n; ++v)
        if (!s.recurrent[s.scc[v] - 1]) T_.push_back(v);
    const std::size_t K = closed.size(), nT = T_.size();
    std::vector<std::vector<T>> H(K, std::vector<T>(nT, zero));  // H[k][a]
    if (nT > 0) {
        Matrix<T> M(nT, nT, zero);
        for (std::size_t a = 0; a < nT; ++a) {
            T out = zero;
            for (std::size_t j = 0; j < n; ++j) out += A(T_[a], j);
            for (std::size_t b = 0; b < nT; ++b) M(a, b) = zero - A(T_[a], T_[b]);
            M(a, a) += out;
        }
        for (std::size_t k = 0; k < K; ++k) {
            std::vector<T> rhs(nT, zero);
            for (std::size_t a = 0; a < nT; ++a)
                for (std::size_t j : s.members[closed[k]]) rhs[a] += A(T_[a], j);
            H[k] = solve(M, rhs);
        }
        for (std::size_t a = 0; a < nT; ++a) {
            T tot = zero;
            for (std::size_t k = 0; k < K; ++k) {
                if (H[k][a] < zero) H[k][a] = zero;
                tot += H[k][a];
            }
            for (std::size_t k = 0; k < K; ++k) H[k][a] /= tot;
        }
    }
    for (std::size_t k = 0; k < K; ++k) {
        std::vector<std::size_t> Ck = s.members[closed[k]];
        std::sort(Ck.begin(), Ck.end());
        bool ok = true;
        const std::vector<T> pik = detail::gth_irreducible(detail::gth_sub(A, Ck, Ck), ok);
        if (!ok)
            throw NumericError(ErrorKind{"Degenerate"},
                               "ctmc_cesaro: GTH elimination met a zero pivot on a closed class: the rates underflowed");
        for (std::size_t i : Ck)
            for (std::size_t j = 0; j < Ck.size(); ++j) J(i, Ck[j]) = pik[j];
        for (std::size_t a = 0; a < nT; ++a)
            for (std::size_t j = 0; j < Ck.size(); ++j) J(T_[a], Ck[j]) = H[k][a] * pik[j];
    }
    return J;
}
 
/** @brief Cesaro limit of P^k for any finite DTMC; equals ctmc_cesaro(P - I). */
template <class T>
Matrix<T> dtmc_cesaro(const Matrix<T>& P) {
    return ctmc_cesaro(P);
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_GTH_H