ctmc_passage.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_PASSAGE_H
#define LINE_API_MC_CTMC_PASSAGE_H
 
/**
 * @file
 * @ingroup api_mc
 * First passage times into a target STATE SET, for Markov and semi-Markov
 * chains.
 *
 * Port of matlab/src/api/mc/ctmc_passage_*.m and smp_passage_*.m, which are the
 * reference. Source: P. G. Harrison and W. J. Knottenbelt, "Passage Time
 * Distributions in Large Markov Chains", 2002 -- Eqs. 1-3 for the Markov case,
 * Eqs. 4-8 for the semi-Markov one.
 *
 * THE IDENTITY THE WHOLE FAMILY RESTS ON. The first passage time from an
 * initial law pi0 into a target set B is PHASE-TYPE. With A the complement,
 *
 *     S     = Q(A,A)        sub-generator: the passage has not completed
 *     s0    = -S*1          exit vector, equal to Q(A,B)*1
 *     alpha = pi0(A)        UNNORMALIZED, see below
 *     atom  = sum pi0(B)
 *
 * so L(s) = alpha (sI-S)^{-1} s0 + atom and F(t) = 1 - alpha exp(St) 1. The
 * paper writes the same system as n scalar equations with L_i = 1 on B.
 *
 * ALPHA IS DELIBERATELY NOT NORMALIZED. Its mass is 1 - atom; the missing mass
 * is the ATOM AT ZERO carried by initial states already inside the target. A
 * caller that normalizes alpha and forgets the atom reports F(0) = 0 for a
 * passage that has already completed with probability atom.
 *
 * THIS IS NOT THE SPLIT `solver_ctmc_cdf.h` USES. That one is by EVENT (the
 * tagged job arriving at or departing from a station, through the filtration),
 * this one is by STATE SET. The two are complementary and must not be merged.
 *
 * ARITHMETIC: needs `expm` and transcendentals, so exact arithmetic is refused
 * by name rather than silently producing a wrong type.
 */
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/lti/laplace_invert.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/api/mc/dtmc_solve.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace mc {
 
/** The phase-type form of a first passage time. */
template <class T>
struct PassagePh {
    std::vector<T> alpha;        ///< pi0 restricted to the non-target block, UNNORMALIZED
    Matrix<T> S;                 ///< sub-generator Q(A,A)
    std::vector<T> s0;           ///< exit vector -S*1
    std::vector<std::size_t> keep;  ///< row of S -> state index of Q
    T atom;                      ///< mass of pi0 already inside the target: F(0)
};
 
/** Moments of a first passage time, `[Mall, m]` of the reference. */
template <class T>
struct PassageMoments {
    Matrix<T> mall;      ///< (nstates x nmax), zero on the target, inf where unreachable
    std::vector<T> m;    ///< the pi0-weighted moment vector
};
 
namespace passage_detail {
 
inline std::vector<std::size_t> unique_target(const std::vector<std::size_t>& target,
                                              std::size_t n, const char* fn) {
    std::vector<std::size_t> t = target;
    std::sort(t.begin(), t.end());
    t.erase(std::unique(t.begin(), t.end()), t.end());
    if (t.empty())
        throw InputError(std::string(fn) +
                         ": the target state set is empty: a first passage time into no state "
                         "is undefined");
    if (t.back() >= n)
        throw InputError(std::string(fn) + ": a target state index is outside the state space");
    return t;
}
 
/**
 * Gaussian elimination with partial pivoting for a COMPLEX system.
 *
 * `line::solve` cannot serve here: its pivot search compares magnitudes with
 * `operator>`, which std::complex does not have, so instantiating it on
 * complex is a compile error rather than a silent wrong answer. The transform
 * routes need complex coefficients, so they carry their own solve.
 */
inline std::vector<std::complex<double>> solve_cplx(Matrix<std::complex<double>> A,
                                                    std::vector<std::complex<double>> b,
                                                    const char* fn) {
    using C = std::complex<double>;
    const std::size_t n = A.rows();
    if (A.cols() != n || b.size() != n)
        throw InputError(std::string(fn) + ": complex system is not square");
    for (std::size_t k = 0; k < n; ++k) {
        std::size_t piv = k;
        double best = std::abs(A(k, k));
        for (std::size_t i = k + 1; i < n; ++i) {
            const double m = std::abs(A(i, k));
            if (m > best) {
                best = m;
                piv = i;
            }
        }
        if (!(best > 0.0))
            throw NumericError(ErrorKind{"Singular"}, std::string(fn) + ": singular transform matrix");
        if (piv != k) {
            for (std::size_t j = 0; j < n; ++j) std::swap(A(k, j), A(piv, j));
            std::swap(b[k], b[piv]);
        }
        for (std::size_t i = k + 1; i < n; ++i) {
            const C f = A(i, k) / A(k, k);
            if (f == C(0.0, 0.0)) continue;
            for (std::size_t j = k; j < n; ++j) A(i, j) -= f * A(k, j);
            b[i] -= f * b[k];
        }
    }
    std::vector<C> x(n, C(0.0, 0.0));
    for (std::size_t i = n; i-- > 0;) {
        C s = b[i];
        for (std::size_t j = i + 1; j < n; ++j) s -= A(i, j) * x[j];
        x[i] = s / A(i, i);
    }
    return x;
}
 
/**
 * The reference's `m = pi0 * Mall` of smp_passage_moments: empty when pi0 is
 * not given, refused when it is not a law over the state space.
 */
template <class T>
void smp_weighted_moments(PassageMoments<T>& out, const std::vector<T>& pi0, std::size_t n,
                          std::size_t nmax, const char* fn) {
    out.m.clear();
    if (pi0.empty()) return;
    if (pi0.size() != n)
        throw InputError(std::string(fn) +
                         ": PI0 must be a distribution over the state space, one entry per state");
    out.m.assign(nmax, T(0));
    for (std::size_t q = 0; q < nmax; ++q) {
        T acc = T(0);
        for (std::size_t i = 0; i < n; ++i) acc += pi0[i] * out.mall(i, q);
        out.m[q] = acc;
    }
}
 
/** Backward reachability closure over the transition graph. */
template <class T>
std::vector<bool> reaches_target(const Matrix<T>& Q, const std::vector<std::size_t>& keep,
                                 const std::vector<std::size_t>& target) {
    const std::size_t n = Q.rows();
    std::vector<bool> seen(n, false);
    std::vector<std::size_t> frontier;
    for (std::size_t k : target) {
        seen[k] = true;
        frontier.push_back(k);
    }
    while (!frontier.empty()) {
        std::vector<std::size_t> next;
        for (std::size_t j : frontier)
            for (std::size_t i = 0; i < n; ++i)
                if (i != j && !seen[i] && Q(i, j) != T(0)) {
                    seen[i] = true;
                    next.push_back(i);
                }
        frontier.swap(next);
    }
    std::vector<bool> out(keep.size(), false);
    for (std::size_t a = 0; a < keep.size(); ++a) out[a] = seen[keep[a]];
    return out;
}
 
}  // namespace passage_detail
 
/**
 * Phase-type representation of the first passage time from `pi0` into `target`.
 *
 * @param pi0 empty selects the conditional stationary law on the complement
 * @param target 0-based state indices (1-based in the MATLAB reference)
 */
template <class T>
PassagePh<T> ctmc_passage_ph(const Matrix<T>& Q, const std::vector<T>& pi0,
                             const std::vector<std::size_t>& target) {
    const std::size_t n = Q.rows();
    if (Q.cols() != n) throw InputError("ctmc_passage_ph: the generator must be square");
    const std::vector<std::size_t> tgt =
        passage_detail::unique_target(target, n, "ctmc_passage_ph");
 
    double scale = 1.0;
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            scale = std::max(scale, std::abs(num_traits<T>::to_double(Q(i, j))));
    for (std::size_t i = 0; i < n; ++i) {
        T rs = T(0);
        for (std::size_t j = 0; j < n; ++j) rs += Q(i, j);
        if (std::abs(num_traits<T>::to_double(rs)) > 1e-8 * scale)
            throw InputError(
                "ctmc_passage_ph: Q is not an infinitesimal generator: its rows do not sum to "
                "zero. Pass it through ctmc_makeinfgen first");
    }
 
    std::vector<bool> is_target(n, false);
    for (std::size_t k : tgt) is_target[k] = true;
    PassagePh<T> out;
    for (std::size_t i = 0; i < n; ++i)
        if (!is_target[i]) out.keep.push_back(i);
 
    const std::size_t nA = out.keep.size();
    out.S = Matrix<T>(nA, nA);
    for (std::size_t a = 0; a < nA; ++a)
        for (std::size_t c = 0; c < nA; ++c) out.S(a, c) = Q(out.keep[a], out.keep[c]);
    out.s0.assign(nA, T(0));
    for (std::size_t a = 0; a < nA; ++a) {
        T r = T(0);
        for (std::size_t c = 0; c < nA; ++c) r += out.S(a, c);
        out.s0[a] = -r;
    }
 
    out.alpha.assign(nA, T(0));
    out.atom = T(0);
    if (pi0.empty()) {
        // An empty initial law selects the conditional stationary one on the
        // complement of the target set, the contract of the MATLAB reference
        // and of the python and JAR api twins.
        const std::vector<T> p = ctmc_solve(Q);
        T mass = T(0);
        for (std::size_t a = 0; a < nA; ++a) mass += p[out.keep[a]];
        if (!(num_traits<T>::to_double(mass) > 0.0))
            throw InputError(
                "ctmc_passage_ph: the stationary law puts no mass outside the target set, so "
                "there is no passage to time");
        for (std::size_t a = 0; a < nA; ++a) out.alpha[a] = p[out.keep[a]] / mass;
        return out;
    }
    if (pi0.size() != n)
        throw InputError(
            "ctmc_passage_ph: pi0 must be a distribution over the state space, one entry per "
            "state");
    for (std::size_t a = 0; a < nA; ++a) out.alpha[a] = pi0[out.keep[a]];
    for (std::size_t k : tgt) out.atom += pi0[k];
    return out;
}
 
/**
 * L(s) = alpha (sI-S)^{-1} s0 + atom at the (complex) points `s`. Eqs. 1-2: one
 * linear system per value of s.
 *
 * ONE SOLVE PER s, NOT PER (s,t) PAIR. The saving over a dense matrix
 * exponential is that the solves are sparse, so this route reaches chains a
 * dense expm cannot hold. It is NOT a saving in the number of time points:
 * every Abate-Whitt inverter places its nodes at s = beta/t, so a grid of T
 * points costs T*|beta| solves.
 */
template <class T>
std::vector<std::complex<double>> ctmc_passage_lst(const Matrix<T>& Q, const std::vector<T>& pi0,
                                                   const std::vector<std::size_t>& target,
                                                   const std::vector<std::complex<double>>& s) {
    const PassagePh<T> ph = ctmc_passage_ph(Q, pi0, target);
    const std::size_t nA = ph.S.rows();
    using C = std::complex<double>;
    std::vector<C> out(s.size(), C(0.0, 0.0));
    for (std::size_t is = 0; is < s.size(); ++is) {
        Matrix<C> A(nA, nA);
        for (std::size_t i = 0; i < nA; ++i)
            for (std::size_t j = 0; j < nA; ++j)
                A(i, j) = (i == j ? s[is] : C(0.0, 0.0)) -
                          C(num_traits<T>::to_double(ph.S(i, j)), 0.0);
        std::vector<C> b(nA);
        for (std::size_t i = 0; i < nA; ++i) b[i] = C(num_traits<T>::to_double(ph.s0[i]), 0.0);
        const std::vector<C> x = passage_detail::solve_cplx(A, b, "ctmc_passage_lst");
        C acc(num_traits<T>::to_double(ph.atom), 0.0);
        for (std::size_t i = 0; i < nA; ++i)
            acc += C(num_traits<T>::to_double(ph.alpha[i]), 0.0) * x[i];
        out[is] = acc;
    }
    return out;
}
 
/**
 * Moments of order 1..nmax of the first passage time into `target`.
 *
 * This is Eq. 3, -q_ii M_i(n) = sum_{k not in B} q_ik M_k(n) + n M_i(n-1),
 * i.e. (-S) M(n) = n M(n-1) with M(0) = 1: nmax linear solves and no transform
 * inversion at all. The equivalent closed form n! alpha (-S)^{-n} 1 is NOT how
 * it is evaluated -- forming the inverse of the sub-generator destroys the
 * sparsity the recursion preserves.
 */
template <class T>
PassageMoments<T> ctmc_passage_moments(const Matrix<T>& Q, const std::vector<T>& pi0,
                                       const std::vector<std::size_t>& target,
                                       std::size_t nmax = 1) {
    static_assert(num_traits<T>::has_transcendental,
                  "ctmc_passage_moments marks unreachable states with an infinity and therefore "
                  "requires an arithmetic that has one");
    if (nmax == 0) throw InputError("ctmc_passage_moments: nmax must be positive");
    const PassagePh<T> ph = ctmc_passage_ph(Q, pi0, target);
    const std::size_t n = Q.rows();
    const std::vector<std::size_t> tgt =
        passage_detail::unique_target(target, n, "ctmc_passage_moments");
    const std::size_t nA = ph.keep.size();
 
    PassageMoments<T> out;
    out.mall = Matrix<T>(n, nmax);
    out.m.assign(nmax, T(0));
    if (nA == 0) return out;
 
    // A state that cannot reach the target has an infinite passage time; the
    // sub-generator is singular on that block, and a solve that ignored this
    // would return a finite number instead of saying so.
    const std::vector<bool> reach = passage_detail::reaches_target(Q, ph.keep, tgt);
    bool all_reach = true;
    for (bool r : reach) all_reach = all_reach && r;
 
    Matrix<T> A(nA, nA);
    for (std::size_t i = 0; i < nA; ++i)
        for (std::size_t j = 0; j < nA; ++j) A(i, j) = -ph.S(i, j);
 
    const T inf = num_traits<T>::from_double(std::numeric_limits<double>::infinity());
    std::vector<T> x(nA, T(1));
    for (std::size_t k = 1; k <= nmax; ++k) {
        std::vector<T> rhs(nA);
        for (std::size_t i = 0; i < nA; ++i) rhs[i] = num_traits<T>::from_double(double(k)) * x[i];
        if (all_reach) {
            x = solve(A, rhs);
        } else {
            // Restrict to the reachable block: the unreachable rows are exactly
            // the singular ones, and they are reported as infinite rather than
            // regularized away.
            std::vector<std::size_t> idx;
            for (std::size_t i = 0; i < nA; ++i)
                if (reach[i]) idx.push_back(i);
            Matrix<T> Ar(idx.size(), idx.size());
            std::vector<T> br(idx.size());
            for (std::size_t a = 0; a < idx.size(); ++a) {
                for (std::size_t c = 0; c < idx.size(); ++c) Ar(a, c) = A(idx[a], idx[c]);
                br[a] = rhs[idx[a]];
            }
            const std::vector<T> xr = solve(Ar, br);
            x.assign(nA, inf);
            for (std::size_t a = 0; a < idx.size(); ++a) x[idx[a]] = xr[a];
        }
        for (std::size_t i = 0; i < nA; ++i)
            out.mall(ph.keep[i], k - 1) = reach[i] ? x[i] : inf;
        if (!all_reach)
            for (std::size_t i = 0; i < nA; ++i)
                if (!reach[i]) x[i] = T(1);  // keeps the recursion finite on the reachable block
    }
 
    bool unreachable_start = false;
    for (std::size_t i = 0; i < nA; ++i)
        if (!reach[i] && num_traits<T>::to_double(ph.alpha[i]) > 0.0) unreachable_start = true;
    for (std::size_t k = 0; k < nmax; ++k) {
        if (unreachable_start) {
            out.m[k] = inf;
            continue;
        }
        T acc = T(0);
        for (std::size_t i = 0; i < nA; ++i)
            if (reach[i]) acc += ph.alpha[i] * out.mall(ph.keep[i], k);
        out.m[k] = acc;
    }
    return out;
}
 
/**
 * Mean time to reach any state in `target` from each state of a CTMC.
 *
 * Continuous-time twin of `dtmc_hitting_time` and the first-moment special case
 * of `ctmc_passage_moments`: (-S) h = 1 on the non-target block, where
 * `dtmc_hitting_time` solves (I - P_NT) h = 1. Unreachable states give infinity.
 */
template <class T>
std::vector<T> ctmc_hitting_time(const Matrix<T>& Q, const std::vector<std::size_t>& target) {
    const std::size_t n = Q.rows();
    // mall does not depend on the initial law, so a uniform one is passed
    // rather than requiring the caller to invent one.
    std::vector<T> pi0(n, num_traits<T>::from_double(1.0 / double(n)));
    const PassageMoments<T> pm = ctmc_passage_moments(Q, pi0, target, 1);
    std::vector<T> h(n);
    for (std::size_t i = 0; i < n; ++i) h[i] = pm.mall(i, 0);
    return h;
}
 
/** A passage-time law on a grid. */
template <class T>
struct PassageCurve {
    std::vector<double> t;
    std::vector<double> F;
    std::vector<double> f;
    double atom = 0.0;
};
 
/**
 * CDF and density of the first passage time on the grid `tset`:
 * F(t) = 1 - alpha exp(St) 1 and f(t) = alpha exp(St) s0.
 *
 * `method` is "expm" (default) or "lt". The transform route exists for chains
 * whose non-target block is too large for a dense exp(St), not because it needs
 * fewer time points; on a small chain "expm" is both faster and more accurate,
 * which is why it is the default.
 */
template <class T>
PassageCurve<T> ctmc_passage_time(const Matrix<T>& Q, const std::vector<T>& pi0,
                                  const std::vector<std::size_t>& target,
                                  const std::vector<double>& tset,
                                  const std::string& method = "expm",
                                  const std::string& lti_method = "euler") {
    static_assert(num_traits<T>::has_transcendental,
                  "ctmc_passage_time takes a matrix exponential and therefore requires an "
                  "arithmetic with transcendental functions");
    const PassagePh<T> ph = ctmc_passage_ph(Q, pi0, target);
    const std::size_t nA = ph.S.rows();
    PassageCurve<T> out;
    out.t = tset;
    out.F.assign(tset.size(), 0.0);
    out.f.assign(tset.size(), 0.0);
    out.atom = num_traits<T>::to_double(ph.atom);
 
    if (method == "expm") {
        bool uniform = tset.size() > 2;
        const double dt = tset.size() > 1 ? tset[1] - tset[0] : 0.0;
        if (uniform && !(dt > 0.0)) uniform = false;
        for (std::size_t i = 1; uniform && i + 1 < tset.size(); ++i)
            if (std::abs((tset[i + 1] - tset[i]) - dt) > 1e-12 * std::max(1.0, std::abs(dt)))
                uniform = false;
 
        std::vector<double> v(nA, 0.0);
        Matrix<double> Sd(nA, nA);
        for (std::size_t i = 0; i < nA; ++i)
            for (std::size_t j = 0; j < nA; ++j) Sd(i, j) = num_traits<T>::to_double(ph.S(i, j));
        std::vector<double> s0d(nA), ad(nA);
        for (std::size_t i = 0; i < nA; ++i) {
            s0d[i] = num_traits<T>::to_double(ph.s0[i]);
            ad[i] = num_traits<T>::to_double(ph.alpha[i]);
        }
 
        auto step = [&](const Matrix<double>& E, std::vector<double>& w) {
            std::vector<double> z(nA, 0.0);
            for (std::size_t j = 0; j < nA; ++j) {
                double acc = 0.0;
                for (std::size_t i = 0; i < nA; ++i) acc += w[i] * E(i, j);
                z[j] = acc;
            }
            w.swap(z);
        };
 
        if (uniform) {
            // One exponential, then propagate: recomputing expm(S*t) at every
            // grid point is the same answer at a cost linear in the grid.
            Matrix<double> Sdt = Sd;
            for (std::size_t i = 0; i < nA; ++i)
                for (std::size_t j = 0; j < nA; ++j) Sdt(i, j) = Sd(i, j) * dt;
            const Matrix<double> E = expm(Sdt);
            Matrix<double> S0 = Sd;
            for (std::size_t i = 0; i < nA; ++i)
                for (std::size_t j = 0; j < nA; ++j) S0(i, j) = Sd(i, j) * tset[0];
            const Matrix<double> E0 = expm(S0);
            v.assign(nA, 0.0);
            for (std::size_t j = 0; j < nA; ++j) {
                double acc = 0.0;
                for (std::size_t i = 0; i < nA; ++i) acc += ad[i] * E0(i, j);
                v[j] = acc;
            }
            for (std::size_t i = 0; i < tset.size(); ++i) {
                if (i > 0) step(E, v);
                double sF = 0.0, sf = 0.0;
                for (std::size_t j = 0; j < nA; ++j) {
                    sF += v[j];
                    sf += v[j] * s0d[j];
                }
                out.F[i] = 1.0 - sF;
                out.f[i] = sf;
            }
        } else {
            for (std::size_t i = 0; i < tset.size(); ++i) {
                if (tset[i] < 0.0) continue;
                Matrix<double> St = Sd;
                for (std::size_t a = 0; a < nA; ++a)
                    for (std::size_t b = 0; b < nA; ++b) St(a, b) = Sd(a, b) * tset[i];
                const Matrix<double> E = expm(St);
                double sF = 0.0, sf = 0.0;
                for (std::size_t j = 0; j < nA; ++j) {
                    double acc = 0.0;
                    for (std::size_t a = 0; a < nA; ++a) acc += ad[a] * E(a, j);
                    sF += acc;
                    sf += acc * s0d[j];
                }
                out.F[i] = 1.0 - sF;
                out.f[i] = sf;
            }
        }
    } else if (method == "lt") {
        const double atom = out.atom;
        const lti::LaplaceFn L = [&](std::complex<double> s) {
            std::vector<std::complex<double>> sv(1, s);
            return ctmc_passage_lst(Q, pi0, target, sv)[0];
        };
        const lti::LaplaceMethod lm = lti::laplace_method(lti_method);
        out.F = lti::laplace_invert_cdf(L, tset, lm);
        const lti::LaplaceFn Ld = [&](std::complex<double> s) { return L(s) - atom; };
        out.f = lti::laplace_invert_pdf(Ld, tset, lm);
    } else {
        throw InputError("ctmc_passage_time: unknown method '" + method +
                         "', expected expm or lt");
    }
 
    for (std::size_t i = 0; i < out.F.size(); ++i) {
        out.F[i] = std::min(1.0, std::max(0.0, out.F[i]));
        out.f[i] = std::max(0.0, out.f[i]);
    }
    return out;
}
 
// ---------------------------------------------------------------------------
// First passage under a taboo set
// ---------------------------------------------------------------------------
 
/**
 * First passage into a target set B under a TABOO set C, `[F, f, out]` of
 * `ctmc_passage_taboo.m`: the passage completes by entering B and is lost for
 * good by entering C first.
 *
 * The chain is cut into four parts: B and C, made absorbing; N, the states that
 * can reach neither, also absorbing since they can only hold the passage up
 * forever; and the transient states K, on which S = Q(K,K) is nonsingular by
 * construction. A state in both B and C counts as B.
 */
template <class T>
struct PassageTaboo {
    std::vector<double> t;       ///< the grid the curve was asked on
    std::vector<double> F;       ///< DEFECTIVE CDF Pr{B entered by t, C not before}, clamped to [0, 1]
    std::vector<double> f;       ///< its density alpha exp(St) sB
    T atom;                      ///< mass of pi0 already in B, i.e. F(0)
    T reachprob;                 ///< Pr{B entered before C}, the limit of F
    T tabooprob;                 ///< Pr{C entered before B}
    std::vector<T> pB;           ///< (nstates) reachprob per starting state, one on B
    std::vector<T> pC;           ///< (nstates) tabooprob per starting state, one on C
    std::vector<T> alpha;        ///< pi0 restricted to K, UNNORMALIZED
    Matrix<T> S;                 ///< Q(K,K)
    std::vector<T> sB;           ///< Q(K,B)*1, the exit rates into the target
    std::vector<T> sC;           ///< Q(K,C)*1, the exit rates into the taboo set
    std::vector<std::size_t> keep;    ///< row of S -> state index of Q (the set K)
    std::vector<std::size_t> target;  ///< B, sorted and unique
    std::vector<std::size_t> taboo;   ///< C minus B, sorted and unique
    std::string method;
};
 
/** Conditional moments under a taboo set, `[mall, m, pB]` of `ctmc_passage_taboo_moments.m`. */
template <class T>
struct PassageTabooMoments {
    Matrix<T> mall;      ///< (nstates x nmax): zero on B, NaN where B cannot be entered first
    std::vector<T> m;    ///< the pi0-weighted conditional moments, NaN if B is never entered first
    std::vector<T> pB;   ///< (nstates) probability of entering B before C
};
 
namespace passage_detail {
 
/**
 * States with a path into `from` that enters no `stop` state before its last
 * step. A stop state is absorbing for the passage, so it is never a predecessor.
 */
template <class T>
std::vector<bool> backward_closure(const Matrix<T>& Q, const std::vector<bool>& from,
                                   const std::vector<bool>& stop) {
    const std::size_t n = Q.rows();
    std::vector<bool> seen = from;
    std::vector<std::size_t> frontier;
    for (std::size_t i = 0; i < n; ++i)
        if (from[i]) frontier.push_back(i);
    while (!frontier.empty()) {
        std::vector<std::size_t> next;
        for (std::size_t j : frontier)
            for (std::size_t i = 0; i < n; ++i)
                if (i != j && !seen[i] && !stop[i] && Q(i, j) != T(0)) {
                    seen[i] = true;
                    next.push_back(i);
                }
        frontier.swap(next);
    }
    return seen;
}
 
/** The partition, the sub-generator and the ultimate probabilities: everything but the curve. */
template <class T>
PassageTaboo<T> taboo_partition(const Matrix<T>& Q, const std::vector<T>& pi0,
                                const std::vector<std::size_t>& target,
                                const std::vector<std::size_t>& taboo, const char* fn) {
    const std::size_t n = Q.rows();
    if (Q.cols() != n) throw InputError(std::string(fn) + ": the generator must be square");
    PassageTaboo<T> out;
    out.target = unique_target(target, n, fn);
    std::vector<std::size_t> tb = taboo;
    std::sort(tb.begin(), tb.end());
    tb.erase(std::unique(tb.begin(), tb.end()), tb.end());
    if (!tb.empty() && tb.back() >= n)
        throw InputError(std::string(fn) + ": a target or taboo state index is outside the state space");
 
    double scale = 1.0;
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j)
            scale = std::max(scale, std::abs(num_traits<T>::to_double(Q(i, j))));
    for (std::size_t i = 0; i < n; ++i) {
        T rs = T(0);
        for (std::size_t j = 0; j < n; ++j) rs += Q(i, j);
        if (std::abs(num_traits<T>::to_double(rs)) > 1e-8 * scale)
            throw InputError(std::string(fn) +
                             ": Q is not an infinitesimal generator: its rows do not sum to zero. "
                             "Pass it through ctmc_makeinfgen first");
    }
 
    std::vector<bool> isB(n, false), isC(n, false);
    for (std::size_t k : out.target) isB[k] = true;
    for (std::size_t k : tb)
        if (!isB[k]) {
            isC[k] = true;
            out.taboo.push_back(k);
        }
 
    std::vector<T> p0;
    if (pi0.empty()) {
        // The conditional stationary law on the complement of B, as in ctmc_passage_ph.
        p0 = ctmc_solve(Q);
        T mass = T(0);
        for (std::size_t i = 0; i < n; ++i) {
            if (isB[i]) p0[i] = T(0);
            mass += p0[i];
        }
        if (!(num_traits<T>::to_double(mass) > 0.0))
            throw InputError(std::string(fn) +
                             ": the stationary law puts no mass outside the target set, so there "
                             "is no passage to time");
        for (std::size_t i = 0; i < n; ++i) p0[i] = T(p0[i] / mass);
    } else {
        if (pi0.size() != n)
            throw InputError(std::string(fn) +
                             ": PI0 must be a distribution over the state space, one entry per "
                             "state");
        p0 = pi0;
    }
 
    std::vector<bool> isBC(n, false);
    for (std::size_t i = 0; i < n; ++i) isBC[i] = isB[i] || isC[i];
    const std::vector<bool> toB = backward_closure(Q, isB, isBC);
    const std::vector<bool> toBC = backward_closure(Q, isBC, isBC);
 
    std::vector<std::size_t> absorb, Nset;
    for (std::size_t i = 0; i < n; ++i)
        if (isB[i]) absorb.push_back(i);
    for (std::size_t i = 0; i < n; ++i)
        if (isC[i]) absorb.push_back(i);
    for (std::size_t i = 0; i < n; ++i) {
        if (!toBC[i]) Nset.push_back(i);
        else if (!isBC[i]) out.keep.push_back(i);
    }
    absorb.insert(absorb.end(), Nset.begin(), Nset.end());
    const std::size_t nB = out.target.size(), nC = out.taboo.size(), nK = out.keep.size();
 
    out.pB.assign(n, T(0));
    out.pC.assign(n, T(0));
    for (std::size_t i = 0; i < n; ++i) {
        if (isB[i]) out.pB[i] = T(1);
        if (isC[i]) out.pC[i] = T(1);
    }
    out.S = Matrix<T>(nK, nK);
    if (nK > 0) {
        // The ultimate probabilities are the exit-probability matrix of the
        // stochastic complement onto B, C and N, with every non-K row zeroed so
        // that the three sets are absorbing.
        Matrix<T> Qabs = Q;
        std::vector<bool> isK(n, false);
        for (std::size_t k : out.keep) isK[k] = true;
        for (std::size_t i = 0; i < n; ++i)
            if (!isK[i])
                for (std::size_t j = 0; j < n; ++j) Qabs(i, j) = T(0);
        const StochCompResult<T> sc = ctmc_stochcomp(Qabs, absorb);
        out.S = sc.Q22;
        for (std::size_t a = 0; a < nK; ++a) {
            T b = T(0), c = T(0);
            for (std::size_t j = 0; j < nB; ++j) b += sc.X(a, j);
            for (std::size_t j = 0; j < nC; ++j) c += sc.X(a, nB + j);
            if (toB[out.keep[a]]) out.pB[out.keep[a]] = b;
            out.pC[out.keep[a]] = c;
        }
    }
 
    out.alpha.assign(nK, T(0));
    out.sB.assign(nK, T(0));
    out.sC.assign(nK, T(0));
    for (std::size_t a = 0; a < nK; ++a) {
        const std::size_t i = out.keep[a];
        out.alpha[a] = p0[i];
        for (std::size_t j : out.target) out.sB[a] += Q(i, j);
        for (std::size_t j : out.taboo) out.sC[a] += Q(i, j);
    }
    out.atom = T(0);
    for (std::size_t j : out.target) out.atom += p0[j];
    out.reachprob = T(0);
    out.tabooprob = T(0);
    for (std::size_t i = 0; i < n; ++i) {
        out.reachprob += T(p0[i] * out.pB[i]);
        out.tabooprob += T(p0[i] * out.pC[i]);
    }
    return out;
}
 
}  // namespace passage_detail
 
/**
 * First passage into `target` under the taboo set `taboo`, the reference's
 * `ctmc_passage_taboo`. F is the DEFECTIVE distribution
 *
 *     F(t) = Pr{ B is entered by time t, and C is not entered before B }
 *          = atom + alpha (-S)^{-1} (I - exp(St)) sB
 *
 * on the grid `tset`, and f its density alpha exp(St) sB; F tends to
 * `reachprob` < 1 rather than to one. In temporal-logic terms this is the
 * time-bounded until Pr{ not C U<=t B }.
 *
 * The ultimate probabilities come from the exit-probability matrix X of
 * `ctmc_stochcomp`, so they are one linear solve and need no time grid. An empty
 * `tset` computes only those; an infinite entry returns `reachprob` and a
 * negative one zero. With an empty `taboo`, F equals the CDF of
 * `ctmc_passage_time`.
 *
 * @param pi0 empty selects the conditional stationary law on the complement of B
 * @param target 0-based state indices of B (1-based in the MATLAB reference)
 * @param taboo 0-based state indices of C; a state also in B counts as B
 * @param method "expm" (default) or "lt", as in `ctmc_passage_time`
 *
 * Reference: C. Baier, B. Haverkort, H. Hermanns, J.-P. Katoen, "Model-checking
 * algorithms for continuous-time Markov chains", IEEE TSE 29(6), 2003.
 */
template <class T>
PassageTaboo<T> ctmc_passage_taboo(const Matrix<T>& Q, const std::vector<T>& pi0,
                                   const std::vector<std::size_t>& target,
                                   const std::vector<std::size_t>& taboo,
                                   const std::vector<double>& tset = std::vector<double>(),
                                   const std::string& method = "expm",
                                   const std::string& lti_method = "euler") {
    static_assert(num_traits<T>::has_transcendental,
                  "ctmc_passage_taboo takes a matrix exponential and therefore requires an "
                  "arithmetic with transcendental functions");
    PassageTaboo<T> out =
        passage_detail::taboo_partition(Q, pi0, target, taboo, "ctmc_passage_taboo");
    out.method = method;
    out.t = tset;
    out.F.assign(tset.size(), 0.0);
    out.f.assign(tset.size(), 0.0);
    if (tset.empty()) return out;
 
    const double inf = std::numeric_limits<double>::infinity();
    const double atom = num_traits<T>::to_double(out.atom);
    std::vector<std::size_t> pos;  // finite, nonnegative entries
    for (std::size_t i = 0; i < tset.size(); ++i) {
        if (tset[i] == inf) out.F[i] = num_traits<T>::to_double(out.reachprob);
        else if (std::isfinite(tset[i]) && tset[i] >= 0.0) pos.push_back(i);
    }
    std::vector<double> tf(pos.size());
    for (std::size_t i = 0; i < pos.size(); ++i) tf[i] = tset[pos[i]];
    std::vector<double> Ff(tf.size(), atom), ff(tf.size(), 0.0);
 
    const std::size_t nK = out.keep.size();
    if (nK > 0 && !tf.empty()) {
        std::vector<double> ad(nK), sBd(nK);
        for (std::size_t a = 0; a < nK; ++a) {
            ad[a] = num_traits<T>::to_double(out.alpha[a]);
            sBd[a] = num_traits<T>::to_double(out.sB[a]);
        }
        if (method == "expm") {
            // One absorbing column for B appended to S: its mass at time t is the
            // probability of having entered B by t without leaving K elsewhere first.
            const std::size_t m = nK + 1;
            Matrix<double> G(m, m, 0.0);
            for (std::size_t a = 0; a < nK; ++a) {
                for (std::size_t c = 0; c < nK; ++c) G(a, c) = num_traits<T>::to_double(out.S(a, c));
                G(a, nK) = sBd[a];
            }
            auto scaled = [&](double s) {
                Matrix<double> Gs(m, m);
                for (std::size_t a = 0; a < m; ++a)
                    for (std::size_t c = 0; c < m; ++c) Gs(a, c) = G(a, c) * s;
                return expm(Gs);
            };
            auto times = [&](const std::vector<double>& w, const Matrix<double>& E) {
                std::vector<double> z(m, 0.0);
                for (std::size_t c = 0; c < m; ++c) {
                    double acc = 0.0;
                    for (std::size_t a = 0; a < m; ++a) acc += w[a] * E(a, c);
                    z[c] = acc;
                }
                return z;
            };
            auto record = [&](std::size_t i, const std::vector<double>& v) {
                Ff[i] = atom + v[nK];
                double acc = 0.0;
                for (std::size_t a = 0; a < nK; ++a) acc += v[a] * sBd[a];
                ff[i] = acc;
            };
            std::vector<double> v0(m, 0.0);
            for (std::size_t a = 0; a < nK; ++a) v0[a] = ad[a];
 
            bool uniform = tf.size() > 2;
            const double dt = tf.size() > 1 ? tf[1] - tf[0] : 0.0;
            if (uniform && !(dt > 0.0)) uniform = false;
            for (std::size_t i = 1; uniform && i + 1 < tf.size(); ++i)
                if (!(std::abs((tf[i + 1] - tf[i]) - dt) < 1e-12 * std::max(1.0, std::abs(dt))))
                    uniform = false;
            if (uniform) {
                const Matrix<double> E = scaled(dt);
                std::vector<double> v = times(v0, scaled(tf[0]));
                for (std::size_t i = 0; i < tf.size(); ++i) {
                    if (i > 0) v = times(v, E);
                    record(i, v);
                }
            } else {
                for (std::size_t i = 0; i < tf.size(); ++i) record(i, times(v0, scaled(tf[i])));
            }
        } else if (method == "lt") {
            using C = std::complex<double>;
            const Matrix<T>& S = out.S;
            const lti::LaplaceFn L = [&](C s) {
                Matrix<C> A(nK, nK);
                for (std::size_t a = 0; a < nK; ++a)
                    for (std::size_t c = 0; c < nK; ++c)
                        A(a, c) = (a == c ? s : C(0.0, 0.0)) -
                                  C(num_traits<T>::to_double(S(a, c)), 0.0);
                std::vector<C> b(nK);
                for (std::size_t a = 0; a < nK; ++a) b[a] = C(sBd[a], 0.0);
                const std::vector<C> x = passage_detail::solve_cplx(A, b, "ctmc_passage_taboo");
                C acc(atom, 0.0);
                for (std::size_t a = 0; a < nK; ++a) acc += C(ad[a], 0.0) * x[a];
                return acc;
            };
            const lti::LaplaceMethod lm = lti::laplace_method(lti_method);
            Ff = lti::laplace_invert_cdf(L, tf, lm);
            const lti::LaplaceFn Ld = [&](C s) { return L(s) - atom; };
            ff = lti::laplace_invert_pdf(Ld, tf, lm);
        } else {
            throw InputError("ctmc_passage_taboo: unknown passage-time method '" + method +
                             "', expected expm or lt");
        }
    }
    for (std::size_t i = 0; i < pos.size(); ++i) {
        out.F[pos[i]] = Ff[i];
        out.f[pos[i]] = ff[i];
    }
    for (std::size_t i = 0; i < out.F.size(); ++i) {
        out.F[i] = std::min(1.0, std::max(0.0, out.F[i]));
        out.f[i] = std::max(0.0, out.f[i]);
    }
    return out;
}
 
/**
 * Moments of order 1..nmax of the first passage time into `target`, CONDITIONAL
 * on entering it before `taboo`, the reference's `ctmc_passage_taboo_moments`:
 *
 *     E[T^k | B before C] = k! alpha (-S)^{-(k+1)} sB / Pr{B before C}
 *
 * evaluated as (-S) y_k = k y_(k-1) from y_0 = pB, one solve per order on one
 * factorization, with S, pB and the partition of `ctmc_passage_taboo`.
 *
 * `mall` is zero on B and NaN where B cannot be entered before C, since the
 * conditional moment is then undefined; `m` is NaN when B is never entered
 * first. With an empty `taboo` and B entered surely, the moments equal those of
 * `ctmc_passage_moments`.
 */
template <class T>
PassageTabooMoments<T> ctmc_passage_taboo_moments(const Matrix<T>& Q, const std::vector<T>& pi0,
                                                  const std::vector<std::size_t>& target,
                                                  const std::vector<std::size_t>& taboo,
                                                  std::size_t nmax = 1) {
    static_assert(num_traits<T>::has_transcendental,
                  "ctmc_passage_taboo_moments marks undefined conditional moments with a NaN "
                  "and therefore requires an arithmetic that has one");
    if (nmax == 0) throw InputError("ctmc_passage_taboo_moments: nmax must be positive");
    const PassageTaboo<T> pt =
        passage_detail::taboo_partition(Q, pi0, target, taboo, "ctmc_passage_taboo_moments");
    const std::size_t n = Q.rows(), nK = pt.keep.size();
    const T nan = num_traits<T>::from_double(std::numeric_limits<double>::quiet_NaN());
 
    PassageTabooMoments<T> out;
    out.pB = pt.pB;
    out.mall = Matrix<T>(n, nmax, nan);
    for (std::size_t i : pt.target)
        for (std::size_t k = 0; k < nmax; ++k) out.mall(i, k) = T(0);
 
    Matrix<T> Y(nK, nmax, T(0));
    if (nK > 0) {
        Matrix<T> LU(nK, nK);
        for (std::size_t a = 0; a < nK; ++a)
            for (std::size_t c = 0; c < nK; ++c) LU(a, c) = T(-pt.S(a, c));
        const std::vector<std::size_t> piv = lu_factor(LU);
        std::vector<T> y(nK);
        for (std::size_t a = 0; a < nK; ++a) y[a] = pt.pB[pt.keep[a]];
        for (std::size_t k = 1; k <= nmax; ++k) {
            const T kk = num_traits<T>::from_double(static_cast<double>(k));
            for (std::size_t a = 0; a < nK; ++a) y[a] = T(kk * y[a]);
            lu_solve(LU, piv, y);
            for (std::size_t a = 0; a < nK; ++a) Y(a, k - 1) = y[a];
        }
        for (std::size_t a = 0; a < nK; ++a) {
            const T p = pt.pB[pt.keep[a]];
            if (!(num_traits<T>::to_double(p) > 0.0)) continue;
            for (std::size_t k = 0; k < nmax; ++k) out.mall(pt.keep[a], k) = T(Y(a, k) / p);
        }
    }
 
    T PB = pt.atom;
    for (std::size_t a = 0; a < nK; ++a) PB += T(pt.alpha[a] * pt.pB[pt.keep[a]]);
    out.m.assign(nmax, nan);
    if (num_traits<T>::to_double(PB) > 0.0)
        for (std::size_t k = 0; k < nmax; ++k) {
            T acc = T(0);
            for (std::size_t a = 0; a < nK; ++a) acc += T(pt.alpha[a] * Y(a, k));
            out.m[k] = T(acc / PB);
        }
    return out;
}
 
// ---------------------------------------------------------------------------
// Semi-Markov chains (Sec. 3)
// ---------------------------------------------------------------------------
 
/**
 * Moments of the semi-Markov first passage time from the per-state holding
 * moments m_i(r), Eq. 7 with the u_i(r) recurrence of Eq. 8:
 *
 *     u_i(r) = -sum_{j=1..r} C(r,j) m_i(j) u_i(r-j),   u_i(0) = 1,
 *
 * which are the derivatives at the origin of 1/h*_i(s). Cheaper than Eq. 6
 * because it needs no per-pair moments.
 *
 * @param hmom (nstates x nmax): hmom(i,r-1) is the r-th moment of the sojourn
 *             in state i
 */
template <class T>
PassageMoments<T> smp_passage_moments(const Matrix<T>& P, const Matrix<T>& hmom,
                                      const std::vector<T>& pi0,
                                      const std::vector<std::size_t>& target,
                                      std::size_t nmax = 1) {
    const std::size_t n = P.rows();
    if (P.cols() != n)
        throw InputError("smp_passage_moments: the embedded transition matrix must be square");
    if (nmax == 0) throw InputError("smp_passage_moments: nmax must be positive");
    for (std::size_t i = 0; i < n; ++i) {
        T rs = T(0);
        for (std::size_t j = 0; j < n; ++j) rs += P(i, j);
        if (std::abs(num_traits<T>::to_double(rs) - 1.0) > 1e-8)
            throw InputError(
                "smp_passage_moments: the embedded transition matrix rows must sum to one");
    }
    if (hmom.rows() != n)
        throw InputError("smp_passage_moments: hmom must carry one row per state");
    if (hmom.cols() < nmax)
        throw InputError(
            "smp_passage_moments: hmom must carry at least nmax holding-time moments per state");
    const std::vector<std::size_t> tgt =
        passage_detail::unique_target(target, n, "smp_passage_moments");
 
    std::vector<bool> is_target(n, false);
    for (std::size_t k : tgt) is_target[k] = true;
    std::vector<std::size_t> A;
    for (std::size_t i = 0; i < n; ++i)
        if (!is_target[i]) A.push_back(i);
    const std::size_t nA = A.size();
 
    PassageMoments<T> out;
    out.mall = Matrix<T>(n, nmax);
    out.m.clear();  // the reference's `m = []` when pi0 is not given
    if (nA == 0) {
        out.m.assign(nmax, T(0));  // the reference's `m = zeros(1,nmax)` here, pi0 or not
        return out;
    }
 
    // Eq. 8, per state.
    Matrix<T> u(nA, nmax);
    for (std::size_t r = 1; r <= nmax; ++r) {
        for (std::size_t a = 0; a < nA; ++a) {
            T acc = T(0);
            for (std::size_t j = 1; j <= r; ++j) {
                const T base = (r - j == 0) ? T(1) : u(a, r - j - 1);
                double c = 1.0;
                for (std::size_t q = 0; q < j; ++q)
                    c = c * double(r - q) / double(q + 1);
                acc += num_traits<T>::from_double(c) * hmom(A[a], j - 1) * base;
            }
            u(a, r - 1) = -acc;
        }
    }
 
    Matrix<T> IPAA(nA, nA);
    for (std::size_t a = 0; a < nA; ++a)
        for (std::size_t c = 0; c < nA; ++c)
            IPAA(a, c) = (a == c ? T(1) : T(0)) - P(A[a], A[c]);
 
    Matrix<T> M(nA, nmax);
    for (std::size_t q = 1; q <= nmax; ++q) {
        std::vector<T> b(nA, T(0));
        for (std::size_t r = 1; r <= q; ++r) {
            double c = 1.0;
            for (std::size_t k = 0; k < r; ++k) c = c * double(q - k) / double(k + 1);
            for (std::size_t a = 0; a < nA; ++a) {
                const T base = (r < q) ? M(a, q - r - 1) : T(1);
                b[a] -= num_traits<T>::from_double(c) * u(a, r - 1) * base;
            }
        }
        const std::vector<T> x = solve(IPAA, b);
        for (std::size_t a = 0; a < nA; ++a) M(a, q - 1) = x[a];
    }
 
    for (std::size_t a = 0; a < nA; ++a)
        for (std::size_t q = 0; q < nmax; ++q) out.mall(A[a], q) = M(a, q);
    passage_detail::smp_weighted_moments(out, pi0, n, nmax, "smp_passage_moments");
    return out;
}
 
/**
 * Moments of the semi-Markov first passage time from the PER-PAIR kernel
 * moments, Eq. 6: hmom[i][k] holds the moments of the sojourn in state i when
 * the next state is k (an empty entry is a zero vector), so the recursion runs
 * over the kernel r_ik(t) = P(i,k) h_ik(t) directly,
 *
 *     M_i(q) = sum_{k not in B} sum_{r=1..q} C(q,r) m_ik(r) M_k(q-r)
 *            + sum_{k in B} m_ik(q),        M_k(0) = 1.
 *
 * This is the reference's `iscell(hmom)` branch; the (nstates x nmax) matrix
 * form above is its `else` branch.
 *
 * @param hmom (nstates x nstates) table of moment vectors; hmom[i][k][r-1] is
 *             the r-th moment of the sojourn in i given the next state k, and
 *             the entries are expected to already carry the P(i,k) weight, as in
 *             the reference
 */
template <class T>
PassageMoments<T> smp_passage_moments(const Matrix<T>& P,
                                      const std::vector<std::vector<std::vector<T>>>& hmom,
                                      const std::vector<T>& pi0,
                                      const std::vector<std::size_t>& target,
                                      std::size_t nmax = 1) {
    const std::size_t n = P.rows();
    if (P.cols() != n)
        throw InputError("smp_passage_moments: the embedded transition matrix must be square");
    if (nmax == 0) throw InputError("smp_passage_moments: nmax must be positive");
    for (std::size_t i = 0; i < n; ++i) {
        T rs = T(0);
        for (std::size_t j = 0; j < n; ++j) rs += P(i, j);
        if (std::abs(num_traits<T>::to_double(rs) - 1.0) > 1e-8)
            throw InputError(
                "smp_passage_moments: the embedded transition matrix rows must sum to one");
    }
    if (hmom.size() != n)
        throw InputError(
            "smp_passage_moments: a cell HMOM must be (nstates x nstates), entry (i,k) holding "
            "the moments of the sojourn in i when the next state is k");
    for (std::size_t i = 0; i < n; ++i)
        if (hmom[i].size() != n)
            throw InputError(
                "smp_passage_moments: a cell HMOM must be (nstates x nstates), entry (i,k) "
                "holding the moments of the sojourn in i when the next state is k");
    const std::vector<std::size_t> tgt =
        passage_detail::unique_target(target, n, "smp_passage_moments");
 
    std::vector<bool> is_target(n, false);
    for (std::size_t k : tgt) is_target[k] = true;
    std::vector<std::size_t> A;
    for (std::size_t i = 0; i < n; ++i)
        if (!is_target[i]) A.push_back(i);
    const std::size_t nA = A.size();
 
    PassageMoments<T> out;
    out.mall = Matrix<T>(n, nmax, T(0));
    out.m.clear();
    if (nA == 0) {
        out.m.assign(nmax, T(0));  // the reference's `m = zeros(1,nmax)` here, pi0 or not
        return out;
    }
 
    // the reference's local_kernel_moments: an empty entry is zeros(1, nmax)
    const auto mik = [&](std::size_t i, std::size_t k, std::size_t r) {
        const std::vector<T>& v = hmom[i][k];
        return r - 1 < v.size() ? v[r - 1] : T(0);
    };
 
    Matrix<T> IPAA(nA, nA);
    for (std::size_t a = 0; a < nA; ++a)
        for (std::size_t c = 0; c < nA; ++c)
            IPAA(a, c) = (a == c ? T(1) : T(0)) - P(A[a], A[c]);
 
    Matrix<T> M(nA, nmax, T(0));
    for (std::size_t q = 1; q <= nmax; ++q) {
        std::vector<T> b(nA, T(0));
        for (std::size_t a = 0; a < nA; ++a) {
            const std::size_t i = A[a];
            T acc = T(0);
            for (std::size_t c = 0; c < nA; ++c) {
                const std::size_t k = A[c];
                for (std::size_t r = 1; r <= q; ++r) {
                    double bin = 1.0;
                    for (std::size_t t = 0; t < r; ++t) bin = bin * double(q - t) / double(t + 1);
                    const T base = (r < q) ? M(c, q - r - 1) : T(1);  // M_k(0) = 1
                    acc += num_traits<T>::from_double(bin) * mik(i, k, r) * base;
                }
            }
            for (std::size_t k : tgt) acc += mik(i, k, q);
            b[a] = acc;
        }
        const std::vector<T> x = solve(IPAA, b);
        for (std::size_t a = 0; a < nA; ++a) M(a, q - 1) = x[a];
    }
 
    for (std::size_t a = 0; a < nA; ++a)
        for (std::size_t q = 0; q < nmax; ++q) out.mall(A[a], q) = M(a, q);
    passage_detail::smp_weighted_moments(out, pi0, n, nmax, "smp_passage_moments");
    return out;
}
 
/**
 * L(s) of the semi-Markov first passage time, Eqs. 4-5:
 *
 *     L_i(s) = sum_{k not in B} r*_ik(s) L_k(s) + sum_{k in B} r*_ik(s),
 *
 * so (I - R*_AA(s)) L_A(s) = R*_AB(s) 1, one linear system per value of s.
 *
 * @param hlst per-state sojourn transforms h*_i(s), so r*_ik(s) = P(i,k) h*_i(s)
 *             and the complex numbers stay on the DIAGONAL of the system
 */
template <class T>
std::vector<std::complex<double>> smp_passage_lst(
    const Matrix<T>& P, const std::vector<std::function<std::complex<double>(std::complex<double>)>>& hlst,
    const std::vector<T>& pi0, const std::vector<std::size_t>& target,
    const std::vector<std::complex<double>>& s) {
    const std::size_t n = P.rows();
    const std::vector<std::size_t> tgt = passage_detail::unique_target(target, n, "smp_passage_lst");
    if (hlst.size() != n)
        throw InputError("smp_passage_lst: hlst must carry one transform per state");
    std::vector<bool> is_target(n, false);
    for (std::size_t k : tgt) is_target[k] = true;
    std::vector<std::size_t> A;
    for (std::size_t i = 0; i < n; ++i)
        if (!is_target[i]) A.push_back(i);
    const std::size_t nA = A.size();
 
    using C = std::complex<double>;
    double atom = 0.0;
    if (pi0.size() == n)
        for (std::size_t k : tgt) atom += num_traits<T>::to_double(pi0[k]);
 
    std::vector<C> out(s.size(), C(0.0, 0.0));
    for (std::size_t is = 0; is < s.size(); ++is) {
        std::vector<C> h(nA);
        for (std::size_t a = 0; a < nA; ++a) h[a] = hlst[A[a]](s[is]);
        Matrix<C> M(nA, nA);
        std::vector<C> b(nA, C(0.0, 0.0));
        for (std::size_t a = 0; a < nA; ++a) {
            for (std::size_t c = 0; c < nA; ++c)
                M(a, c) = (a == c ? C(1.0, 0.0) : C(0.0, 0.0)) -
                          h[a] * C(num_traits<T>::to_double(P(A[a], A[c])), 0.0);
            double pb = 0.0;
            for (std::size_t k : tgt) pb += num_traits<T>::to_double(P(A[a], k));
            b[a] = h[a] * C(pb, 0.0);
        }
        const std::vector<C> x = passage_detail::solve_cplx(M, b, "smp_passage_lst");
        C acc(atom, 0.0);
        for (std::size_t a = 0; a < nA; ++a)
            acc += C(pi0.size() == n ? num_traits<T>::to_double(pi0[A[a]]) : 1.0 / double(n), 0.0) *
                   x[a];
        out[is] = acc;
    }
    return out;
}
 
/**
 * L(s) of the semi-Markov first passage time from the (n x n) KERNEL of
 * transforms, the reference's non-`perState` branch of Eq. 5: hlst[i][k] is
 * r*_ik(s), the transform of the sojourn in i jointly with the next state being
 * k (so it already carries P(i,k)); an empty entry is zero. Then
 *
 *     (I - R*_AA(s)) L_A(s) = R*_AB(s) 1,
 *
 * with R*(s) the matrix of those transforms, and P is read only for its size.
 *
 * @param hlst (nstates x nstates) table of transforms, empty where the kernel
 *             has no mass
 */
template <class T>
std::vector<std::complex<double>> smp_passage_lst(
    const Matrix<T>& P,
    const std::vector<std::vector<std::function<std::complex<double>(std::complex<double>)>>>& hlst,
    const std::vector<T>& pi0, const std::vector<std::size_t>& target,
    const std::vector<std::complex<double>>& s) {
    const std::size_t n = P.rows();
    const std::vector<std::size_t> tgt = passage_detail::unique_target(target, n, "smp_passage_lst");
    if (hlst.size() != n)
        throw InputError("smp_passage_lst: a kernel hlst must be (nstates x nstates)");
    for (std::size_t i = 0; i < n; ++i)
        if (hlst[i].size() != n)
            throw InputError("smp_passage_lst: a kernel hlst must be (nstates x nstates)");
    if (!pi0.empty() && pi0.size() != n)
        throw InputError(
            "smp_passage_lst: PI0 must be a distribution over the state space, one entry per "
            "state");
    std::vector<bool> is_target(n, false);
    for (std::size_t k : tgt) is_target[k] = true;
    std::vector<std::size_t> A;
    for (std::size_t i = 0; i < n; ++i)
        if (!is_target[i]) A.push_back(i);
    const std::size_t nA = A.size();
 
    using C = std::complex<double>;
    // the reference's `if isempty(pi0), pi0 = ones(1,n)/n`
    std::vector<double> p0(n, 1.0 / double(n));
    if (pi0.size() == n)
        for (std::size_t i = 0; i < n; ++i) p0[i] = num_traits<T>::to_double(pi0[i]);
    double atom = 0.0;
    for (std::size_t k : tgt) atom += p0[k];
 
    std::vector<C> out(s.size(), C(0.0, 0.0));
    for (std::size_t is = 0; is < s.size(); ++is) {
        Matrix<C> M(nA, nA);
        std::vector<C> b(nA, C(0.0, 0.0));
        for (std::size_t a = 0; a < nA; ++a) {
            for (std::size_t c = 0; c < nA; ++c) {
                const C r = hlst[A[a]][A[c]] ? hlst[A[a]][A[c]](s[is]) : C(0.0, 0.0);
                M(a, c) = (a == c ? C(1.0, 0.0) : C(0.0, 0.0)) - r;
            }
            for (std::size_t k : tgt)
                if (hlst[A[a]][k]) b[a] += hlst[A[a]][k](s[is]);
        }
        const std::vector<C> x = passage_detail::solve_cplx(M, b, "smp_passage_lst");
        C acc(atom, 0.0);
        for (std::size_t a = 0; a < nA; ++a) acc += C(p0[A[a]], 0.0) * x[a];
        out[is] = acc;
    }
    return out;
}
 
/**
 * CDF and density of the semi-Markov first passage time, by inverting
 * `smp_passage_lst` through api/lti.
 *
 * There is no matrix-exponential route here: a semi-Markov chain has no
 * generator to exponentiate, which is exactly the case uniformization does not
 * reach and the transform does.
 *
 * `lti_method` defaults to "euler" RATHER THAN "weeks". Semi-Markov passage
 * densities are the case Sec. 4.2 singles out as slow-converging for a Laguerre
 * series, and `laplace_weeks_scaling` then refuses by name rather than
 * returning noise.
 */
template <class T>
PassageCurve<T> smp_passage_time(
    const Matrix<T>& P, const std::vector<std::function<std::complex<double>(std::complex<double>)>>& hlst,
    const std::vector<T>& pi0, const std::vector<std::size_t>& target,
    const std::vector<double>& tset, const std::string& lti_method = "euler") {
    const std::size_t n = P.rows();
    const std::vector<std::size_t> tgt =
        passage_detail::unique_target(target, n, "smp_passage_time");
    double atom = 0.0;
    if (pi0.size() == n)
        for (std::size_t k : tgt) atom += num_traits<T>::to_double(pi0[k]);
 
    const lti::LaplaceFn L = [&](std::complex<double> s) {
        std::vector<std::complex<double>> sv(1, s);
        return smp_passage_lst(P, hlst, pi0, target, sv)[0];
    };
    const lti::LaplaceMethod lm = lti::laplace_method(lti_method);
    PassageCurve<T> out;
    out.t = tset;
    out.atom = atom;
    out.F = lti::laplace_invert_cdf(L, tset, lm);
    const lti::LaplaceFn Ld = [&](std::complex<double> s) { return L(s) - atom; };
    out.f = lti::laplace_invert_pdf(Ld, tset, lm);
    return out;
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_PASSAGE_H