ctmc_state_space_logsize.h

Go to the documentation of this file.

#ifndef LINE_API_MC_CTMC_STATE_SPACE_LOGSIZE_H
#define LINE_API_MC_CTMC_STATE_SPACE_LOGSIZE_H
 
/**
 * @file
 * @ingroup api_mc
 * Worst-case log-size of the CTMC state space induced by a NetworkStruct.
 *
 * Port of `matlab/src/api/mc/ctmc_state_space_logsize.m`. The estimate is the
 * product of five factors, summed in log space:
 *
 *   1. job placements: stars-and-bars C(n_k+M-1, M-1) per class over the
 *      stations that keep no ordered buffer, with an open class holding AT MOST
 *      cutoff jobs over the non-Source stations;
 *   2. buffer orderings: a station outside the share family keeps the CLASS
 *      SEQUENCE of the jobs it holds, so with K>1 classes one occupancy vector
 *      is as many states as its sequences;
 *   3. service phases: the phase count raised to the number of jobs that can be
 *      in service concurrently at the station;
 *   4. routing state: one pointer over the outgoing links per (node, class)
 *      routed RROBIN or WRROBIN;
 *   5. cache content: the ordered items held by each Cache node, times its
 *      retrieval-system subsets; its in-transit job slot is summed with the
 *      placements, taking a closed job out of its chain.
 *
 * It is computed in LOG space throughout because the quantity it exists to
 * detect overflows a double: the `intractableCTMC` fixture (8 PS stations,
 * N=400, Erlang-5) sits at exp(200), and a linear-space estimator would report
 * `inf` for everything above exp(709) and lose the ability to rank one
 * intractable model against another.
 *
 * This is the quantity fed to `ctmc_memory_gate`. It is separate from the gate
 * so a caller such as SolverAUTO can screen CTMC out of a ranking without
 * building the chain.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <limits>
#include <utility>
#include <vector>
 
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
 
namespace line {
namespace mc {
 
 
 
/** Largest (m_1..m_K) box the exact ordered-buffer DP will walk. */
inline constexpr double kOrderGridMax = 1.0e6;
 
/**
 * Log count of (placement, ordering) configurations over ALL order-preserving
 * stations at once, POPULATION CONSERVED. `caps_per[a][k]` bounds buffered class
 * k at ordered station a, `cap_tot[a]` bounds the buffer TOTAL there (a finite
 * station capacity is a slot count, not a per-class bound).
 *
 * The DP state is the population still to place per POPULATION UNIT, not per
 * class: `nj_u[u]` is unit u's population, `b_unit_of[k]` the unit of buffered
 * class k, and unit u spreads what the ordered stations leave over `m_rem_u[u]`
 * share (station, class) slots. A unit is one class, or a closed chain whose
 * jobs switch class, whose N jobs may then take any class mix at a station.
 *
 * Cutoff truncates an OPEN class's population in the network exactly as the
 * plain stars-and-bars term treats it, so open classes are conserved too.
 */
inline double log_ordered_joint(const std::vector<std::vector<int>>& caps_per,
                                const std::vector<double>& cap_tot,
                                const std::vector<int>& nj_u,
                                const std::vector<std::size_t>& b_unit_of,
                                const std::vector<double>& m_rem_u,
                                const std::vector<std::vector<double>>& log_ph_per = {},
                                const std::vector<std::vector<int>>& lev_per = {}) {
    const std::size_t Kb = b_unit_of.size();
    const std::size_t U = nj_u.size();
    std::vector<std::size_t> dims(U);
    std::size_t nstate = 1;
    for (std::size_t u = 0; u < U; ++u) { dims[u] = static_cast<std::size_t>(nj_u[u]) + 1; nstate *= dims[u]; }
    const double NEG = -std::numeric_limits<double>::infinity();
    auto idx_of = [&](const std::vector<int>& v) {
        std::size_t ix = 0, mult = 1;
        for (std::size_t u = 0; u < U; ++u) { ix += static_cast<std::size_t>(v[u]) * mult; mult *= dims[u]; }
        return ix;
    };
    auto sub_of = [&](std::size_t ix) {
        std::vector<int> v(U);
        for (std::size_t u = 0; u < U; ++u) { v[u] = static_cast<int>(ix % dims[u]); ix /= dims[u]; }
        return v;
    };
    std::vector<double> L(nstate, NEG);
    L[idx_of(nj_u)] = 0.0;
    for (std::size_t a = 0; a < caps_per.size(); ++a) {
        std::vector<double> Ln(nstate, NEG);
        // lev_per[a][k] is the priority level of buffered class k at a station that orders jobs only WITHIN a level
        // (LCFSPRIO, State priority-canonical buffers); an empty row marks a station that orders all its jobs.
        const bool by_level = a < lev_per.size() && lev_per[a].size() == Kb;
        std::vector<int> ulev;
        if (by_level) {
            ulev = lev_per[a];
            std::sort(ulev.begin(), ulev.end());
            ulev.erase(std::unique(ulev.begin(), ulev.end()), ulev.end());
        }
        for (std::size_t si = 0; si < nstate; ++si) {
            if (!std::isfinite(L[si])) continue;
            const std::vector<int> rem = sub_of(si);
            std::vector<int> av(Kb);
            for (std::size_t k = 0; k < Kb; ++k) av[k] = std::min(caps_per[a][k], rem[b_unit_of[k]]);
            std::vector<int> m(Kb, 0);
            for (;;) {
                int t = 0;
                for (std::size_t k = 0; k < Kb; ++k) t += m[k];
                if (!(std::isfinite(cap_tot[a]) && t > cap_tot[a])) {
                    double v = L[si];
                    if (!by_level) {
                        v += std::lgamma(static_cast<double>(t) + 1.0);
                        for (std::size_t k = 0; k < Kb; ++k) v -= std::lgamma(static_cast<double>(m[k]) + 1.0);
                    } else if (t > 0) {
                        // One job of SOME present class c is in service and the other t-1 wait in the canonical
                        // layout, ordered only within a level: sum over c of prod_l (t_l(m-e_c))! / prod_k (m-e_c)_k!.
                        // Counting the waiting line alone (no in-service choice) gave 49 against 85 enumerated on a
                        // two-class open M/M/1 at cutoff 6 (MATLAB ctmc_logOrderedJoint).
                        std::vector<double> lc;
                        for (std::size_t c = 0; c < Kb; ++c) {
                            if (m[c] <= 0) continue;
                            std::vector<int> mc = m;
                            mc[c] -= 1;
                            double x = 0.0;
                            for (std::size_t k = 0; k < Kb; ++k) x -= std::lgamma(static_cast<double>(mc[k]) + 1.0);
                            for (int l : ulev) {
                                int tl = 0;
                                for (std::size_t k = 0; k < Kb; ++k)
                                    if (lev_per[a][k] == l) tl += mc[k];
                                x += std::lgamma(static_cast<double>(tl) + 1.0);
                            }
                            lc.push_back(x);
                        }
                        double top = lc[0];
                        for (double x : lc) top = std::max(top, x);
                        double acc = 0.0;
                        for (double x : lc) acc += std::exp(x - top);
                        v += top + std::log(acc);
                    }
                    if (a < log_ph_per.size())
                        for (std::size_t k = 0; k < Kb && k < log_ph_per[a].size(); ++k)
                            v += m[k] * log_ph_per[a][k];
                    // a unit's classes draw on ONE remaining population, so a class mix may overdraw it
                    std::vector<int> nx = rem;
                    bool fits = true;
                    for (std::size_t k = 0; k < Kb; ++k) nx[b_unit_of[k]] -= m[k];
                    for (std::size_t u = 0; u < U; ++u) if (nx[u] < 0) fits = false;
                    if (fits) {
                        const std::size_t di = idx_of(nx);
                        if (!std::isfinite(Ln[di])) Ln[di] = v;
                        else { const double mx = std::max(Ln[di], v);
                               Ln[di] = mx + std::log(std::exp(Ln[di] - mx) + std::exp(v - mx)); }
                    }
                }
                std::size_t pos = 0;
                while (pos < Kb && m[pos] == av[pos]) { m[pos] = 0; ++pos; }
                if (pos == Kb) break;
                ++m[pos];
            }
        }
        L.swap(Ln);
    }
    std::vector<double> terms;
    for (std::size_t si = 0; si < nstate; ++si) {
        if (!std::isfinite(L[si])) continue;
        const std::vector<int> rem = sub_of(si);
        double v = L[si];
        bool stranded = false;
        for (std::size_t u = 0; u < U; ++u) {
            const double r = rem[u];
            if (m_rem_u[u] >= 1.0)
                v += std::lgamma(r + m_rem_u[u]) - std::lgamma(r + 1.0) - std::lgamma(m_rem_u[u]);
            else if (rem[u] > 0)
                stranded = true;  // leftovers with no share slot to go to
        }
        if (stranded) continue;
        terms.push_back(v);
    }
    if (terms.empty()) return NEG;
    double top = terms[0];
    for (double v : terms) top = std::max(top, v);
    double acc = 0.0;
    for (double v : terms) acc += std::exp(v - top);
    return top + std::log(acc);
}
 
/** Options the estimator reads; only the cutoff matters. */
struct CtmcSizeOptions {
    double cutoff = -1.0;  ///< < 0 or non-finite = not given, take the solver default
};
 
/**
 * Worst-case log state-space size of `sn`.
 *
 * @param sn  the struct, ideally after `sn_nonmarkov_toph`, since a
 *            non-Markovian service becomes phases the raw struct does not carry
 * @return    natural log of the worst-case state count
 */
template <class T>
double ctmc_state_space_logsize(const qn::NetworkStruct<T>& sn,
                                const CtmcSizeOptions& opt = CtmcSizeOptions()) {
    using lang::RoutingStrategy;
    using lang::SchedStrategy;
 
    const std::size_t M = sn.nstations;
    const std::size_t K = sn.nclasses;
    if (M == 0 || K == 0) return 0.0;
 
    // The analyzer's own default for an open or mixed model, so the estimate is
    // taken against the space that would actually be built.
    double cutoff = opt.cutoff;
    if (!(cutoff > 0.0) || !std::isfinite(cutoff))
        cutoff = std::ceil(std::pow(6000.0, 1.0 / static_cast<double>(M * K)));
 
    // ORDERED BUFFERS, computed EXACTLY -- see the python twin for the measured
    // numbers. Omitting it under-priced gallery_mmap1_multiclass 311x; BOUNDING
    // it instead of computing it double counts and refused a working model.
    const auto is_share = [](SchedStrategy s) {
        return s == SchedStrategy::INF || s == SchedStrategy::PS || s == SchedStrategy::DPS ||
               s == SchedStrategy::GPS || s == SchedStrategy::PSPRIO ||
               s == SchedStrategy::DPSPRIO || s == SchedStrategy::GPSPRIO ||
               s == SchedStrategy::LPS;
    };
    std::vector<bool> is_buffered(K, true);
    std::size_t Kb = 0;
    for (std::size_t k = 0; k < K; ++k) {
        if (k < sn.issignal.size() && sn.issignal[k]) is_buffered[k] = false;
        if (is_buffered[k]) ++Kb;
    }
    // PREEMPT FAMILY. A preempt-family station stores the phase of EVERY job it
    // holds, the buffered ones as [class, phase] pairs, so a class-k job there is
    // phasessz(i,k) states on its own. Pricing only the jobs in service read a
    // FCFSPRPRIO queue with an Erlang-2 class as 6862 states at cutoff 6 against
    // 148863 enumerated. Such a station enters the ordered DP even with a single
    // buffered class, carrying the per-job phase weight there.
    const auto is_preempt = [](SchedStrategy s) {
        return s == SchedStrategy::FCFSPR || s == SchedStrategy::FCFSPI ||
               s == SchedStrategy::LCFSPR || s == SchedStrategy::LCFSPI ||
               s == SchedStrategy::FCFSPRPRIO || s == SchedStrategy::FCFSPIPRIO ||
               s == SchedStrategy::LCFSPRPRIO || s == SchedStrategy::LCFSPIPRIO;
    };
    const auto is_ordered = [&](std::size_t i) {
        const SchedStrategy s = sn.stations[i].sched;
        return is_preempt(s) || (Kb > 1 && !(s == SchedStrategy::EXT || is_share(s)));
    };
    const auto log_phases = [&](std::size_t i, std::size_t k) {
        const double p = static_cast<double>(sn.phasessz_of(i + 1, k + 1));
        return (std::isfinite(p) && p > 1.0) ? std::log(p) : 0.0;
    };
    std::size_t n_ord = 0;
    for (std::size_t i = 0; i < M; ++i)
        if (is_ordered(i)) ++n_ord;
 
    double log_nstates = 0.0;
    std::vector<double> nk_eff(K, 0.0);
    const std::vector<double>& njobs = sn.njobs();
    for (std::size_t k = 0; k < K; ++k)
        nk_eff[k] = (k < njobs.size() && std::isfinite(njobs[k])) ? njobs[k] : cutoff;
 
    const auto place = [&](double nk, double ms) {
        return std::lgamma(1.0 + nk + ms - 1.0) - std::lgamma(1.0 + ms - 1.0) -
               std::lgamma(1.0 + nk);
    };
    // A ZERO per-class capacity means the class is DISABLED at that station, so it
    // never occupies a slot and the placement term must spread it over the
    // stations that admit it, not over all M. ld_whittle_bandwidth disables each
    // of its three PS routes for the other two classes; counting all M=4 priced it
    // at C(9,6)^3 = 592704 states, 7852 GB under the quadratic byte model, and the
    // gate refused a model whose true space is 7^3 = 343 and solves at once.
    const auto disabled_at = [&](std::size_t i, std::size_t k) {
        return i < sn.classcap.size() && k < sn.classcap[i].size() && sn.classcap[i][k] == 0.0;
    };
    // ZERO VISITS IS THE SECOND WAY IN, and it is the state space generator's own
    // rule: a (station, class) pair whose chain visit ratio is zero gets capacity
    // 0 there, because a class the routing never delivers to a station cannot be
    // there. Only an SPN PLACE escapes, keeping its declared initial marking, so a
    // Place is left alone here too. Reading only classcap missed all of it, since
    // that zeroes a pair whose SERVICE RATE is NaN -- the class is disabled at the
    // station, not unrouted to it. A class that exists only as a switch target on
    // some departure arc has a good service rate everywhere and a visit ratio of
    // zero everywhere: a Source -> FCFS Queue -> Sink model with such a class was
    // priced at 184755 states against the 11 the generator enumerates, predicted
    // 17.62 GB, and had the gate refuse a model that solves in 0.05 s.
    // Visit rows are STATEFUL-node indexed, not station indexed. A pair the chain
    // never enters carries floating-point noise (~1e-17) rather than an exact 0,
    // which is the tolerance below, shared with the generator.
    const double kVisitZeroTol = 1.0e-12;
    std::vector<std::vector<bool>> never(M, std::vector<bool>(K, false));
    if (!sn.visits.empty() && !sn.chains.empty()) {
        for (std::size_t i = 0; i < M; ++i) {
            if (sn.stations[i].nodetype == lang::NodeType::Place) continue;
            const std::size_t nd = sn.node_of_station(i + 1);
            if (nd == 0) continue;
            const std::size_t isf1 = sn.stateful_index(nd);
            if (isf1 == 0) continue;
            const std::size_t isf = isf1 - 1;
            for (std::size_t k = 0; k < K; ++k) {
                std::size_t c = sn.nchains;
                for (std::size_t ci = 0; ci < sn.chains.size(); ++ci)
                    if (k < sn.chains[ci].size() && sn.chains[ci][k]) { c = ci; break; }
                if (c >= sn.visits.size()) continue;
                const auto& V = sn.visits[c];
                if (isf >= V.rows() || k >= V.cols()) continue;
                const double v = num_traits<T>::to_double(V(isf, k));
                if (std::isfinite(v) && std::fabs(v) < kVisitZeroTol) never[i][k] = true;
            }
        }
    }
    const auto admits = [&](std::size_t i, std::size_t k) {
        return !disabled_at(i, k) && !never[i][k];
    };
    const auto admitting_count = [&](std::size_t k) {
        std::size_t n = 0;
        for (std::size_t i = 0; i < M; ++i)
            if (admits(i, k)) ++n;
        return n;
    };
    const auto admitting_all = [&](std::size_t k) {
        const std::size_t n = admitting_count(k);
        return n < 1 ? std::size_t(1) : n;
    };
    const auto admitting_rem = [&](std::size_t k) {
        std::size_t n = 0;
        for (std::size_t i = 0; i < M; ++i) {
            if (is_ordered(i)) continue;
            if (admits(i, k)) ++n;
        }
        return n;
    };
    // A class no station admits holds no job anywhere, so its population is zero
    // rather than the cutoff; leaving it at the cutoff made the class occupy
    // ordered buffers it can never enter.
    for (std::size_t k = 0; k < K; ++k)
        if (nk_eff[k] > 0.0 && admitting_count(k) == 0) nk_eff[k] = 0.0;
    // AN OPEN CLASS HOLDS AT MOST nk_eff JOBS, not exactly nk_eff: the generator enumerates 0..cutoff of them. The
    // Source played that slack bin only for a class it emits; a class that exists only as a switch target downstream
    // (a cache's hit and miss classes) has zero visits at the Source and was priced at ONE placement: an open LRU
    // cache feeding a Delay and a PS queue read 24 states against 11616. So an open class takes C(n+M', M') over the
    // M' non-Source stations admitting it, a closed one C(n+M-1, M-1) over all of them (MATLAB ctmc_placement_term).
    const auto placement_term = [&](std::size_t k, double n) {
        const bool is_open = !(k < njobs.size() && std::isfinite(njobs[k]));
        if (!is_open) return place(n, static_cast<double>(admitting_all(k)));
        std::size_t mopen = 0;
        for (std::size_t i = 0; i < M; ++i)
            if (sn.stations[i].sched != SchedStrategy::EXT && admits(i, k)) ++mopen;
        const double mo = static_cast<double>(mopen);
        return std::lgamma(1.0 + n + mo) - std::lgamma(1.0 + mo) - std::lgamma(1.0 + n);
    };
 
    // CLOSED CHAINS THAT SWITCH CLASS conserve their population per CHAIN, not per class: LINE gives the reference
    // class all N jobs and the switch targets none, so per-class terms placed N jobs over the reference class's
    // stations and NONE over the targets'. A closed cache loop (Delay -> Cache -> Delay, hit/miss classes back, 10
    // items capacity 3) read 2880 states against 4320, a closed switching tandem 4 against 10. Such a chain is ONE
    // population unit: its N jobs spread over every admitted (station, class) pair, C(N+P-1, P-1), and at an ordered
    // station the DP sums over the chain's class mixes there. Each class of the chain may hold up to N jobs, which
    // also bounds its phase and capacity terms (MATLAB ctmc_state_space_logsize.m, _kb/06).
    std::vector<std::size_t> unit_of(K, 0);
    std::vector<std::vector<std::size_t>> unit_cls;
    for (std::size_t ci = 0; ci < sn.chains.size(); ++ci) {
        std::vector<std::size_t> cls;
        for (std::size_t k = 0; k < K && k < sn.chains[ci].size(); ++k)
            if (sn.chains[ci][k] && k < njobs.size() && std::isfinite(njobs[k]) && admitting_count(k) > 0)
                cls.push_back(k);
        if (cls.size() > 1) {
            double ntot = 0.0;
            for (std::size_t k : cls) ntot += njobs[k];
            unit_cls.push_back(cls);
            for (std::size_t k : cls) { nk_eff[k] = ntot; unit_of[k] = unit_cls.size(); }
        }
    }
    for (std::size_t k = 0; k < K; ++k)
        if (unit_of[k] == 0) { unit_cls.push_back(std::vector<std::size_t>{k}); unit_of[k] = unit_cls.size(); }
    // unit_of is 1-based from here on; unit_cls[unit_of[k] - 1] holds class k
    const auto unit_closed = [&](std::size_t u) {
        const std::size_t k = unit_cls[u][0];
        return k < njobs.size() && std::isfinite(njobs[k]);
    };
 
    // Placements and buffer orderings with per-class populations nk.
    const auto placement_block = [&](const std::vector<double>& nk) -> double {
        double lp = 0.0;
        if (n_ord == 0) {
            for (std::size_t u = 0; u < unit_cls.size(); ++u) {
                const std::vector<std::size_t>& cls = unit_cls[u];
                if (cls.size() == 1) { lp += placement_term(cls[0], nk[cls[0]]); continue; }
                const double n = nk[cls[0]];
                double P = 0.0;
                for (std::size_t k : cls) P += static_cast<double>(admitting_count(k));
                P = std::max(P, 1.0);
                lp += std::lgamma(n + P) - std::lgamma(P) - std::lgamma(n + 1.0);
            }
            return lp;
        }
        for (std::size_t k = 0; k < K; ++k)
            if (!is_buffered[k]) lp += placement_term(k, nk[k]);
        std::vector<int> njb;
        std::vector<std::size_t> bunits;  // 0-based unit ids, first-seen order over the buffered classes
        std::vector<std::size_t> b_unit_of;
        for (std::size_t k = 0; k < K; ++k) {
            if (!is_buffered[k]) continue;
            njb.push_back(static_cast<int>(std::floor(nk[k])));
            const std::size_t u = unit_of[k] - 1;
            std::size_t bu = 0;
            while (bu < bunits.size() && bunits[bu] != u) ++bu;
            if (bu == bunits.size()) bunits.push_back(u);
            b_unit_of.push_back(bu);
        }
        std::vector<int> nj_u;
        std::vector<double> m_rem_u;
        double grid = 1.0;
        for (std::size_t u : bunits) {
            const int n = static_cast<int>(std::floor(nk[unit_cls[u][0]]));
            nj_u.push_back(n);
            grid *= static_cast<double>(n) + 1.0;
            double mr = 0.0;
            for (std::size_t k : unit_cls[u]) mr += static_cast<double>(admitting_rem(k));
            m_rem_u.push_back(mr);
        }
        if (grid <= kOrderGridMax) {
            std::vector<std::vector<int>> caps_per;
            std::vector<std::vector<double>> log_ph_per;
            std::vector<std::vector<int>> lev_per;
            std::vector<double> cap_tot;
            for (std::size_t i = 0; i < M; ++i) {
                if (!is_ordered(i)) continue;
                std::vector<double> lph;
                for (std::size_t k = 0; k < K; ++k)
                    if (is_buffered[k]) lph.push_back(is_preempt(sn.stations[i].sched) ? log_phases(i, k) : 0.0);
                log_ph_per.push_back(lph);
                std::vector<int> lev;
                if (sn.stations[i].sched == SchedStrategy::LCFSPRIO)
                    for (std::size_t k = 0; k < K; ++k)
                        if (is_buffered[k]) lev.push_back(k < sn.classes.size() ? sn.classes[k].prio : 0);
                lev_per.push_back(lev);
                std::vector<int> per;
                std::size_t bi = 0;
                for (std::size_t k = 0; k < K; ++k)
                    if (is_buffered[k]) {
                        int c = njb[bi];
                        if (i < sn.classcap.size() && k < sn.classcap[i].size() &&
                            std::isfinite(sn.classcap[i][k]))
                            c = std::min(c, static_cast<int>(std::floor(sn.classcap[i][k])));
                        if (!admits(i, k)) c = 0;
                        per.push_back(c);
                        ++bi;
                    }
                caps_per.push_back(per);
                const double c = sn.stations[i].cap;
                cap_tot.push_back(std::isfinite(c) && c >= 0
                                      ? std::floor(c) : std::numeric_limits<double>::infinity());
            }
            return lp + log_ordered_joint(caps_per, cap_tot, nj_u, b_unit_of, m_rem_u, log_ph_per, lev_per);
        }
        double total = 0.0;
        for (int c : njb) total += c;
        const std::size_t nb = njb.size();
        if (nb > 1) {  // one class has a single sequence; only a preempt station gets here with it
            const double lkb = std::log(static_cast<double>(nb));
            lp += static_cast<double>(n_ord) *
                  ((total + 1.0) * lkb - std::log(static_cast<double>(nb) - 1.0) +
                   std::log1p(-std::exp(-(total + 1.0) * lkb)));
        }
        for (std::size_t i = 0; i < M; ++i)
            if (is_preempt(sn.stations[i].sched))
                for (std::size_t k = 0; k < K; ++k)
                    if (is_buffered[k]) lp += nk[k] * log_phases(i, k);
        for (std::size_t k = 0; k < K; ++k)
            if (is_buffered[k]) {
                const std::size_t mk = admitting_rem(k);
                if (mk >= 1) lp += place(nk[k], static_cast<double>(mk));
            }
        return lp;
    };
 
    // A CACHE'S IN-TRANSIT SLOT holds no job or one job of a class visiting the cache. A closed job held there is
    // NOT at a station, so it is taken out of its chain's population and the placements are SUMMED over the slot
    // states; multiplying them by (1 + nin) counted that job twice and read the closed cache loop at 2x truth. An
    // open job's truncated population is not reduced, as the generator enumerates it.
    std::vector<std::vector<std::size_t>> transit;  // per Cache node, the visiting classes
    for (const auto& kv : sn.nodeparam) {
        const std::size_t ind = kv.first;
        if (ind == 0 || ind > sn.nodes.size() || sn.nodes[ind - 1].nodetype != lang::NodeType::Cache) continue;
        if (kv.second.nitems == 0) continue;
        std::vector<std::size_t> cls;
        const std::size_t isf1 = sn.stateful_index(ind);
        if (isf1 > 0 && !sn.visits.empty()) {
            for (std::size_t k = 0; k < K; ++k) {
                std::size_t c = sn.chains.size();
                for (std::size_t ci = 0; ci < sn.chains.size(); ++ci)
                    if (k < sn.chains[ci].size() && sn.chains[ci][k]) { c = ci; break; }
                if (c >= sn.visits.size()) continue;
                const auto& V = sn.visits[c];
                if (isf1 - 1 >= V.rows() || k >= V.cols()) continue;
                if (num_traits<T>::to_double(V(isf1 - 1, k)) > kVisitZeroTol) cls.push_back(k);
            }
        }
        transit.push_back(cls);
    }
    // per cache: (unit + 1, weight) options, unit 0 = empty slot of weight 1; weight = #visiting classes of the unit
    std::vector<std::vector<std::pair<std::size_t, double>>> topts;
    double combos = 1.0;
    for (const auto& cls : transit) {
        std::vector<std::pair<std::size_t, double>> o{{0, 1.0}};
        for (std::size_t k : cls) {
            const std::size_t u1 = unit_of[k];
            bool found = false;
            for (auto& e : o)
                if (e.first == u1) { e.second += 1.0; found = true; }
            if (!found) o.push_back({u1, 1.0});
        }
        combos *= static_cast<double>(o.size());
        topts.push_back(o);
    }
    if (transit.empty() || combos > 4096.0) {
        log_nstates += placement_block(nk_eff);
        for (const auto& cls : transit) log_nstates += std::log(1.0 + static_cast<double>(cls.size()));
    } else {
        std::vector<double> vals;
        std::vector<std::size_t> pick(topts.size(), 0);
        for (;;) {
            std::vector<double> nk = nk_eff;
            double lw = 0.0;
            bool ok = true;
            for (std::size_t t = 0; t < topts.size(); ++t) {
                const std::size_t u1 = topts[t][pick[t]].first;
                if (u1 == 0) continue;
                lw += std::log(topts[t][pick[t]].second);
                const std::size_t u = u1 - 1;
                if (unit_closed(u)) {
                    for (std::size_t k : unit_cls[u]) nk[k] -= 1.0;
                    if (nk[unit_cls[u][0]] < 0.0) ok = false;
                }
            }
            if (ok) vals.push_back(lw + placement_block(nk));
            std::size_t pos = 0;
            while (pos < topts.size() && pick[pos] + 1 == topts[pos].size()) { pick[pos] = 0; ++pos; }
            if (pos == topts.size()) break;
            ++pick[pos];
        }
        double top = vals[0];
        for (double v : vals) top = std::max(top, v);
        double acc = 0.0;
        for (double v : vals) acc += std::exp(v - top);
        log_nstates += top + std::log(acc);
    }
 
    // A sharing discipline can hold every job in service at once; a queueing one
    // holds at most its server count.
    for (std::size_t i = 1; i <= M; ++i) {
        const SchedStrategy sched = sn.stations[i - 1].sched;
        const bool share = sched == SchedStrategy::INF || sched == SchedStrategy::PS ||
                           sched == SchedStrategy::DPS || sched == SchedStrategy::GPS ||
                           sched == SchedStrategy::PSPRIO || sched == SchedStrategy::DPSPRIO ||
                           sched == SchedStrategy::GPSPRIO || sched == SchedStrategy::LPS;
        if (is_preempt(sched)) continue;  // every job's phase is already in the ordered term
        for (std::size_t r = 1; r <= K; ++r) {
            const double p = static_cast<double>(sn.phasessz_of(i, r));
            if (!std::isfinite(p) || p <= 1.0 || !admits(i - 1, r - 1)) continue;  // no job, no phase to keep
            double m;
            if (sched == SchedStrategy::EXT)
                m = 1.0;
            else if (share)
                m = nk_eff[r - 1];
            else
                m = std::min(nk_eff[r - 1], sn.stations[i - 1].nservers);
            if (!std::isfinite(m)) m = nk_eff[r - 1];
            log_nstates += std::lgamma(1.0 + m + p - 1.0) - std::lgamma(1.0 + p - 1.0) -
                           std::lgamma(1.0 + m);
        }
    }
 
    // Round-robin routing is stateful: the pointer over the outgoing links is
    // part of the state. The reference reads `sn.connmatrix`; this port has no
    // such field and takes the out-degree from `rtnodes`, the same graph after
    // the refresh resolved the routing strategies.
    for (std::size_t ind = 1; ind <= sn.nodes.size(); ++ind) {
        const std::vector<RoutingStrategy>& rt_i = sn.nodes[ind - 1].routing;
        std::size_t nrr = 0;
        for (std::size_t r = 0; r < K && r < rt_i.size(); ++r)
            if (rt_i[r] == RoutingStrategy::RROBIN || rt_i[r] == RoutingStrategy::WRROBIN) ++nrr;
        if (nrr == 0) continue;
        const std::size_t nout = sn.downstream_stations(ind).size();
        if (nout <= 1) continue;
        log_nstates += static_cast<double>(nrr) * std::log(static_cast<double>(nout));
    }
 
    // 5. cache content: a Cache node keeps the ORDERED placement of C = sum(itemcap) distinct items out of n,
    //    n!/(n-C)! local states, times the subsets of at most retrieval_capacity of the other n-C items in its
    //    retrieval system. Leaving it out priced an open LRU cache (10 items, capacity 3) as if it had one local
    //    state and let a chain 720 times larger run.
    for (const auto& kv : sn.nodeparam) {
        const std::size_t ind = kv.first;
        if (ind == 0 || ind > sn.nodes.size() || sn.nodes[ind - 1].nodetype != lang::NodeType::Cache) continue;
        const qn::CacheParam<T>& cp = kv.second;
        if (cp.nitems == 0) continue;
        const double n = static_cast<double>(cp.nitems);
        double csum = 0.0;
        for (int c : cp.itemcap) csum += c;
        const double C = std::min(csum, n);
        log_nstates += std::lgamma(n + 1.0) - std::lgamma(n - C + 1.0);
        // the job in transit is priced with the placements above, summed over the slot states
        if (cp.retrieval_capacity > 0) {
            const double rest = n - C;
            const double smax = std::min(static_cast<double>(cp.retrieval_capacity), rest);
            std::vector<double> lt;
            for (double sv = 0.0; sv <= smax; sv += 1.0)
                lt.push_back(std::lgamma(rest + 1.0) - std::lgamma(sv + 1.0) - std::lgamma(rest - sv + 1.0));
            double top = lt[0];
            for (double v : lt) top = std::max(top, v);
            double acc = 0.0;
            for (double v : lt) acc += std::exp(v - top);
            log_nstates += top + std::log(acc);
        }
    }
 
    return log_nstates;
}
 
}  // namespace mc
}  // namespace line
 
#endif  // LINE_API_MC_CTMC_STATE_SPACE_LOGSIZE_H