ctmc_solve_reducible_blkdecomp.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_SOLVE_REDUCIBLE_BLKDECOMP_H
#define LINE_API_MC_CTMC_SOLVE_REDUCIBLE_BLKDECOMP_H
 
/**
 * @file
 * @ingroup api_mc
 * Limiting distribution of a reducible CTMC by direct block decomposition of
 * the generator.
 *
 * Templated port of matlab/src/api/mc/ctmc_solve_reducible_blkdecomp.m and
 * jar/src/main/java/jline/api/mc/Ctmc_solve_reducible_blkdecomp.java. Unlike
 * ctmc_solve_reducible it never uniformizes: the states are split into
 * transient and recurrent classes by strong connectivity, the expected sojourn
 * of the transient part is obtained from sojourn Q_tt = -p0_t (Q_tt is Hurwitz,
 * so this is a plain non-singular solve), the absorption probabilities follow
 * as hit = sojourn Q_ta + p0_r, and each recurrent class contributes its own
 * stationary vector scaled by the probability of reaching it.
 *
 * EXACT, AND THE EXACTNESS IS THE POINT. Every step is a finite linear solve
 * over the field of the rates: strong connectivity is combinatorial, Q_tt is
 * inverted once per starting class, and each recurrent class goes through
 * ctmc_solve. At Rational the absorption probabilities are exact rationals,
 * where the reference computes them in double precision on a matrix that is
 * ill-conditioned precisely when the transient class is nearly closed -- the
 * regime the routine exists to handle.
 *
 * The two thresholds MATLAB uses are structural tests on the input, not
 * convergence criteria, and are exposed as parameters: a reachability
 * probability below 1e-15 is treated as unreachable, and a state whose column
 * of |Q| sums below 1e-12 is treated as having no incoming rate.
 *
 * With an initial vector `pin` the answer is the limiting law reached FROM
 * `pin`, `absorb_limiting(pin)` as the reference computes it; `pis` and `pi0`
 * still report the uniform-start law of every component.
 */
 
#include <algorithm>
#include <cstddef>
#include <vector>
 
#include "line/api/mc/ctmc_solve.h"
#include "line/api/mc/stronglyconncomp.h"
#include "line/lang/lang_types.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 BlkDecompResult {
    std::vector<T> pi;             ///< limiting distribution, length N
    Matrix<T> pis;                 ///< numSCC x N, one row per starting component
    Matrix<T> pi0;                 ///< numSCC x N, the uniform-in-component starting vectors
    std::vector<std::size_t> scc;  ///< component index of each state, 1-based
    std::vector<bool> isrec;       ///< recurrence flag per component
};
 
/**
 * @brief Limiting distribution of a reducible CTMC by direct block
 *        decomposition of the generator.
 *
 * @param Qin generator; the diagonal is recomputed
 * @param pin initial distribution; empty if not available
 * @param reachTol probability below which a recurrent class counts as unreached
 * @param zeroColTol column-sum threshold for a state with no incoming rate
 */
template <class T>
BlkDecompResult<T> ctmc_solve_reducible_blkdecomp(const Matrix<T>& Qin, const std::vector<T>& pin,
                                                  double reachTol = 1e-15,
                                                  double zeroColTol = 1e-12) {
    const std::size_t N = Qin.rows();
    if (Qin.cols() != N) throw InputError("ctmc_solve_reducible_blkdecomp: generator is not square");
    if (!pin.empty() && pin.size() != N)
        throw InputError("ctmc_solve_reducible_blkdecomp: initial vector has the wrong length");
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
 
    const Matrix<T> Q = ctmc_makeinfgen(Qin);
 
    // Adjacency from the off-diagonal entries, `abs(Adj) > GlobalConstants.ArcTol`
    // as the reference tests it: an entry at or below the arc tolerance is not
    // an edge.
    const T arctol = num_traits<T>::from_double(lang::GlobalConstants::ArcTol);
    Matrix<T> Adj = Q;
    for (std::size_t i = 0; i < N; ++i) {
        Adj(i, i) = zero;
        for (std::size_t j = 0; j < N; ++j)
            if (i != j && !(num_abs(T(Adj(i, j))) > arctol)) Adj(i, j) = zero;
    }
    const SccResult s = stronglyconncomp(Adj);
    const std::size_t numSCC = s.numSCC();
 
    BlkDecompResult<T> r;
    r.scc = s.scc;
    r.isrec = s.recurrent;
 
    if (numSCC == 1) {
        r.pi = ctmc_solve(Q);
        r.pis = Matrix<T>(1, N);
        for (std::size_t j = 0; j < N; ++j) r.pis(0, j) = r.pi[j];
        r.pi0 = Matrix<T>();
        return r;
    }
 
    std::vector<std::size_t> transStates, recStates, transSccIds, recSccIds;
    for (std::size_t c = 0; c < numSCC; ++c) {
        if (s.recurrent[c]) {
            recSccIds.push_back(c);
            recStates.insert(recStates.end(), s.members[c].begin(), s.members[c].end());
        } else {
            transSccIds.push_back(c);
            transStates.insert(transStates.end(), s.members[c].begin(), s.members[c].end());
        }
    }
    std::sort(transStates.begin(), transStates.end());
    std::sort(recStates.begin(), recStates.end());
    const std::size_t nt = transStates.size(), nr = recStates.size();
    if (nr == 0)
        throw NumericError(ErrorKind{"Infeasible"},
            "ctmc_solve_reducible_blkdecomp: no recurrent class, the chain admits no limiting "
            "distribution");
 
    // Q_tt (transposed, ready for the solve) and Q_ta.
    Matrix<T> Q_ttT(nt, nt), Q_ta(nt, nr);
    std::vector<std::size_t> lupiv;
    Matrix<T> LU;
    bool lu_ready = false;
    if (nt > 0) {
        for (std::size_t a = 0; a < nt; ++a) {
            for (std::size_t b = 0; b < nt; ++b) Q_ttT(b, a) = Q(transStates[a], transStates[b]);
            for (std::size_t b = 0; b < nr; ++b) Q_ta(a, b) = Q(transStates[a], recStates[b]);
        }
    }
 
    const T rtol = num_traits<T>::from_double(reachTol);
    // Position of each state within recStates, for the per-class gather below.
    std::vector<std::size_t> recPos(N, static_cast<std::size_t>(-1));
    for (std::size_t k = 0; k < nr; ++k) recPos[recStates[k]] = k;
    // The stationary law of a recurrent class does not depend on the starting component, so it
    // is solved once and reused: re-solving it per component cost one dense LU per SCC.
    std::vector<std::vector<T>> pirec(numSCC);
 
    // The expected sojourn in the transient block from p0_t, the reference's
    // `sojourn = (-p0_t) / Q_tt`: above GMRES_MIN_STATES the Krylov route is
    // tried first (ctmc_gmres, then ctmc_bicgstab) and the direct solve is the
    // fallback, as in the reference; below it, and under exact arithmetic, the
    // LU of Q_tt' is factored once and reused for every starting vector.
    auto transient_sojourn = [&](const std::vector<T>& rhs) {
        std::vector<T> sojourn;
        if constexpr (num_traits<T>::has_transcendental) {
            if (nt > GMRES_MIN_STATES) {
                const GmresResult<T> g = ctmc_gmres(Q_ttT, rhs);
                if (g.flag == 0) {
                    sojourn = g.x;
                } else {
                    const BicgstabResult<T> bb = ctmc_bicgstab(Q_ttT, rhs);
                    if (bb.flag == 0) sojourn = bb.x;
                }
            }
        }
        if (sojourn.empty()) {
            if (!lu_ready) {
                LU = Q_ttT;
                lupiv = lu_factor(LU);
                lu_ready = true;
            }
            sojourn = rhs;
            lu_solve(LU, lupiv, sojourn);
        }
        return sojourn;
    };
 
    // The reference's absorb_limiting(p0): the limiting law reached from p0.
    auto absorb_limiting = [&](const std::vector<T>& p0) {
        std::vector<T> piv(N, zero);
        std::vector<T> hit(nr, zero);
        if (nt > 0) {
            bool anyT = false;
            std::vector<T> rhs(nt);
            for (std::size_t a = 0; a < nt; ++a) {
                rhs[a] = -p0[transStates[a]];
                if (rhs[a] != zero) anyT = true;
            }
            if (anyT) {
                // Solve sojourn * Q_tt = -p0_t, i.e. Q_tt' * sojourn' = -p0_t'.
                const std::vector<T> sojourn = transient_sojourn(rhs);
                for (std::size_t b = 0; b < nr; ++b) {
                    T acc = zero;
                    for (std::size_t a = 0; a < nt; ++a) acc += sojourn[a] * Q_ta(a, b);
                    hit[b] = acc;
                }
            }
        }
        for (std::size_t k = 0; k < nr; ++k) hit[k] += p0[recStates[k]];
 
        for (std::size_t cr : recSccIds) {
            const std::vector<std::size_t>& idx = s.members[cr];
            T reach = zero;
            for (std::size_t a : idx) reach += hit[recPos[a]];
            if (reach < rtol) continue;
            if (idx.size() == 1) {
                piv[idx[0]] = reach;
            } else {
                if (pirec[cr].empty()) pirec[cr] = ctmc_solve(detail::submatrix(Q, idx));
                const std::vector<T>& pi_c = pirec[cr];
                for (std::size_t k = 0; k < idx.size(); ++k) piv[idx[k]] = pi_c[k] * reach;
            }
        }
        return piv;
    };
 
    r.pis = Matrix<T>(numSCC, N, zero);
    r.pi0 = Matrix<T>(numSCC, N, zero);
    for (std::size_t c = 0; c < numSCC; ++c) {
        std::vector<T> p0(N, zero);
        const T w = one / num_traits<T>::from_int(static_cast<long>(s.members[c].size()));
        for (std::size_t a : s.members[c]) p0[a] = w;
        for (std::size_t j = 0; j < N; ++j) r.pi0(c, j) = p0[j];
        const std::vector<T> piv = absorb_limiting(p0);
        for (std::size_t j = 0; j < N; ++j) r.pis(c, j) = piv[j];
    }
 
    r.pi.assign(N, zero);
    if (pin.empty()) {
        // Probability of starting in each component: uniform over the components
        // that have some incoming rate.
        std::vector<T> pinl(numSCC, one);
        const T ztol = num_traits<T>::from_double(zeroColTol);
        for (std::size_t j = 0; j < N; ++j) {
            T cs = zero;
            for (std::size_t i = 0; i < N; ++i) cs += num_abs(T(Q(i, j)));
            if (cs < ztol) pinl[s.scc[j] - 1] = zero;
        }
        T tot = zero;
        for (const T& v : pinl) tot += v;
        if (tot > zero) {
            for (T& v : pinl) v /= tot;
        } else {
            for (std::size_t c = 0; c < numSCC; ++c)
                pinl[c] = one / num_traits<T>::from_int(static_cast<long>(numSCC));
        }
        for (std::size_t c = 0; c < numSCC; ++c) {
            if (!(pinl[c] > zero)) continue;
            for (std::size_t j = 0; j < N; ++j) r.pi[j] += r.pis(c, j) * pinl[c];
        }
        if (transSccIds.size() == 1)
            for (std::size_t j = 0; j < N; ++j) r.pi[j] = r.pis(transSccIds[0], j);
    } else {
        // The reference: the limiting law of the GIVEN initial vector, computed
        // directly as absorb_limiting(pin). Lumping pin onto the components and
        // mixing the uniform-start rows is NOT the same thing when pin puts its
        // mass unevenly inside a multi-state transient component.
        std::vector<T> p0 = pin;
        T tot0 = zero;
        for (const T& v : p0) tot0 += v;
        if (tot0 > zero)
            for (T& v : p0) v /= tot0;
        r.pi = absorb_limiting(p0);
    }
 
    T tot = zero;
    for (const T& v : r.pi) tot += v;
    if (tot > zero)
        for (T& v : r.pi) v /= tot;
    return r;
}
 
/** Overload without an initial vector. */
template <class T>
BlkDecompResult<T> ctmc_solve_reducible_blkdecomp(const Matrix<T>& Q, double reachTol = 1e-15,
                                                  double zeroColTol = 1e-12) {
    return ctmc_solve_reducible_blkdecomp(Q, std::vector<T>(), reachTol, zeroColTol);
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_SOLVE_REDUCIBLE_BLKDECOMP_H