ldqbd_censor.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_MAM_LDQBD_CENSOR_H
#define LINE_API_MAM_LDQBD_CENSOR_H
 
/**
 * @file
 * @ingroup api_mam
 * Stationary distribution of a chain of connected LDQBDs by successive
 * censoring.
 *
 * Templated port of matlab/src/api/mam/ldqbd_censor.m, implementing the
 * algorithm of N. Baer, A. Al Hanbali, R. J. Boucherie, J. C. W. van Ommeren,
 * "A successive censoring algorithm for a system of connected
 * LDQBD-processes", Annals of Operations Research 310:389-410, 2022.
 *
 * The state space is split into S sets, each a level-dependent QBD on its
 * own; set s has levels 1..L_s and level l of set s has phases[s][l] phases.
 * States are ordered by (set, level, phase). Transitions between sets are
 * arbitrary. Sets 1..S-1 are censored out in turn,
 *
 *     Q^{k+1}_{ij} = Q^k_{ij} + Q^k_{ik} (-Q^k_{kk})^{-1} Q^k_{kj},   i, j > k,
 *
 * the last LDQBD is solved by level reduction (Gaver), and the vector is
 * expanded back with pi_k = (sum_{i>k} pi_i Q^k_{ik}) (-Q^k_{kk})^{-1}. Each
 * transient LDQBD is factored once by block-tridiagonal LU (block Thomas,
 * S_l = D_l - W_l S_{l-1}^{-1} U_{l-1}); its inverse is never formed.
 *
 * STRUCTURE. The method applies exactly when every censored diagonal block
 * Q^k_{kk} is block tridiagonal by level. The update is a sum of nonnegative
 * terms, so the numerical pattern equals the structural one and the check is
 * made directly on Q^k_{kk}; a failure raises LINE:unsupported:Structure
 * (another order of the sets may succeed, Remark 2 of the paper).
 *
 * ARITHMETIC. The algorithm is finite (no tolerance, no iteration), so it
 * instantiates at Rational and returns the exact stationary law. Exact T
 * raises NumericError on an exactly singular Schur complement; inexact T also
 * raises when the 1-norm reciprocal condition of a block is below machine
 * epsilon, matching the MATLAB rcond test. The row-sum check is exact at
 * Rational and relative (tol) otherwise.
 */
 
#include <cstddef>
#include <map>
#include <string>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mam {
 
/** Stationary distribution of a chain of connected LDQBDs. */
template <class T>
struct LdqbdCensorResult {
    std::vector<T> pi;                                 ///< in the ordering of Q
    std::vector<std::vector<std::vector<T>>> pi_cells;  ///< pi_cells[s][l], level l of set s
};
 
namespace ldqbd_censor_detail {
 
template <class T>
T zero() { return num_traits<T>::from_int(0); }
 
template <class T>
double norm1(const Matrix<T>& A) {
    double best = 0.0;
    for (std::size_t j = 0; j < A.cols(); ++j) {
        double s = 0.0;
        for (std::size_t i = 0; i < A.rows(); ++i) s += num_traits<T>::to_double(num_abs(T(A(i, j))));
        if (s > best) best = s;
    }
    return best;
}
 
/** Inverse of a Schur complement, refusing a singular one. */
template <class T>
Matrix<T> checked_inverse(const Matrix<T>& D, std::size_t level, std::size_t set) {
    const std::string msg = "ldqbd_censor: the Schur complement of level " + std::to_string(level + 1) +
                            " of set " + std::to_string(set + 1) + " is singular; the chain may be reducible.";
    Matrix<T> inv;
    try {
        inv = inverse(D);
    } catch (const NumericError&) {
        throw NumericError(ErrorKind{"Singular"}, msg);
    }
    if (!num_traits<T>::is_exact) {
        const double rc = 1.0 / (norm1(D) * norm1(inv));
        if (!(rc >= 2.220446049250313e-16)) throw NumericError(ErrorKind{"Singular"}, msg);
    }
    return inv;
}
 
template <class T>
Matrix<T> block(const Matrix<T>& M, std::size_t r0, std::size_t nr, std::size_t c0, std::size_t nc) {
    Matrix<T> B(nr, nc, zero<T>());
    for (std::size_t i = 0; i < nr; ++i)
        for (std::size_t j = 0; j < nc; ++j) B(i, j) = M(r0 + i, c0 + j);
    return B;
}
 
/** Block-tridiagonal LU of M = -Q_kk: S on the diagonal, G below, U above. */
template <class T>
struct Factor {
    std::vector<std::size_t> P, off;
    std::vector<Matrix<T>> S, Sinv, G, U;
 
    Factor(const Matrix<T>& M, const std::vector<std::size_t>& P_, std::size_t set, bool last_singular)
        : P(P_), off(P_.size() + 1, 0) {
        const std::size_t L = P.size();
        for (std::size_t l = 0; l < L; ++l) off[l + 1] = off[l] + P[l];
        S.resize(L);
        Sinv.resize(L);
        G.resize(L);
        U.resize(L);
        for (std::size_t l = 0; l < L; ++l) {
            Matrix<T> D = block(M, off[l], P[l], off[l], P[l]);
            if (l > 0) {
                const Matrix<T> W = block(M, off[l], P[l], off[l - 1], P[l - 1]);
                G[l] = matmul(W, Sinv[l - 1]);
                const Matrix<T> GU = matmul(G[l], U[l - 1]);
                for (std::size_t i = 0; i < D.rows(); ++i)
                    for (std::size_t j = 0; j < D.cols(); ++j) D(i, j) -= GU(i, j);
            }
            if (l + 1 < L) U[l] = block(M, off[l], P[l], off[l + 1], P[l + 1]);
            S[l] = D;
            if (l + 1 < L || !last_singular) Sinv[l] = checked_inverse(D, l, set);
        }
    }
 
    std::size_t levels() const { return P.size(); }
 
    /** X = M^{-1} Y. */
    Matrix<T> colsolve(const Matrix<T>& Y) const {
        const std::size_t L = levels(), m = Y.cols();
        std::vector<Matrix<T>> Z(L), X(L);
        for (std::size_t l = 0; l < L; ++l) {
            Z[l] = block(Y, off[l], P[l], 0, m);
            if (l > 0) {
                const Matrix<T> GZ = matmul(G[l], Z[l - 1]);
                for (std::size_t i = 0; i < P[l]; ++i)
                    for (std::size_t j = 0; j < m; ++j) Z[l](i, j) -= GZ(i, j);
            }
        }
        X[L - 1] = matmul(Sinv[L - 1], Z[L - 1]);
        for (std::size_t l = L - 1; l-- > 0;) {
            const Matrix<T> UX = matmul(U[l], X[l + 1]);
            for (std::size_t i = 0; i < P[l]; ++i)
                for (std::size_t j = 0; j < m; ++j) Z[l](i, j) -= UX(i, j);
            X[l] = matmul(Sinv[l], Z[l]);
        }
        Matrix<T> out(Y.rows(), m, zero<T>());
        for (std::size_t l = 0; l < L; ++l)
            for (std::size_t i = 0; i < P[l]; ++i)
                for (std::size_t j = 0; j < m; ++j) out(off[l] + i, j) = X[l](i, j);
        return out;
    }
 
    /** x = y M^{-1}. */
    std::vector<T> rowsolve(const std::vector<T>& y) const {
        const std::size_t L = levels();
        std::vector<std::vector<T>> w(L);
        for (std::size_t l = 0; l < L; ++l) {
            std::vector<T> r(y.begin() + off[l], y.begin() + off[l + 1]);
            if (l > 0) {
                const std::vector<T> wu = vecmul(w[l - 1], U[l - 1]);
                for (std::size_t i = 0; i < r.size(); ++i) r[i] -= wu[i];
            }
            w[l] = vecmul(r, Sinv[l]);
        }
        std::vector<T> x(y.size(), zero<T>());
        for (std::size_t i = 0; i < P[L - 1]; ++i) x[off[L - 1] + i] = w[L - 1][i];
        for (std::size_t l = L - 1; l-- > 0;) {
            const std::vector<T> xn(x.begin() + off[l + 1], x.begin() + off[l + 2]);
            const std::vector<T> xg = vecmul(xn, G[l + 1]);
            for (std::size_t i = 0; i < P[l]; ++i) x[off[l] + i] = w[l][i] - xg[i];
        }
        return x;
    }
};
 
/** Throws LINE:unsupported:Structure if M is not block tridiagonal by level. */
template <class T>
void check(const Matrix<T>& M, const std::vector<std::size_t>& P, std::size_t set) {
    std::vector<std::size_t> lev;
    for (std::size_t l = 0; l < P.size(); ++l)
        for (std::size_t i = 0; i < P[l]; ++i) lev.push_back(l + 1);
    for (std::size_t i = 0; i < M.rows(); ++i)
        for (std::size_t j = 0; j < M.cols(); ++j) {
            const std::size_t d = lev[i] > lev[j] ? lev[i] - lev[j] : lev[j] - lev[i];
            if (d > 1 && !(M(i, j) == zero<T>()))
                throw UnsupportedError(
                    ErrorKind{"Structure"},
                    "ldqbd_censor: set " + std::to_string(set + 1) +
                        " is not block tridiagonal by level after censoring out the sets before it "
                        "(transition from level " + std::to_string(lev[i]) + " to level " +
                        std::to_string(lev[j]) + "); successive censoring does not apply in this set order.");
        }
}
 
}  // namespace ldqbd_censor_detail
 
/**
 * Stationary distribution of a chain of connected LDQBDs (ldqbd_censor.m).
 *
 * @param Q      generator ordered by (set, level, phase)
 * @param phases phases[s][l] is the number of phases of level l of set s
 * @param tol    relative tolerance on the row sums (ignored at exact T,
 *               where the rows must sum to exactly zero)
 */
template <class T>
LdqbdCensorResult<T> ldqbd_censor(const Matrix<T>& Q, const std::vector<std::vector<std::size_t>>& phases,
                                  double tol = 1e-10) {
    using namespace ldqbd_censor_detail;
    if (phases.empty()) throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: phases must give one vector of phase counts per set.");
    const std::size_t S = phases.size();
    std::vector<std::size_t> nset(S, 0), setOff(S + 1, 0);
    for (std::size_t s = 0; s < S; ++s) {
        if (phases[s].empty())
            throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: phases{" + std::to_string(s + 1) +
                                                        "} must be a non-empty vector of positive integers.");
        for (std::size_t p : phases[s]) {
            if (p < 1)
                throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: phases{" + std::to_string(s + 1) +
                                                            "} must be a non-empty vector of positive integers.");
            nset[s] += p;
        }
        setOff[s + 1] = setOff[s] + nset[s];
    }
    const std::size_t n = setOff[S];
    if (Q.rows() != Q.cols()) throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: Q must be square.");
    if (Q.rows() != n)
        throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: Q has order " + std::to_string(Q.rows()) +
                                                    " but the phase counts add up to " + std::to_string(n) + ".");
 
    // sparse row store of the working generator, indexed by global state
    std::vector<std::map<std::size_t, T>> rows(n);
    double scale = 1.0;
    for (std::size_t i = 0; i < n; ++i) {
        const double d = num_traits<T>::to_double(num_abs(T(Q(i, i))));
        if (d > scale) scale = d;
    }
    for (std::size_t i = 0; i < n; ++i) {
        T rs = zero<T>();
        for (std::size_t j = 0; j < n; ++j) {
            const T v = Q(i, j);
            if (v == zero<T>()) continue;
            if (i != j && v < zero<T>())
                throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: the off-diagonal entries of Q must be nonnegative.");
            rows[i][j] = v;
            rs += v;
        }
        const bool bad = num_traits<T>::is_exact ? !(rs == zero<T>())
                                                 : num_traits<T>::to_double(num_abs(rs)) > tol * scale;
        if (bad) throw InputError(ErrorKind{"Argument"}, "ldqbd_censor: the rows of Q must sum to zero.");
    }
 
    // successive censoring of sets 1..S-1
    std::vector<Factor<T>> fac;
    std::vector<std::vector<std::size_t>> bRows(S);
    std::vector<std::vector<std::vector<T>>> bVals(S);
    for (std::size_t k = 0; k + 1 < S; ++k) {
        const std::size_t o0 = setOff[k], o1 = setOff[k + 1], nk = nset[k];
        Matrix<T> M(nk, nk, zero<T>());
        std::map<std::size_t, std::size_t> colPos;
        for (std::size_t i = 0; i < nk; ++i)
            for (const auto& e : rows[o0 + i]) {
                if (e.first >= o0 && e.first < o1) M(i, e.first - o0) = -e.second;
                else if (e.first >= o1) colPos[e.first] = 0;
            }
        check(M, phases[k], k);
        fac.emplace_back(M, phases[k], k, false);
        std::vector<std::size_t> cols;
        for (auto& c : colPos) {
            c.second = cols.size();
            cols.push_back(c.first);
        }
        Matrix<T> QkB(nk, cols.size(), zero<T>());
        for (std::size_t i = 0; i < nk; ++i)
            for (const auto& e : rows[o0 + i])
                if (e.first >= o1) QkB(i, colPos[e.first]) = e.second;
        for (std::size_t r = o1; r < n; ++r) {
            std::vector<T> v;
            for (const auto& e : rows[r]) {
                if (e.first >= o0 && e.first < o1 && !(e.second == zero<T>())) {
                    if (v.empty()) v.assign(nk, zero<T>());
                    v[e.first - o0] = e.second;
                }
            }
            if (!v.empty()) {
                bRows[k].push_back(r);
                bVals[k].push_back(v);
            }
        }
        if (!cols.empty() && !bRows[k].empty()) {
            const Matrix<T> X = fac.back().colsolve(QkB);
            for (std::size_t i = 0; i < bRows[k].size(); ++i) {
                const std::vector<T> u = vecmul(bVals[k][i], X);
                auto& row = rows[bRows[k][i]];
                for (std::size_t j = 0; j < cols.size(); ++j) {
                    if (u[j] == zero<T>()) continue;
                    auto it = row.find(cols[j]);
                    if (it == row.end()) row[cols[j]] = u[j];
                    else it->second += u[j];
                }
            }
        }
    }
 
    // intermediate step: level reduction on the last LDQBD
    const std::size_t oS = setOff[S - 1], nS = nset[S - 1];
    Matrix<T> MS(nS, nS, zero<T>());
    for (std::size_t i = 0; i < nS; ++i)
        for (const auto& e : rows[oS + i])
            if (e.first >= oS) MS(i, e.first - oS) = -e.second;
    check(MS, phases[S - 1], S - 1);
    const Factor<T> f(MS, phases[S - 1], S - 1, true);
    const std::size_t L = f.levels();
    Matrix<T> A = f.S[L - 1];
    for (std::size_t i = 0; i < A.rows(); ++i) A(i, 0) = num_traits<T>::from_int(1);
    // x A = e1  <=>  A' x' = e1'
    Matrix<T> At(A.cols(), A.rows(), zero<T>());
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j) At(j, i) = A(i, j);
    std::vector<T> e1(A.rows(), zero<T>());
    e1[0] = num_traits<T>::from_int(1);
    std::vector<T> xL;
    try {
        xL = solve(At, e1);
    } catch (const NumericError&) {
        throw NumericError(ErrorKind{"Singular"}, "ldqbd_censor: the last level of set " + std::to_string(S) +
                                                      " has no unique stationary vector; the chain may be reducible.");
    }
    std::vector<T> pB(n, zero<T>());
    for (std::size_t i = 0; i < xL.size(); ++i) pB[oS + f.off[L - 1] + i] = xL[i];
    for (std::size_t l = L - 1; l-- > 0;) {
        const std::vector<T> xn(pB.begin() + oS + f.off[l + 1], pB.begin() + oS + f.off[l + 2]);
        const std::vector<T> xg = vecmul(xn, f.G[l + 1]);
        for (std::size_t i = 0; i < f.P[l]; ++i) pB[oS + f.off[l] + i] = -xg[i];
    }
 
    // expansion back through sets S-1..1
    for (std::size_t k = S - 1; k-- > 0;) {
        std::vector<T> y(nset[k], zero<T>());
        for (std::size_t i = 0; i < bRows[k].size(); ++i) {
            const T p = pB[bRows[k][i]];
            if (p == zero<T>()) continue;
            for (std::size_t t = 0; t < nset[k]; ++t) y[t] += p * bVals[k][i][t];
        }
        const std::vector<T> x = fac[k].rowsolve(y);
        for (std::size_t t = 0; t < nset[k]; ++t) pB[setOff[k] + t] = x[t];
    }
    T tot = zero<T>();
    for (const T& v : pB) tot += v;
    LdqbdCensorResult<T> res;
    res.pi.resize(n);
    for (std::size_t i = 0; i < n; ++i) res.pi[i] = pB[i] / tot;
    std::size_t pos = 0;
    res.pi_cells.resize(S);
    for (std::size_t s = 0; s < S; ++s)
        for (std::size_t p : phases[s]) {
            res.pi_cells[s].emplace_back(res.pi.begin() + pos, res.pi.begin() + pos + p);
            pos += p;
        }
    return res;
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_API_MAM_LDQBD_CENSOR_H