ctmc_ssd.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_SSD_H
#define LINE_API_MC_CTMC_SSD_H
 
/**
 * @file
 * @ingroup api_mc
 * Stochastically stable distribution of a perturbed Markov chain (Wicks 2009).
 *
 * Port of matlab/src/api/mc/ctmc_ssd.m and dtmc_ssd.m; JAR twin
 * jline.api.mc.Ctmc_ssd, Python twin line_solver.api.mc.gth.ctmc_ssd. The
 * off-diagonal rates behave as Q_eps(i,j) ~ C(i,j) eps^R(i,j) as eps -> 0+
 * (rows are sources; a zero rate is R = Inf or C = 0), and the result is the
 * exact limit of the stationary vector, no eps ever substituted. Q_eps must be
 * irreducible for eps > 0. The recursion (see ctmc_ssd.m for the full account):
 *   1. shift R to minimum 0 and take G0, the zero-resistance graph;
 *   2. one closed class in G0: GTH on the constant chain over G0;
 *   3. else a class of G0 with two or more states: drop its positive-resistance
 *      internal edges (Lemma 7.22) and censor all its states but one, with the
 *      constant N0 = (-Q0_ss)^{-1} (Theorem 7.31), in (C,R) arithmetic;
 *   4. else divide the rates out of the closed singletons by eps^smin (Lemma
 *      7.16) and keep the scaled solution on those states.
 * (C,R) arithmetic: a sum keeps the least resistance and adds the tied costs, a
 * product adds resistances and multiplies costs. Resistances are compared with
 * a relative tolerance 1e-10, so the computation is double-valued.
 * (J. R. Wicks, PhD thesis, Brown University, 2009, Ch. 7.)
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/mc/ctmc_gth.h"
#include "line/api/mc/stronglyconncomp.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
namespace detail {
 
using SsdMat = Matrix<double>;
 
inline void ssd_mult(const SsdMat& C1, const SsdMat& R1, const SsdMat& C2, const SsdMat& R2, double tol,
                     SsdMat& Cc, SsdMat& Rc) {
    const double inf = std::numeric_limits<double>::infinity();
    const std::size_t a = C1.rows(), b = C1.cols(), c = C2.cols();
    Cc = SsdMat(a, c, 0.0);
    Rc = SsdMat(a, c, inf);
    for (std::size_t t = 0; t < b; ++t)
        for (std::size_t i = 0; i < a; ++i)
            for (std::size_t j = 0; j < c; ++j) {
                const double r = R1(i, t) + R2(t, j);
                const double cc = C1(i, t) * C2(t, j);
                if (!std::isfinite(r) || !(cc > 0.0)) continue;
                if (r < Rc(i, j) - tol) {
                    Cc(i, j) = cc;
                    Rc(i, j) = r;
                } else if (std::fabs(r - Rc(i, j)) <= tol) {
                    Cc(i, j) += cc;
                }
            }
}
 
inline void ssd_add_inplace(SsdMat& C1, SsdMat& R1, const SsdMat& C2, const SsdMat& R2, double tol) {
    for (std::size_t i = 0; i < C1.rows(); ++i)
        for (std::size_t j = 0; j < C1.cols(); ++j) {
            if (!std::isfinite(R2(i, j))) continue;
            if (R2(i, j) < R1(i, j) - tol) {
                R1(i, j) = R2(i, j);
                C1(i, j) = C2(i, j);
            } else if (std::fabs(R2(i, j) - R1(i, j)) <= tol) {
                C1(i, j) += C2(i, j);
            }
        }
}
 
inline std::vector<double> ssd_recurse(SsdMat C, SsdMat R, double tol, int depth) {
    const double inf = std::numeric_limits<double>::infinity();
    const std::size_t n = C.rows();
    if (n == 1) return std::vector<double>{1.0};
    if (depth <= 0)
        throw NumericError(ErrorKind{"NonConvergence"},
                           "ctmc_ssd: the SSD recursion did not terminate; please report the input matrices");
    double rmin = inf;
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            if (std::isfinite(R(i, j)) && R(i, j) < rmin) rmin = R(i, j);
    SsdMat Z(n, n, 0.0);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            if (std::isfinite(R(i, j))) {
                R(i, j) -= rmin;
                if (R(i, j) <= tol) {
                    R(i, j) = 0.0;
                    Z(i, j) = 1.0;
                }
            }
    const SccResult s = stronglyconncomp(Z);
    std::vector<std::size_t> closed;
    for (std::size_t c = 0; c < s.recurrent.size(); ++c)
        if (s.recurrent[c]) closed.push_back(c);
    if (closed.size() == 1) {
        std::vector<std::size_t> S = s.members[closed[0]];
        std::sort(S.begin(), S.end());
        SsdMat A(S.size(), S.size(), 0.0);
        for (std::size_t a = 0; a < S.size(); ++a)
            for (std::size_t b = 0; b < S.size(); ++b) A(a, b) = C(S[a], S[b]) * Z(S[a], S[b]);
        bool ok = true;
        const std::vector<double> x = gth_irreducible(A, ok);
        std::vector<double> p(n, 0.0);
        for (std::size_t k = 0; k < S.size(); ++k) p[S[k]] = x[k];
        return p;
    }
    // Components are numbered by decreasing size, so component 0 is the largest.
    if (s.members[0].size() > 1) {
        std::vector<std::size_t> S = s.members[0];
        std::sort(S.begin(), S.end());
        std::vector<bool> inS(n, false);
        for (std::size_t v : S) inS[v] = true;
        for (std::size_t i : S)
            for (std::size_t j : S)
                if (i != j && Z(i, j) == 0.0) {  // Lemma 7.22
                    C(i, j) = 0.0;
                    R(i, j) = inf;
                }
        const std::vector<std::size_t> sv(S.begin() + 1, S.end());
        std::vector<std::size_t> keep;
        for (std::size_t v = 0; v < n; ++v)
            if (!inS[v] || v == S[0]) keep.push_back(v);
        const std::size_t ns = sv.size(), m = keep.size();
        SsdMat C0(n, n, 0.0);
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) C0(i, j) = C(i, j) * Z(i, j);
        SsdMat M(ns, ns, 0.0);  // -Q0_ss
        for (std::size_t a = 0; a < ns; ++a) {
            double out = 0.0;
            for (std::size_t j = 0; j < n; ++j) out += C0(sv[a], j);
            for (std::size_t b = 0; b < ns; ++b) M(a, b) = -C0(sv[a], sv[b]);
            M(a, a) += out;
        }
        SsdMat N0(ns, ns, 0.0), RN(ns, ns, inf);
        for (std::size_t b = 0; b < ns; ++b) {
            std::vector<double> e(ns, 0.0);
            e[b] = 1.0;
            const std::vector<double> col = solve(M, e);
            for (std::size_t a = 0; a < ns; ++a) {
                N0(a, b) = col[a] > 0.0 ? col[a] : 0.0;
                if (N0(a, b) > 0.0) RN(a, b) = 0.0;
            }
        }
        auto sub = [](const SsdMat& X, const std::vector<std::size_t>& r, const std::vector<std::size_t>& c) {
            return gth_sub(X, r, c);
        };
        SsdMat Cks, Rks, Cv, Rv;
        ssd_mult(sub(C, keep, sv), sub(R, keep, sv), N0, RN, tol, Cks, Rks);
        ssd_mult(Cks, Rks, sub(C, sv, keep), sub(R, sv, keep), tol, Cv, Rv);
        SsdMat Cr = sub(C, keep, keep), Rr = sub(R, keep, keep);
        ssd_add_inplace(Cr, Rr, Cv, Rv, tol);
        for (std::size_t i = 0; i < m; ++i) {
            Cr(i, i) = 0.0;
            Rr(i, i) = inf;
        }
        const std::vector<double> x = ssd_recurse(Cr, Rr, tol, depth - 1);
        std::vector<double> p(n, 0.0);
        for (std::size_t k = 0; k < m; ++k) p[keep[k]] = x[k];
        std::vector<double> y(ns, 0.0);  // x * C0(keep, s)
        for (std::size_t a = 0; a < ns; ++a)
            for (std::size_t k = 0; k < m; ++k) y[a] += x[k] * C0(keep[k], sv[a]);
        for (std::size_t b = 0; b < ns; ++b) {
            double v = 0.0;
            for (std::size_t a = 0; a < ns; ++a) v += y[a] * N0(a, b);
            p[sv[b]] = v;
        }
        double tot = 0.0;
        for (double v : p) tot += v;
        for (double& v : p) v /= tot;
        return p;
    }
    // Every class of G0 is a single state and several are absorbing in G0.
    std::vector<bool> isClosed(n, false);
    for (std::size_t v = 0; v < n; ++v) isClosed[v] = s.recurrent[s.scc[v] - 1];
    double smin = inf;
    for (std::size_t i = 0; i < n; ++i)
        if (isClosed[i])
            for (std::size_t j = 0; j < n; ++j)
                if (std::isfinite(R(i, j)) && R(i, j) < smin) smin = R(i, j);
    for (std::size_t i = 0; i < n; ++i)
        if (isClosed[i])
            for (std::size_t j = 0; j < n; ++j)
                if (std::isfinite(R(i, j))) R(i, j) -= smin;
    std::vector<double> p = ssd_recurse(C, R, tol, depth - 1);
    double tot = 0.0;
    for (std::size_t v = 0; v < n; ++v) {
        if (!isClosed[v]) p[v] = 0.0;
        tot += p[v];
    }
    if (!(tot > 0.0))
        throw NumericError(ErrorKind{"Degenerate"},
                           "ctmc_ssd: internal error: the scaled chain put no mass on the closed states");
    for (double& v : p) v /= tot;
    return p;
}
 
}  // namespace detail
 
/** @brief Stochastically stable distribution, rates C(i,j) eps^R(i,j) (Wicks 2009). */
inline std::vector<double> ctmc_ssd(const Matrix<double>& Cin, const Matrix<double>& Rin) {
    const double inf = std::numeric_limits<double>::infinity();
    const std::size_t n = Cin.rows();
    if (Cin.cols() != n || Rin.rows() != n || Rin.cols() != n)
        throw InputError("ctmc_ssd: C and R must be square matrices of the same size");
    if (n == 0) throw InputError("ctmc_ssd: empty input");
    Matrix<double> C = Cin, R = Rin, adj(n, n, 0.0);
    double rmax = 1.0;
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            if (i == j) {
                C(i, j) = 0.0;
                R(i, j) = inf;
                continue;
            }
            const double c = C(i, j), r = R(i, j);
            if (std::isnan(c) || std::isnan(r) || c < 0.0 || std::isinf(c))
                throw InputError("ctmc_ssd: costs must be finite and nonnegative, and neither C nor R may contain NaN");
            if (r == -inf) throw InputError("ctmc_ssd: a resistance of -Inf is not a valid perturbation order");
            if (c == 0.0 || std::isinf(r)) {
                C(i, j) = 0.0;
                R(i, j) = inf;
            } else {
                adj(i, j) = 1.0;
                rmax = std::max(rmax, std::fabs(r));
            }
        }
    if (n == 1) return std::vector<double>{1.0};
    const SccResult s = stronglyconncomp(adj);
    std::size_t nclosed = 0;
    for (bool b : s.recurrent) nclosed += b ? 1 : 0;
    if (nclosed != 1)
        throw InputError("ctmc_ssd: the perturbed chain has " + std::to_string(nclosed) +
                         " closed classes for eps > 0, so it has no unique stationary distribution and no "
                         "stochastically stable distribution");
    return detail::ssd_recurse(C, R, 1e-10 * rmax, static_cast<int>(4 * n + 8));
}
 
/** @brief SSD of a perturbed transition matrix P_eps(i,j) ~ C(i,j) eps^R(i,j); equals ctmc_ssd. */
inline std::vector<double> dtmc_ssd(const Matrix<double>& C, const Matrix<double>& R) {
    return ctmc_ssd(C, R);
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_SSD_H