solver_ssa_nrm.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_SSA_SOLVER_SSA_NRM_H
#define LINE_SOLVERS_SSA_SOLVER_SSA_NRM_H
 
/**
 * @file
 * @ingroup line_solvers
 * SolverSSA, the `nrm` method: a port of `solver_ssa_nrm.m` and its analyzer
 * `solver_ssa_analyzer_nrm.m`.
 *
 * WHAT THE METHOD IS. The queueing network is rewritten as a REACTION NETWORK
 * over a state vector counting jobs per (node, class, service phase), and the
 * sample path is generated by Gibson and Bruck's Next Reaction Method: every
 * reaction carries an absolute exponential clock `Pk`, the elapsed integrated
 * propensity `Tk`, and fires at `tau = (Pk - Tk)/Ak`; only the reactions in the
 * fired reaction's dependency set have their propensity recomputed. Metrics are
 * accumulated as time integrals along the path rather than from a stored
 * trajectory, so the memory cost is independent of the sample count.
 *
 * NRM IS NOT THE SERIAL SIMULATOR. `method='serial'` is a different engine with
 * a different random stream and a different event granularity. A seed-fixed
 * result from one cannot be reproduced by the other, and a mismatch between
 * them is not evidence of a defect in either.
 *
 * DOUBLE ONLY. A sample path is generated from exponential clocks, which are
 * `-log(u)` of a uniform; there is no exact rational value to compute and no
 * additional information a wider float could carry, since the answer's error is
 * the Monte Carlo error and not the rounding. A non-`double` backend is
 * therefore refused BY NAME rather than silently narrowed, exactly as
 * `solver_fluid` refuses one.
 *
 * WHAT IS PORTED AND WHAT REFUSES lives in `ssa_dispatch.h`, which holds the
 * eligibility gate. This file assumes the gate has already passed and refuses
 * again, by name, for anything it meets that the gate should have caught -- a
 * second gate is cheap and a silently wrong sample path is not.
 *
 * ONE DELIBERATE DIVERGENCE, and it is a reference defect. In
 * `solver_ssa_nrm.m` lines 1511-1521 the finite-capacity loss test is applied
 * to the destination of EVERY single-destination firing, with only the renege
 * and retry columns excluded -- `isPhaseRx` is not in the exclusion list. A
 * PHASE CHANGE is a single-destination firing whose "destination" is another
 * phase slot of the same job at the same station, so at a station that is at
 * its physical capacity with an open class the reference declares the phase
 * change a loss and DESTROYS the job. The condition is reachable (an M/PH/1/K
 * queue under PS, which `phaseNrmOK` admits). This port excludes phase changes
 * from the loss test. The divergence is stated rather than hidden because
 * reproducing it would silently violate job conservation, and it is not fixed
 * in the reference from here because that file is outside this port's scope.
 */
 
#include "line/util/line_console.h"
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <vector>
 
#include "line/lang/distribution.h"
#include "line/lang/qn/network_struct.h"
#include "line/lang/qn/polling_info.h"
#include "line/lang/qn/state_events.h"
#include "line/solvers/ssa/ssa_types.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace ssa {
 
namespace detail {
 
using lang::SchedStrategy;
 
/** True for the non-preemptive disciplines that hold waiting jobs in a buffer. */
inline bool sched_is_buffered(SchedStrategy s) {
    return s == SchedStrategy::FCFS || s == SchedStrategy::LCFS || s == SchedStrategy::SIRO ||
           s == SchedStrategy::HOL || s == SchedStrategy::SEPT || s == SchedStrategy::LEPT ||
           s == SchedStrategy::LCFSPR;
}
 
/**
 * The pass-and-swap family, whose buffer is the FULL ordered job list (oldest
 * first) rather than the waiting jobs: `isListSched`. The reference canonicalizes
 * OI to PAS with an all-zero swap graph; here OI keeps its own tag and
 * `sn.pasparam` carries the zero graph, so both read the same machinery.
 */
inline bool sched_is_list(SchedStrategy s) {
    return s == SchedStrategy::PAS || s == SchedStrategy::OI;
}
 
/**
 * Of those, the ones whose in-service phase composition the NRM can track.
 *
 * LCFSPR is excluded: preempt-resume would need the preempted job's phase
 * remembered in the buffer, which the waiting-only buffer does not record.
 * `phaseNrmOK` in the reference draws the same line.
 */
inline bool sched_is_buf_ph(SchedStrategy s) {
    return s == SchedStrategy::FCFS || s == SchedStrategy::LCFS || s == SchedStrategy::SIRO ||
           s == SchedStrategy::HOL || s == SchedStrategy::SEPT || s == SchedStrategy::LEPT;
}
 
/**
 * The DPS / GPS family, whose rate law reads per-class weights from schedparam.
 *
 * THE PRIORITY VARIANTS BELONG HERE TOO. `dpsprioshare` and `gpsprioshare` call
 * straight through to `dpsshare` / `gpsshare` -- with the priority-restricted
 * population above capacity and the full one below it -- so they read the same
 * weights. Leaving them out left `wnorm_` all zeros at such a station, which
 * makes every share zero and every reaction zero-rate: the engine reports a
 * DEADLOCK rather than a wrong number, which is at least loud.
 */
inline bool sched_is_weighted(SchedStrategy s) {
    return s == SchedStrategy::DPS || s == SchedStrategy::GPS ||
           s == SchedStrategy::DPSPRIO || s == SchedStrategy::GPSPRIO;
}
 
/** One reaction of the network: a departure, or a phase change within a service process. */
struct NrmReaction {
    std::size_t node = 0;   ///< 0-based node the reaction consumes a job at
    std::size_t cls = 0;    ///< 0-based class
    std::size_t from = 0;   ///< 0-based state slot the propensity reads
    double rate = 0.0;      ///< the service-process rate the reaction carries
    bool is_phase = false;  ///< an internal D0 phase change rather than a departure
    std::size_t phase_from = 0, phase_to = 0;  ///< 0-based phases of a phase change
    std::size_t dep_phase = 0;                 ///< 0-based phase a departure absorbs from
    bool is_buf_svc = false;                   ///< departure of a buffered-PH class
    /** A polling SWITCHOVER step: moves no job, reads and advances the controller. */
    bool is_sw = false;
    /** A cache READ: the outcome class and the contents update are drawn at firing. */
    bool is_cache = false;
 
    /** `find(P(:,k))` and its cumulative weights, when the destination is drawn. */
    std::vector<std::size_t> to_slots;
    std::vector<double> cdf;
    /** `find(S(:,k) < 0)`, the slots the firing decrements. */
    std::vector<std::size_t> from_slots;
    std::size_t nnzP = 0;
    /** `find(S(:,k) > 0)` when the destination is deterministic; npos = none. */
    std::size_t det_dest = static_cast<std::size_t>(-1);
 
    /**
     * A STATE-DEPENDENT destination, resolved at firing time rather than drawn
     * from `cdf`.
     *
     * `refresh_routing` expands RROBIN, WRROBIN, JSQ and SQ into a probability
     * split so that the matrix solvers have a matrix to read, and that split is
     * what `cdf` above would carry. It is the right answer for a mean-value
     * solver and the WRONG one for a sample path: a round-robin dispatcher
     * sends the k-th job to the k-th target, not to a uniformly drawn one, and
     * that is precisely the variance the strategy exists to remove. So a
     * reaction leaving a node that declares one of the four keeps the CANDIDATE
     * list instead, in the order `refresh_routing` enumerates it, and picks
     * from the live population.
     */
    lang::RoutingStrategy sd_strategy = lang::RoutingStrategy::PROB;
    /** One candidate per declared out-arc: the destination (node, class). */
    std::vector<std::size_t> sd_node, sd_class;
    /** WRROBIN weights, normalised; empty for the other three. */
    std::vector<double> sd_weight;
    /** Per-candidate entry-phase slots and their `pie` weights. */
    std::vector<std::vector<std::size_t>> sd_slots;
    std::vector<std::vector<double>> sd_pie;
};
 
}  // namespace detail
 
/**
 * The NRM engine: the reaction network built from `sn`, and the sample path.
 *
 * It is a class rather than a chain of free functions because the run loop
 * threads eight pieces of mutable state (the population vector, the per-node
 * buffers, the in-service phase multisets, the clocks, the propensities and the
 * three metric accumulators) through every step, and the reference threads the
 * same eight as arguments to a single 700-line function.
 */
template <class T>
class NrmEngine {
public:
    NrmEngine(const qn::NetworkStruct<T>& sn, const SsaOptions& opt)
        : sn_(sn), opt_(opt), rng_(opt.seed) {
        build_layout();
        check_peaks();
        build_regions();
        build_reactions();
        rr_cursor_.assign(I_ * K_ + K_ + 1, 0);
        build_dependencies();
        build_initial_state();
    }
 
    /** Run `opt.samples` firings and return the time-averaged metrics. */
    SsaSolution run();
 
    /** The state-vector length, for the tests that assert the phase expansion. */
    std::size_t nstates() const { return NS_; }
    /** The reaction count, likewise. */
    std::size_t nreactions() const { return rx_.size(); }
    /** The cache write-back of the last `run()`, one entry per Cache node. */
    const std::vector<SsaCacheRatio>& cache() const { return cache_out_; }
 
private:
    using SchedStrategy = lang::SchedStrategy;
    using Rx = detail::NrmReaction;
 
    // ---- the model ---------------------------------------------------------
    const qn::NetworkStruct<T>& sn_;
    SsaOptions opt_;
    SsaRng rng_;
 
    std::size_t I_ = 0, K_ = 0, M_ = 0, NS_ = 0, maxnph_ = 1;
 
    // ---- the (node, class, phase) slot layout ------------------------------
    std::vector<std::vector<std::size_t>> phoff_, nph_;
    std::vector<std::size_t> slot_node_, slot_class_;
 
    // ---- per-node and per-station data read by the rate laws ---------------
    std::vector<bool> is_station_;
    std::vector<std::size_t> to_station_;  ///< 0-based station, npos when not a station
    std::vector<double> mi_;               ///< servers, infinite off a station
    std::vector<std::vector<double>> rate_;
    std::vector<std::vector<std::vector<double>>> pie_;  ///< entry-phase law per (node, class)
    std::vector<std::vector<bool>> buf_ph_class_;
    std::vector<bool> buf_ph_node_;
    std::vector<SchedStrategy> sched_;
    std::vector<std::vector<double>> lld_;    ///< per station, empty when not load dependent
    std::vector<std::vector<double>> wnorm_;  ///< per station, normalized DPS/GPS weights
    std::vector<double> classprio_;
    std::vector<double> nservers_;
 
    // ---- the reaction network ---------------------------------------------
    std::vector<Rx> rx_;
    Matrix<double> S_;
    std::vector<std::vector<std::size_t>> D_;
    /** Departure reactions of each (station, class), for the throughput integral. */
    std::vector<std::vector<std::vector<std::size_t>>> dep_rx_;
    /**
     * Stations and classes at which a refused arrival must BLOCK rather than be
     * lost, and whether any exist at all. See `capacity_block`.
     */
    std::vector<std::vector<char>> blk_can_;
    bool blk_on_ = false;
 
    // ---- the sample path ---------------------------------------------------
    std::vector<double> nvec0_;
    std::vector<std::vector<std::size_t>> buffers0_;
    std::vector<Matrix<double>> svcph0_;
    /**
     * The derived START and PREEMPT tags, as COUNTS of events along the path;
     * `run` divides them by the simulated time. This engine has no
     * `EventOutcome` to annotate -- it is a reaction network, not a state
     * machine -- so the tags are raised where the buffer moves, which is the
     * only place a job can take or lose a server here. Both engines must carry
     * them or the counters would depend on `options.method`.
     */
    Matrix<double> start_cnt_, preempt_cnt_;
    /**
     * Departures that were BLOCKED, per (station, class). `TN` integrates the
     * PROPENSITY, which counts a departure the station never makes once its
     * successor is full (0.744 against the exact 0.652 on the BUG-81 tandem), so
     * the blocked firings are subtracted from that integral before it is
     * normalized. In expectation the count IS the integral of the blocked share
     * of the rate, so the difference is unbiased -- and unlike recomputing that
     * share it needs no second evaluation of a state-dependent dispatcher, whose
     * draw would otherwise have to be replayed.
     */
    Matrix<double> block_cnt_;
 
    /**
     * One finite capacity region, `fcrPrecompute`: its member NODES, so the gate
     * reads the population straight off the state vector, the caps reduced over
     * the member rows (infinity = unbounded), and the per-class rule.
     */
    struct FcrRegion {
        std::vector<char> member_node;  ///< per node
        std::vector<double> ccap;       ///< per class
        double gcap = std::numeric_limits<double>::infinity();
        double memcap = std::numeric_limits<double>::infinity();
        std::vector<double> sz;         ///< per-class memory footprint
        Matrix<double> A;               ///< optional linear constraint A n <= b
        std::vector<double> b;
        std::vector<char> waitq;        ///< per class: WAITQ parks, DROP loses
    };
    std::vector<FcrRegion> fcr_;
    bool fcr_any_waitq_ = false;
 
    /**
     * The polling controller of each node, the reference's `buffers{ind} =
     * [mode, pos, swk, ctr]`: mode 0 parked, 1 serving `pos`, 2 switching
     * towards `pos` in switchover phase `swk` (0-based). `pos` is 1-based, as
     * `qn::polling_next` reads it. Mutable path state, like the buffers, but a
     * member so the propensity can read it without widening its signature.
     */
    // ---- the Cache sub-engine, `solver_ssa_nrm.m` lines 103-160 / cacheAccess --
    /** 1 at a Cache node; `cache_read_[ind][r]` marks the classes that READ there. */
    std::vector<char> is_cache_node_;
    std::vector<std::vector<char>> cache_read_;
    /** The contents row of each cache, laid out as `after_event_cache` reads it. */
    std::vector<std::vector<T>> cache_var0_, cache_var_;
    /** Jobs PRODUCED per (cache node, class), and merges per (cache node, read class). */
    Matrix<double> cache_prod_, cache_dly_;
    std::vector<SsaCacheRatio> cache_out_;
    void cache_fire(const Rx& x, std::vector<double>& X, std::size_t& dest_pos, bool& have_dest);
    void cache_write_back();
 
    struct PollCtrl {
        int mode = 0;
        std::size_t pos = 1;
        std::size_t swk = 0;
        long ctr = 0;
    };
    std::vector<char> is_poll_;
    std::vector<qn::PollingInfo<T>> poll_info_;
    std::vector<PollCtrl> poll_ctrl_;
    /** The controller row `pollLandCtrl` lands in; raises the START of a visit. */
    PollCtrl poll_land(std::size_t ind, std::size_t q, int mode, long budget, bool raise = true);
    /** Per-class populations at node `ind`, as `polling_next` reads them. */
    std::vector<long> poll_nbuf(const std::vector<double>& X, std::size_t ind) const {
        std::vector<long> nb(K_, 0);
        for (std::size_t r = 0; r < K_; ++r) nb[r] = std::lround(class_pop(X, ind, r));
        return nb;
    }
 
    // ---- construction ------------------------------------------------------
    void build_layout();
    void build_regions();
    void build_reactions();
    void build_dependencies();
    void build_initial_state();
 
    // ---- rate-law helpers, one per reference subfunction -------------------
    double class_pop(const std::vector<double>& X, std::size_t ind, std::size_t r) const {
        double s = 0.0;
        for (std::size_t k = 0; k < nph_[ind][r]; ++k) s += X[phoff_[ind][r] + k];
        return s;
    }
    std::vector<double> class_counts(const std::vector<double>& X, std::size_t ind) const {
        std::vector<double> v(K_, 0.0);
        for (std::size_t r = 0; r < K_; ++r) v[r] = class_pop(X, ind, r);
        return v;
    }
    double kir_frac(const std::vector<double>& X, std::size_t slot, std::size_t ind,
                    std::size_t r) const {
        const double nir = class_pop(X, ind, r);
        return nir <= 0.0 ? 0.0 : X[slot] / nir;
    }
    double lldfac(std::size_t ist, double ntot) const {
        if (ist == npos || lld_[ist].empty() || ntot < 1.0) return 1.0;
        const std::size_t lim = lld_[ist].size();
        std::size_t k = static_cast<std::size_t>(std::llround(ntot));
        if (k > lim) k = lim;
        return lld_[ist][k - 1];
    }
    static double dpsshare(const std::vector<double>& w, const std::vector<double>& n,
                           std::size_t r) {
        double den = 0.0;
        for (std::size_t s = 0; s < n.size(); ++s) den += w[s] * n[s];
        return den <= 0.0 ? 0.0 : w[r] * n[r] / den;
    }
    static double gpsshare(const std::vector<double>& w, const std::vector<double>& n,
                           std::size_t r) {
        if (n[r] <= 0.0) return 0.0;
        double den = 0.0;
        for (std::size_t s = 0; s < n.size(); ++s) den += n[s] > 0.0 ? w[s] : 0.0;
        return den <= 0.0 ? 0.0 : w[r] / den;
    }
    /**
     * The three PRIORITY sharing disciplines, ported 2026-08-15.
     *
     * ALL THREE SPLIT AT THE SERVER COUNT, and that is the whole idea: below
     * capacity every job is in service and the station behaves as its
     * non-priority twin, so priority cannot matter; above capacity only the
     * MOST URGENT NON-EMPTY group shares the servers and everyone else is
     * frozen. LINE orders priorities with the LOWER value more urgent.
     *
     * "Non-empty" is doing real work in `urgent`: an empty class never defines
     * the urgent group, so a station holding only low-priority jobs serves them
     * rather than stalling on a priority class that is not there.
     */
    bool urgent(const std::vector<double>& n, std::size_t r) const {
        double best = 0.0;
        bool any = false;
        for (std::size_t s = 0; s < n.size(); ++s)
            if (n[s] > 0.0 && (!any || classprio_[s] < best)) {
                best = classprio_[s];
                any = true;
            }
        return any && classprio_[r] == best;
    }
    /** The population vector restricted to r's priority group, and its total. */
    std::vector<double> prio_group(const std::vector<double>& n, std::size_t r,
                                   double* total = nullptr) const {
        std::vector<double> act(n.size(), 0.0);
        double tot = 0.0;
        for (std::size_t s = 0; s < n.size(); ++s)
            if (classprio_[s] == classprio_[r]) {
                act[s] = n[s];
                tot += n[s];
            }
        if (total) *total = tot;
        return act;
    }
    double psprioshare(const std::vector<double>& n, std::size_t r, double c) const {
        double ni = 0.0;
        for (double v : n) ni += v;
        if (ni <= 0.0) return 0.0;
        if (ni <= c) return (n[r] / ni) * std::min(ni, c);
        if (!urgent(n, r)) return 0.0;
        double niprio = 0.0;
        prio_group(n, r, &niprio);
        return niprio <= 0.0 ? 0.0 : (n[r] / niprio) * std::min(niprio, c);
    }
    double dpsprioshare(const std::vector<double>& w, const std::vector<double>& n, std::size_t r,
                        double c) const {
        double ni = 0.0;
        for (double v : n) ni += v;
        if (ni <= 0.0) return 0.0;
        if (ni <= c) return dpsshare(w, n, r);
        if (!urgent(n, r)) return 0.0;
        return dpsshare(w, prio_group(n, r), r);
    }
    double gpsprioshare(const std::vector<double>& w, const std::vector<double>& n, std::size_t r,
                        double c) const {
        double ni = 0.0;
        for (double v : n) ni += v;
        if (ni <= 0.0) return 0.0;
        if (ni <= c) return gpsshare(w, n, r);
        if (!urgent(n, r)) return 0.0;
        return gpsshare(w, prio_group(n, r), r);
    }
    /**
     * The population the LOAD-DEPENDENT factor is read at.
     *
     * PSPRIO uses the FULL vector in both branches while DPSPRIO and GPSPRIO
     * use the priority-restricted one above capacity. That asymmetry is
     * inherited from `State.afterEventStation` and is reproduced rather than
     * tidied: the two conventions give different lld factors on the same state,
     * and the reference's numbers were computed with these.
     */
    double prio_pop(const std::vector<double>& n, std::size_t r, double c) const {
        double ni = 0.0;
        for (double v : n) ni += v;
        if (ni <= c || !urgent(n, r)) return ni;
        double niprio = 0.0;
        prio_group(n, r, &niprio);
        return niprio;
    }
    /** `prioVec`: the vector DPSPRIO / GPSPRIO read the class dependence at. */
    std::vector<double> prio_vec(const std::vector<double>& n, std::size_t r, double c) const {
        double ni = 0.0;
        for (double v : n) ni += v;
        if (ni <= c || !urgent(n, r)) return n;
        return prio_group(n, r);
    }
    /**
     * `cdfac`: the class-r entry of eta_i(n) .* beta_{i,r}(n), evaluated on the
     * station's per-class population at firing time, as `State.afterEventInit`
     * folds the joint handle into the class one. 1 when neither is declared.
     */
    /** Utilization at a dependent station is T*S/peak, so a peak must be declared. */
    void check_peaks() const {
        for (std::size_t ist = 0; ist < M_; ++ist) {
            const auto& st = sn_.stations[ist];
            const bool on[2] = {static_cast<bool>(st.cdscaling), static_cast<bool>(st.jdscaling)};
            const std::vector<T>* pks[2] = {&st.cdscalingpeak, &st.jdscalingpeak};
            const char* names[2] = {"setClassDependence", "setJointDependence"};
            for (std::size_t h = 0; h < 2; ++h) {
                if (!on[h]) continue;
                bool ok = pks[h]->size() >= K_;
                for (std::size_t k = 0; ok && k < K_; ++k)
                    ok = num_traits<T>::to_double((*pks[h])[k]) > 0.0;
                if (!ok)
                    throw InputError("SolverSSA(method='nrm'): station '" + st.name +
                                     "' declares a dependent scaling with no declared peak rate. "
                                     "Utilization there is T*E[S]/peak, so pass the peak to " +
                                     names[h]);
            }
        }
    }
 
    double cdfac(std::size_t ist, const std::vector<double>& n, std::size_t r) const {
        if (ist == npos) return 1.0;
        const auto& st = sn_.stations[ist];
        if (!st.cdscaling && !st.jdscaling) return 1.0;
        std::vector<T> nt(n.size());
        for (std::size_t s = 0; s < n.size(); ++s) nt[s] = T(n[s]);
        double f = 1.0;
        for (const auto* h : {&st.cdscaling, &st.jdscaling}) {
            if (!*h) continue;
            const std::vector<T> v = (*h)(nt);
            f *= num_traits<T>::to_double(v[std::min(r, v.size() - 1)]);
        }
        return f;
    }
 
    // ---- the sample path ---------------------------------------------------
    double propensity(std::size_t j, const std::vector<double>& X,
                      const std::vector<std::vector<std::size_t>>& bufs,
                      const std::vector<Matrix<double>>& svc) const;
    std::size_t pick_from_buffer(const std::vector<std::size_t>& buf, std::size_t ist);
    std::size_t draw_entry_phase(std::size_t ind, std::size_t r);
    bool capacity_loss(const std::vector<double>& X, std::size_t slot) const;
    /** The PAS parameters of 0-based station `ist`; refuses a station with no mu(c). */
    const typename qn::NetworkStruct<T>::PasParam& pas_of(std::size_t ist) const;
    /** `oiIncrements`: Delta_mu of each position of the ordered list `c` (0-based classes). */
    std::vector<double> pas_increments(std::size_t ist, const std::vector<std::size_t>& c) const;
    /** `oirate`: the class-r departure rate of the list `c`. */
    double pas_rate(std::size_t ist, const std::vector<std::size_t>& c, std::size_t r) const;
    /** `pasInSvc`: class-r jobs at a served position (Delta_mu > 0). */
    double pas_in_service(std::size_t ist, const std::vector<std::size_t>& c, std::size_t r) const;
    /** `oiStarted`: tag a START for each position served in `cnew` and not in `cold`. */
    void pas_tag_started(std::size_t ind, const std::vector<std::size_t>& cold,
                         const std::vector<std::size_t>& cnew);
    /** `oiDepart`: the pass-and-swap rewrite of a class-r departure. */
    void pas_depart(std::size_t ind, std::vector<std::size_t>& c, std::size_t r);
    std::size_t fcr_refusing_region(const std::vector<double>& X, std::size_t src_node,
                                    std::size_t src_class, std::size_t dst_node,
                                    std::size_t dst_class) const;
    void fcr_outcome(const std::vector<double>& X, const Rx& x, std::size_t dst_node,
                     std::size_t dst_class, std::vector<std::vector<std::size_t>>& fcr_buf,
                     bool& lost, bool& blocked) const;
    std::size_t fcr_release_cascade(std::vector<double>& X,
                                    std::vector<std::vector<std::size_t>>& bufs,
                                    std::vector<std::vector<std::size_t>>& fcr_buf,
                                    std::vector<Matrix<double>>& svc, bool& svc_changed);
    bool capacity_block(const std::vector<double>& X, std::size_t slot,
                        std::size_t src_slot) const;
    std::size_t pick_preempted(const std::vector<double>& X,
                               const std::vector<std::size_t>& buf, std::size_t jnd,
                               std::size_t arr_class);
    void apply_arrival_buffer(std::size_t jnd, std::size_t s, const std::vector<double>& X,
                              std::vector<std::vector<std::size_t>>& bufs,
                              std::vector<Matrix<double>>& svc, bool& svc_changed);
    void update_buffers(std::size_t kfire, const std::vector<double>& X,
                        std::vector<std::vector<std::size_t>>& bufs, std::size_t dest_pos,
                        bool have_dest, std::vector<Matrix<double>>& svc, bool& svc_changed);
 
    /** Fill `sd_*` on a reaction leaving a node that routes by strategy. */
    void build_state_dependent_dest(Rx& x);
    /** Pick the destination slot of a state-dependent firing from the live state. */
    std::size_t resolve_state_dependent_dest(std::size_t kfire, const std::vector<double>& X);
 
    /**
     * Round-robin cursor per (SOURCE NODE, CLASS), keyed `node * K_ + cls`.
     * Not per reaction: a phase-type service splits one dispatcher across
     * several reactions and a pointer each would let the phases cycle
     * independently, which is not a round robin.
     */
    std::vector<std::size_t> rr_cursor_;
 
    /** Raise a START (or, with `preempt`, a PREEMPT) for class R at node IND. */
    void tag(std::size_t ind, std::size_t r, bool preempt) {
        if (!is_station_[ind] || r >= K_) return;
        const std::size_t ist = to_station_[ind];
        if (ist >= M_) return;
        // A Source CREATES jobs rather than admitting them to service, so it
        // seizes nothing; every codebase reports a zero row for it.
        if (sched_[ist] == SchedStrategy::EXT) return;
        (preempt ? preempt_cnt_ : start_cnt_)(ist, r) += 1.0;
    }
 
    static constexpr std::size_t npos = static_cast<std::size_t>(-1);
};
 
// ---------------------------------------------------------------------------
// Construction
// ---------------------------------------------------------------------------
 
/**
 * The phase slot map, `solver_ssa_nrm.m` lines 33-111.
 *
 * The state vector counts jobs per (node, class, PHASE) so phase-type service
 * is represented exactly instead of collapsed onto its mean rate. With one
 * phase everywhere the layout degenerates to the flat (node, class) index, so
 * an exponential model is unaffected -- which is the self-check the reference
 * names for this generalization.
 */
template <class T>
void NrmEngine<T>::build_layout() {
    I_ = sn_.nof_nodes();
    K_ = sn_.nclasses;
    M_ = sn_.nstations;
 
    is_station_.assign(I_, false);
    to_station_.assign(I_, npos);
    for (std::size_t i = 0; i < I_; ++i)
        if (sn_.nodes[i].station != 0) {
            is_station_[i] = true;
            to_station_[i] = sn_.nodes[i].station - 1;
        }
 
    sched_.assign(M_, SchedStrategy::FCFS);
    nservers_.assign(M_, 1.0);
    for (std::size_t i = 0; i < M_; ++i) {
        sched_[i] = sn_.stations[i].sched;
        nservers_[i] = sn_.stations[i].nservers;
    }
 
    // Phase counts. `sn.phasessz` floors at one, so a disabled pair still owns
    // exactly one slot and its zero rate silences the reaction.
    nph_.assign(I_, std::vector<std::size_t>(K_, 1));
    pie_.assign(I_, std::vector<std::vector<double>>(K_));
    for (std::size_t i = 0; i < I_; ++i) {
        for (std::size_t r = 0; r < K_; ++r) {
            pie_[i][r].assign(1, 1.0);
            if (!is_station_[i]) continue;
            const std::size_t ist = to_station_[i];
            if (sn_.disabled[ist][r]) continue;
            const mam::Map<T> m = lang::dist_to_map(sn_.service[ist][r]);
            const std::size_t n = m.D0.rows();
            if (n <= 1) continue;
            nph_[i][r] = n;
            const std::vector<T> p = lang::dist_pie(sn_.service[ist][r]);
            std::vector<double> pd(n, 0.0);
            double tot = 0.0;
            for (std::size_t k = 0; k < n && k < p.size(); ++k) {
                pd[k] = num_traits<T>::to_double(p[k]);
                if (!(pd[k] > 0.0)) pd[k] = 0.0;
                tot += pd[k];
            }
            if (tot > 0.0) {
                for (double& x : pd) x /= tot;
            } else {
                pd.assign(n, 0.0);
                pd[0] = 1.0;
            }
            pie_[i][r] = pd;
        }
    }
 
    phoff_.assign(I_, std::vector<std::size_t>(K_, 0));
    NS_ = 0;
    maxnph_ = 1;
    for (std::size_t i = 0; i < I_; ++i)
        for (std::size_t r = 0; r < K_; ++r) {
            phoff_[i][r] = NS_;
            NS_ += nph_[i][r];
            maxnph_ = std::max(maxnph_, nph_[i][r]);
        }
    slot_node_.assign(NS_, 0);
    slot_class_.assign(NS_, 0);
    for (std::size_t i = 0; i < I_; ++i)
        for (std::size_t r = 0; r < K_; ++r)
            for (std::size_t k = 0; k < nph_[i][r]; ++k) {
                slot_node_[phoff_[i][r] + k] = i;
                slot_class_[phoff_[i][r] + k] = r;
            }
 
    // Buffered phase-type service: only the jobs ACTUALLY in service carry a
    // phase, so their composition lives in svcph and not in the population.
    buf_ph_class_.assign(I_, std::vector<bool>(K_, false));
    buf_ph_node_.assign(I_, false);
    for (std::size_t i = 0; i < I_; ++i) {
        if (!is_station_[i] || !detail::sched_is_buf_ph(sched_[to_station_[i]])) continue;
        for (std::size_t r = 0; r < K_; ++r)
            if (nph_[i][r] > 1) {
                buf_ph_class_[i][r] = true;
                buf_ph_node_[i] = true;
            }
    }
 
    // Rates, server counts and the scaling tables the rate laws read.
    const double inf = std::numeric_limits<double>::infinity();
    mi_.assign(I_, inf);
    rate_.assign(I_, std::vector<double>(K_, 0.0));
    for (std::size_t i = 0; i < I_; ++i) {
        if (is_station_[i]) {
            const std::size_t ist = to_station_[i];
            mi_[i] = sn_.stations[ist].nservers;
            for (std::size_t r = 0; r < K_; ++r)
                if (!sn_.disabled[ist][r])
                    rate_[i][r] = num_traits<T>::to_double(sn_.rates(ist, r));
        } else {
            for (std::size_t r = 0; r < K_; ++r)
                rate_[i][r] = lang::GlobalConstants::Immediate;
        }
    }
 
    lld_.assign(M_, std::vector<double>());
    wnorm_.assign(M_, std::vector<double>(K_, 0.0));
    for (std::size_t i = 0; i < M_; ++i) {
        for (const T& x : sn_.stations[i].lldscaling)
            lld_[i].push_back(num_traits<T>::to_double(x));
        if (!detail::sched_is_weighted(sched_[i])) continue;
        double tot = 0.0;
        for (std::size_t r = 0; r < K_ && r < sn_.stations[i].schedparam.size(); ++r)
            tot += num_traits<T>::to_double(sn_.stations[i].schedparam[r]);
        if (!(tot > 0.0))
            throw InputError("solver_ssa_nrm: station '" + sn_.stations[i].name + "' has " +
                             std::string(lang::sched_to_text(sched_[i])) +
                             " scheduling with non-positive total weight");
        if (sn_.stations[i].nservers > 1.0)
            throw UnsupportedError(
                "solver_ssa_nrm: multi-server " + std::string(lang::sched_to_text(sched_[i])) +
                " at station '" + sn_.stations[i].name +
                "' is not supported; the reference's State.afterEventStation rejects it too");
        for (std::size_t r = 0; r < K_ && r < sn_.stations[i].schedparam.size(); ++r)
            wnorm_[i][r] = num_traits<T>::to_double(sn_.stations[i].schedparam[r]) / tot;
    }
 
    classprio_.assign(K_, 0.0);
    for (std::size_t r = 0; r < K_; ++r) classprio_[r] = sn_.classes[r].prio;
}
 
/**
 * The stoichiometry and the reaction list, `solver_ssa_nrm.m` lines 190-355.
 *
 * A departure is the ABSORPTION of the phase-type service process, so it fires
 * at the D1 row sum of the phase it leaves and the job re-enters its
 * destination in an entry phase drawn from pie -- which is why the destination
 * weight is the product of the routing probability and the entry probability. A
 * phase change is an off-diagonal of D0 and never leaves the node.
 */
/**
 * `fcrPrecompute`: the finite capacity regions, field for field.
 *
 * A region constrains an aggregate of the per-class populations of its member
 * stations, which is a linear function of the NRM state vector, so admission is
 * a 0/1 gate on the drawn destination. DROP loses the refused job at the gate,
 * as a balked arrival is lost; WAITQ parks it in the refusing region's FIFO,
 * which is state the reaction network does not carry and therefore rides beside
 * the buffers (`fcr_buf` in `run`).
 */
template <class T>
void NrmEngine<T>::build_regions() {
    const double inf = std::numeric_limits<double>::infinity();
    fcr_.clear();
    fcr_any_waitq_ = false;
    for (const typename qn::NetworkStruct<T>::Region& rg : sn_.regions) {
        FcrRegion f;
        f.member_node.assign(I_, 0);
        f.ccap.assign(K_, inf);
        f.sz.assign(K_, 1.0);
        f.waitq.assign(K_, 0);
        for (std::size_t i = 0; i < M_ && i < rg.members.size(); ++i) {
            if (!rg.members[i]) continue;
            f.member_node[sn_.station_to_node[i] - 1] = 1;
            for (std::size_t r = 0; r < K_; ++r)
                if (rg.cap[i][r] != -1.0) f.ccap[r] = std::min(f.ccap[r], rg.cap[i][r]);
            if (rg.cap[i][K_] != -1.0) f.gcap = std::min(f.gcap, rg.cap[i][K_]);
            if (i < rg.maxmem.size() && rg.maxmem[i] != -1.0)
                f.memcap = std::min(f.memcap, rg.maxmem[i]);
        }
        for (std::size_t r = 0; r < K_ && r < rg.size.size(); ++r)
            f.sz[r] = num_traits<T>::to_double(rg.size[r]);
        for (std::size_t r = 0; r < K_ && r < rg.rule.size(); ++r) {
            const qn::DropStrategy d = rg.rule[r];
            if (d == qn::DropStrategy::BAS || d == qn::DropStrategy::BBS ||
                d == qn::DropStrategy::RSRD)
                throw UnsupportedError(
                    "solver_ssa_nrm: finite capacity region '" + rg.name + "' declares a "
                    "blocking rule for class '" + sn_.classes[r].name + "'; only DROP and WAITQ "
                    "are region rules in any codebase's NRM");
            f.waitq[r] = d != qn::DropStrategy::DROP;
            if (f.waitq[r]) fcr_any_waitq_ = true;
        }
        if (rg.lincon_A.rows() > 0) {
            f.A = Matrix<double>(rg.lincon_A.rows(), rg.lincon_A.cols(), 0.0);
            for (std::size_t a = 0; a < rg.lincon_A.rows(); ++a)
                for (std::size_t c = 0; c < rg.lincon_A.cols(); ++c)
                    f.A(a, c) = num_traits<T>::to_double(rg.lincon_A(a, c));
            for (const T& v : rg.lincon_b) f.b.push_back(num_traits<T>::to_double(v));
        }
        fcr_.push_back(f);
    }
}
 
/**
 * `fcrRefusingRegion`: the 1-based index of the FIRST region that refuses a
 * class-`dst_class` job entering `dst_node` having just left `src_node` as
 * `src_class`, or 0 when every region admits it.
 *
 * Only a region holding the destination can refuse; a source in the same region
 * frees its slot first. `src_node == npos` is a WAITQ release, whose job already
 * left its source when it was parked, so no slot is freed.
 */
template <class T>
std::size_t NrmEngine<T>::fcr_refusing_region(const std::vector<double>& X, std::size_t src_node,
                                              std::size_t src_class, std::size_t dst_node,
                                              std::size_t dst_class) const {
    for (std::size_t f = 0; f < fcr_.size(); ++f) {
        const FcrRegion& g = fcr_[f];
        if (!g.member_node[dst_node]) continue;
        std::vector<double> x(K_, 0.0);
        for (std::size_t jnd = 0; jnd < I_; ++jnd)
            if (g.member_node[jnd])
                for (std::size_t r = 0; r < K_; ++r) x[r] += class_pop(X, jnd, r);
        if (src_node != npos && g.member_node[src_node]) x[src_class] -= 1.0;
        x[dst_class] += 1.0;
        double tot = 0.0, mem = 0.0;
        bool bad = false;
        for (std::size_t r = 0; r < K_; ++r) {
            tot += x[r];
            mem += x[r] * g.sz[r];
            if (x[r] > g.ccap[r]) bad = true;
        }
        if (tot > g.gcap || mem > g.memcap) bad = true;
        for (std::size_t a = 0; !bad && a < g.A.rows() && a < g.b.size(); ++a) {
            double lhs = 0.0;
            for (std::size_t c = 0; c < K_ && c < g.A.cols(); ++c) lhs += g.A(a, c) * x[c];
            if (lhs > g.b[a]) bad = true;
        }
        if (bad) return f + 1;
    }
    return 0;
}
 
template <class T>
const typename qn::NetworkStruct<T>::PasParam& NrmEngine<T>::pas_of(std::size_t ist) const {
    const auto it = sn_.pasparam.find(ist + 1);
    if (it == sn_.pasparam.end() || !it->second.svc_rate_fun)
        throw InputError("solver_ssa_nrm: PAS/OI station '" + sn_.stations[ist].name +
                         "' has no service rate function mu(c); set one with set_pas");
    return it->second;
}
 
template <class T>
std::vector<double> NrmEngine<T>::pas_increments(std::size_t ist,
                                                 const std::vector<std::size_t>& c) const {
    const auto& mu = pas_of(ist).svc_rate_fun;
    std::vector<double> inc(c.size(), 0.0);
    std::vector<std::size_t> prefix;
    double prev = 0.0;
    for (std::size_t p = 0; p < c.size(); ++p) {
        prefix.push_back(c[p] + 1);  // mu(c) reads 1-based class indices
        const double cur = num_traits<T>::to_double(mu(prefix));
        inc[p] = cur - prev;
        prev = cur;
    }
    return inc;
}
 
template <class T>
double NrmEngine<T>::pas_rate(std::size_t ist, const std::vector<std::size_t>& c,
                              std::size_t r) const {
    if (c.empty()) return 0.0;
    const std::vector<double> inc = pas_increments(ist, c);
    std::vector<std::size_t> c1(c.size());
    for (std::size_t p = 0; p < c.size(); ++p) c1[p] = c[p] + 1;
    const auto& G = pas_of(ist).swap_graph;
    double rt = 0.0;
    for (std::size_t p = 0; p < c.size(); ++p) {
        if (inc[p] <= 0.0) continue;  // position p receives no service
        if (qn::pass_and_swap<T>(c1, p, G).second == r + 1) rt += inc[p];
    }
    return rt;
}
 
template <class T>
double NrmEngine<T>::pas_in_service(std::size_t ist, const std::vector<std::size_t>& c,
                                    std::size_t r) const {
    if (c.empty()) return 0.0;
    const std::vector<double> inc = pas_increments(ist, c);
    double n = 0.0;
    for (std::size_t p = 0; p < c.size(); ++p)
        if (inc[p] > 0.0 && c[p] == r) n += 1.0;
    return n;
}
 
template <class T>
void NrmEngine<T>::pas_tag_started(std::size_t ind, const std::vector<std::size_t>& cold,
                                   const std::vector<std::size_t>& cnew) {
    const std::size_t ist = to_station_[ind];
    const std::vector<double> inc_new = pas_increments(ist, cnew);
    const std::vector<double> inc_old = pas_increments(ist, cold);
    for (std::size_t p = 0; p < inc_new.size(); ++p) {
        if (inc_new[p] <= 0.0) continue;
        if (p < inc_old.size() && inc_old[p] > 0.0) continue;  // already served
        tag(ind, cnew[p], false);
    }
}
 
/**
 * The completing position is drawn among those whose pass-and-swap ejects class
 * r, weighted by that position's own Delta_mu -- the split whose total is the
 * propensity that just fired.
 */
template <class T>
void NrmEngine<T>::pas_depart(std::size_t ind, std::vector<std::size_t>& c, std::size_t r) {
    const std::size_t ist = to_station_[ind];
    const std::vector<double> inc = pas_increments(ist, c);
    std::vector<std::size_t> c1(c.size());
    for (std::size_t p = 0; p < c.size(); ++p) c1[p] = c[p] + 1;
    const auto& G = pas_of(ist).swap_graph;
    std::vector<std::size_t> pos;
    std::vector<double> w;
    double tot = 0.0;
    for (std::size_t p = 0; p < c.size(); ++p) {
        if (inc[p] <= 0.0) continue;
        if (qn::pass_and_swap<T>(c1, p, G).second != r + 1) continue;
        pos.push_back(p);
        w.push_back(inc[p]);
        tot += inc[p];
    }
    if (pos.empty()) return;  // this class cannot depart from the current list
    const double u = rng_.uniform() * tot;
    double acc = 0.0;
    std::size_t pick = pos.back();
    for (std::size_t x = 0; x < pos.size(); ++x) {
        acc += w[x];
        if (u < acc) {
            pick = pos[x];
            break;
        }
    }
    const std::vector<std::size_t> nc = qn::pass_and_swap<T>(c1, pick, G).first;
    c.assign(nc.size(), 0);
    for (std::size_t p = 0; p < nc.size(); ++p) c[p] = nc[p] - 1;
}
 
/**
 * `pollLandCtrl`: SERVING q with the visit budget, SWITCHING into q with the
 * entry phase drawn from the switchover PH, or PARKED. A landing in SERVING is
 * the moment a class-q job takes the server, so it raises that START here: the
 * reference raises none at a polling station at all, which leaves StartN zero
 * there although jobs plainly start service.
 */
template <class T>
typename NrmEngine<T>::PollCtrl NrmEngine<T>::poll_land(std::size_t ind, std::size_t q, int mode,
                                                        long budget, bool raise) {
    PollCtrl c;
    c.pos = q;
    c.mode = mode;
    if (mode == 1) {
        c.ctr = budget;
        if (raise) tag(ind, q - 1, false);
    } else if (mode == 2) {
        const std::vector<T>& pie = poll_info_[ind].sw_pie[q - 1];
        std::vector<double> w(pie.size());
        for (std::size_t k = 0; k < pie.size(); ++k) w[k] = num_traits<T>::to_double(pie[k]);
        c.swk = rng_.draw(w);
    }
    return c;
}
 
/**
 * What a region refusal does to the firing, by the refusing region's rule for
 * the destination class.
 *
 * WAITQ parks the job in that region's FIFO: it leaves its source and joins no
 * station until `fcr_release_cascade` admits it. DROP LOSES an open job, as a
 * balked arrival is lost. A CLOSED job refused under DROP is BLOCKED instead --
 * the firing is cancelled and the job keeps its place -- which is this
 * codebase's reading of the rule (SolverCTMC censors the transition and the
 * serial engine does not take it, so the population is conserved). The
 * reference's NRM, CTMC and serial engines instead all lose the closed job, so
 * the population decays until the region can never refuse again: 1 job of 3 on
 * the capped three-station cycle. That cross-language divergence belongs to the
 * rule, not to this engine; see `_kb/07-cross-language-parity.md`.
 */
template <class T>
void NrmEngine<T>::fcr_outcome(const std::vector<double>& X, const Rx& x, std::size_t dst_node,
                               std::size_t dst_class,
                               std::vector<std::vector<std::size_t>>& fcr_buf, bool& lost,
                               bool& blocked) const {
    const std::size_t fref = fcr_refusing_region(X, x.node, x.cls, dst_node, dst_class);
    if (fref == 0) return;
    if (fcr_[fref - 1].waitq[dst_class]) {
        lost = true;
        fcr_buf[fref - 1].push_back(dst_node * K_ + dst_class);
    } else if (std::isinf(sn_.classes[dst_class].population)) {
        lost = true;
    } else {
        blocked = true;
    }
}
 
/**
 * `fcrReleaseCascade`: strict-FIFO head-of-line release of the parked WAITQ
 * jobs. Each region's head is admitted while its constraints permit, joining the
 * destination as a routed arrival does (entry phase drawn at release, buffer
 * join, START tag); the loop repeats until a full pass frees nothing, so a
 * release that relieves another region cascades. Returns the jobs released.
 */
template <class T>
std::size_t NrmEngine<T>::fcr_release_cascade(std::vector<double>& X,
                                              std::vector<std::vector<std::size_t>>& bufs,
                                              std::vector<std::vector<std::size_t>>& fcr_buf,
                                              std::vector<Matrix<double>>& svc,
                                              bool& svc_changed) {
    std::size_t released = 0;
    bool progress = true;
    while (progress) {
        progress = false;
        for (std::size_t f = 0; f < fcr_buf.size(); ++f) {
            if (fcr_buf[f].empty()) continue;
            const std::size_t tok = fcr_buf[f].front();
            const std::size_t jnd = tok / K_, r = tok % K_;
            if (fcr_refusing_region(X, npos, r, jnd, r) != 0) continue;  // head-of-line
            // A buffered-PH destination counts the job in its class total slot;
            // whether and in which phase it enters service is decided against the
            // server occupancy by `apply_arrival_buffer`, as for a routed arrival.
            const std::size_t ke = buf_ph_node_[jnd] ? 0 : draw_entry_phase(jnd, r);
            X[phoff_[jnd][r] + ke] += 1.0;
            apply_arrival_buffer(jnd, r, X, bufs, svc, svc_changed);
            fcr_buf[f].erase(fcr_buf[f].begin());
            ++released;
            progress = true;
        }
    }
    return released;
}
 
template <class T>
void NrmEngine<T>::build_reactions() {
    rx_.clear();
    // A class READS at a cache iff it has a popularity row and either a hit class
    // or is one of the cache's retrieval classes, which read to COMPLETE a miss.
    is_cache_node_.assign(I_, 0);
    cache_read_.assign(I_, std::vector<char>(K_, 0));
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (sn_.nodes[ind].nodetype != qn::NodeType::Cache) continue;
        const auto ci = sn_.nodeparam.find(ind + 1);
        if (ci == sn_.nodeparam.end())
            throw InputError("solver_ssa_nrm: Cache node '" + sn_.nodes[ind].name +
                             "' has no cache parameters");
        const qn::CacheParam<T>& cp = ci->second;
        is_cache_node_[ind] = 1;
        std::vector<std::size_t> rcl, rci, rco;
        if (cp.retrieval_capacity > 0) qn::cache_retrieval_class_map(cp, rcl, rci, rco);
        for (std::size_t r = 0; r < K_; ++r) {
            if (r >= cp.pread.size() || cp.pread[r].empty()) continue;
            const bool is_retr = std::find(rcl.begin(), rcl.end(), r + 1) != rcl.end();
            const std::size_t hc = r < cp.hitclass.size() ? cp.hitclass[r] : 0;
            const std::size_t mc = r < cp.missclass.size() ? cp.missclass[r] : 0;
            if (hc == 0 && !is_retr) continue;
            // A read class that is its own hit or miss class would read again on
            // the spot: the reference NRM livelocks there, so it is refused.
            if (hc == r + 1 || mc == r + 1)
                throw UnsupportedError(
                    "solver_ssa_nrm: class '" + sn_.classes[r].name + "' at Cache '" +
                    sn_.nodes[ind].name + "' is its own hit or miss class, so the NRM's cache "
                    "read would fire again at once; ask for method='serial'");
            cache_read_[ind][r] = 1;
        }
    }
    // Departures, one per (node, class, phase).
    for (std::size_t ind = 0; ind < I_; ++ind) {
        for (std::size_t r = 0; r < K_; ++r) {
            for (std::size_t kk = 0; kk < nph_[ind][r]; ++kk) {
                Rx x;
                x.node = ind;
                x.cls = r;
                x.dep_phase = kk;
                x.is_buf_svc = buf_ph_class_[ind][r];
                // At a buffered-PH source the population lives entirely in the
                // first phase slot; which phase completed is carried in
                // dep_phase and consumed from svcph at firing.
                x.from = buf_ph_class_[ind][r] ? phoff_[ind][r] : phoff_[ind][r] + kk;
                x.rate = rate_[ind][r];
                x.is_cache = is_cache_node_[ind] && cache_read_[ind][r];
                if (is_station_[ind] && nph_[ind][r] > 1) {
                    const std::size_t ist = to_station_[ind];
                    const mam::Map<T> m = lang::dist_to_map(sn_.service[ist][r]);
                    double s = 0.0;
                    for (std::size_t c = 0; c < m.D1.cols(); ++c)
                        s += num_traits<T>::to_double(m.D1(kk, c));
                    x.rate = s;
                }
                rx_.push_back(x);
            }
        }
    }
    const std::size_t ndep = rx_.size();
 
    // Phase transitions, one per (node, class, ka -> kb) with D0(ka,kb) > 0.
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (!is_station_[ind]) continue;
        const std::size_t ist = to_station_[ind];
        for (std::size_t r = 0; r < K_; ++r) {
            if (nph_[ind][r] <= 1 || sn_.disabled[ist][r]) continue;
            const mam::Map<T> m = lang::dist_to_map(sn_.service[ist][r]);
            for (std::size_t ka = 0; ka < nph_[ind][r]; ++ka)
                for (std::size_t kb = 0; kb < nph_[ind][r]; ++kb) {
                    if (ka == kb) continue;
                    const double d = num_traits<T>::to_double(m.D0(ka, kb));
                    if (!(d > 0.0)) continue;
                    Rx x;
                    x.node = ind;
                    x.cls = r;
                    x.is_phase = true;
                    x.phase_from = ka;
                    x.phase_to = kb;
                    x.rate = d;
                    // A buffered-PH phase change moves a job inside svcph and
                    // leaves the population alone, so its stoichiometry column
                    // is all zeros and its propensity reads the first slot.
                    x.from = buf_ph_class_[ind][r] ? phoff_[ind][r] : phoff_[ind][r] + ka;
                    rx_.push_back(x);
                }
        }
    }
 
    // Polling switchover reactions, `solver_ssa_nrm.m` lines 419-482. A walk
    // between buffers is a timed event that moves no job, so it is one reaction
    // per polling node with a timed switchover: an all-zero stoichiometry column
    // whose propensity reads the controller. Only EXPONENTIAL service is taken at
    // a polling station: the controller tracks no in-service phase, and spreading
    // the class over phases would serve several fictitious phase-jobs at once.
    is_poll_.assign(I_, 0);
    poll_info_.assign(I_, qn::PollingInfo<T>());
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (!is_station_[ind] || sched_[to_station_[ind]] != SchedStrategy::POLLING) continue;
        is_poll_[ind] = 1;
        poll_info_[ind] = qn::polling_info(sn_, ind + 1);
        for (std::size_t r = 0; r < K_; ++r)
            if (nph_[ind][r] > 1)
                throw UnsupportedError(
                    "solver_ssa_nrm: polling station '" + sn_.stations[to_station_[ind]].name +
                    "' serves class '" + sn_.classes[r].name + "' with phase-type service. The "
                    "NRM polling controller supports exponential service only, as the reference's "
                    "does; ask for method='serial'");
        bool anysw = false;
        for (bool h : poll_info_[ind].has_sw) anysw = anysw || h;
        if (!anysw) continue;  // every switchover immediate: the server never dwells
        Rx x;
        x.node = ind;
        x.cls = 0;
        x.from = phoff_[ind][0];
        x.is_sw = true;
        rx_.push_back(x);
    }
 
    // The stoichiometry matrix, states x reactions.
    S_ = Matrix<double>(NS_, rx_.size(), 0.0);
    for (std::size_t k = 0; k < rx_.size(); ++k) {
        Rx& x = rx_[k];
        if (x.is_sw) continue;  // a switchover moves no job
        if (x.is_cache) {
            // The read consumes the job; where it reappears is drawn at firing.
            S_(x.from, k) -= 1.0;
            continue;
        }
        if (x.is_phase) {
            if (!buf_ph_class_[x.node][x.cls]) {
                S_(phoff_[x.node][x.cls] + x.phase_from, k) -= 1.0;
                S_(phoff_[x.node][x.cls] + x.phase_to, k) += 1.0;
            }
            continue;
        }
        S_(x.from, k) -= 1.0;
        for (std::size_t jnd = 0; jnd < I_; ++jnd)
            for (std::size_t s = 0; s < K_; ++s) {
                const double p =
                    num_traits<T>::to_double(sn_.rtnodes(x.node * K_ + x.cls, jnd * K_ + s));
                if (!(p > 0.0)) continue;
                if (buf_ph_class_[jnd][s]) {
                    // The arrival lands in the total-population slot; whether it
                    // enters service, and in which phase, is decided at firing
                    // from the server occupancy and pie.
                    S_(phoff_[jnd][s], k) += p;
                } else {
                    for (std::size_t ke = 0; ke < nph_[jnd][s]; ++ke) {
                        const double pe = pie_[jnd][s][ke];
                        if (!(pe > 0.0)) continue;
                        S_(phoff_[jnd][s] + ke, k) += p * pe;
                    }
                }
            }
    }
 
    // `P = S; P(P<0) = P(P<0)+1` selects the destination draw: nnz(P) > 1 means
    // the firing must pick among several destination slots, and the weights are
    // read off P so a class that routes back into its own slot competes with the
    // others on equal footing.
    for (std::size_t k = 0; k < rx_.size(); ++k) {
        Rx& x = rx_[k];
        std::vector<double> Pcol(NS_, 0.0);
        for (std::size_t i = 0; i < NS_; ++i) {
            const double v = S_(i, k);
            Pcol[i] = v < 0.0 ? v + 1.0 : v;
            if (v < 0.0) x.from_slots.push_back(i);
            if (v > 0.0 && x.det_dest == npos) x.det_dest = i;
        }
        for (std::size_t i = 0; i < NS_; ++i)
            if (Pcol[i] != 0.0) ++x.nnzP;
        if (x.nnzP > 1) {
            double acc = 0.0;
            for (std::size_t i = 0; i < NS_; ++i)
                if (Pcol[i] != 0.0) {
                    x.to_slots.push_back(i);
                    acc += Pcol[i];
                    x.cdf.push_back(acc);
                }
        }
        if (!x.is_phase && !x.is_sw) build_state_dependent_dest(x);
    }
 
    // Departure reactions of each (station, class): the reference rescans all
    // reactions inside the metric loop, which is the same set every step.
    dep_rx_.assign(M_, std::vector<std::vector<std::size_t>>(K_));
    for (std::size_t k = 0; k < ndep; ++k)
        if (is_station_[rx_[k].node]) dep_rx_[to_station_[rx_[k].node]][rx_[k].cls].push_back(k);
 
    // Stations and classes where a refusal BLOCKS; see `capacity_block`.
    // A cap that CANNOT BIND is not a blocking site. Every closed model carries
    // cap[ist] = N by default, and a station that can hold the whole population
    // never refuses one: the arriving job is itself one of the N, so the
    // pre-arrival count is at most N-1. Excluding those keeps `blk_on_` false --
    // and the per-firing test unpaid -- on ordinary models. An open class present
    // anywhere makes a finite station cap binding again, its jobs not being in N.
    double closed_total = 0.0;
    bool any_open = false;
    for (std::size_t r = 0; r < K_; ++r) {
        const double nj = sn_.classes[r].population;
        if (std::isinf(nj))
            any_open = true;
        else
            closed_total += nj;
    }
    blk_can_.assign(M_, std::vector<char>(K_, 0));
    blk_on_ = false;
    for (std::size_t ist = 0; ist < M_; ++ist)
        for (std::size_t r = 0; r < K_; ++r) {
            const double nj = sn_.classes[r].population;
            if (std::isinf(nj)) continue;  // open: refusal LOSES
            const double ccap = sn_.classcap[ist][r];
            const bool binds_st =
                !std::isinf(sn_.cap[ist]) && (any_open || sn_.cap[ist] < closed_total);
            const bool binds_cl = ccap > 0.0 && !std::isinf(ccap) && ccap < nj;
            if (binds_st || binds_cl) {
                blk_can_[ist][r] = 1;
                blk_on_ = true;
            }
        }
}
 
/**
 * The candidate destinations of a reaction whose node routes BY STRATEGY.
 *
 * The candidates are the DECLARED out-arcs, `P(r,s)(i,j) > 0`, enumerated in
 * the order `refresh_routing` uses, and not the probability split it wrote into
 * `rtnodes` from them: the split is what a matrix solver needs and it is not
 * what the dispatcher does. Nothing is filled for PROB, RAND or DISABLED, which
 * are draws or refusals and stay on the `cdf` path.
 */
template <class T>
void NrmEngine<T>::build_state_dependent_dest(Rx& x) {
    const qn::NodeDef& nd = sn_.nodes[x.node];
    const lang::RoutingStrategy rs = nd.routing.size() > x.cls
                                         ? nd.routing[x.cls]
                                         : lang::RoutingStrategy::PROB;
    if (rs != lang::RoutingStrategy::RROBIN && rs != lang::RoutingStrategy::WRROBIN &&
        rs != lang::RoutingStrategy::JSQ && rs != lang::RoutingStrategy::SQ)
        return;
    for (std::size_t s = 1; s <= K_; ++s)
        for (std::size_t j = 1; j <= I_; ++j) {
            if (!(num_traits<T>::to_double(sn_.get_route(x.cls + 1, s, x.node + 1, j)) > 0.0))
                continue;
            x.sd_node.push_back(j - 1);
            x.sd_class.push_back(s - 1);
        }
    if (x.sd_node.size() < 2) {  // one arc is not a dispatch decision
        x.sd_node.clear();
        x.sd_class.clear();
        return;
    }
    x.sd_strategy = rs;
    // The entry-phase law of each candidate, so the phase is drawn AFTER the
    // destination is chosen rather than folded into one joint distribution.
    x.sd_slots.resize(x.sd_node.size());
    x.sd_pie.resize(x.sd_node.size());
    for (std::size_t d = 0; d < x.sd_node.size(); ++d) {
        const std::size_t jnd = x.sd_node[d], s = x.sd_class[d];
        if (buf_ph_class_[jnd][s]) {
            x.sd_slots[d].push_back(phoff_[jnd][s]);
            x.sd_pie[d].push_back(1.0);
            continue;
        }
        for (std::size_t ke = 0; ke < nph_[jnd][s]; ++ke) {
            const double pe = pie_[jnd][s][ke];
            if (pe > 0.0) {
                x.sd_slots[d].push_back(phoff_[jnd][s] + ke);
                x.sd_pie[d].push_back(pe);
            }
        }
    }
    if (rs == lang::RoutingStrategy::WRROBIN) {
        const std::map<std::size_t, double>* w =
            nd.routing_weights.size() > x.cls ? &nd.routing_weights[x.cls] : NULL;
        double total = 0.0;
        x.sd_weight.assign(x.sd_node.size(), 0.0);
        if (w != NULL)
            for (std::size_t d = 0; d < x.sd_node.size(); ++d) {
                const std::map<std::size_t, double>::const_iterator it =
                    w->find(x.sd_node[d] + 1);
                x.sd_weight[d] = (it == w->end()) ? 0.0 : it->second;
                total += x.sd_weight[d];
            }
        if (!(total > 0.0)) {
            // No weights declared is plain round robin, which is what the
            // reference's WRROBIN degenerates to with a uniform weight vector.
            x.sd_weight.clear();
            x.sd_strategy = lang::RoutingStrategy::RROBIN;
        } else {
            for (std::size_t d = 0; d < x.sd_weight.size(); ++d) x.sd_weight[d] /= total;
        }
    }
}
 
/**
 * The destination slot of a state-dependent firing, from the LIVE population.
 *
 * Each rule is the reference's (`solver_ssa_nrm.m` lines 1162-1200, 1390-1430):
 *
 *   JSQ      the candidate holding the smallest TOTAL population, each
 *            evaluated on its own queue and never on the routing node's, with
 *            TIES SPLIT UNIFORMLY. The tie split is not a detail: on a fleet of
 *            identical servers the system starts with every queue equal, so a
 *            first-candidate tie-break sends a run of jobs to candidate one and
 *            reports a skew the policy does not have.
 *   SQ(d)    the power-of-d-choices rule, NOT "shortest queue": sample d
 *            candidates uniformly WITH replacement and join the smallest of
 *            those, ties by first occurrence in the sampled tuple. `d` defaults
 *            to 2, this port carrying no per-node override for it.
 *   RROBIN   a pointer advanced once per firing, kept per (SOURCE NODE, CLASS)
 *            and NOT per reaction: a phase-type service splits one dispatcher
 *            across several reactions, and a pointer each would let the phases
 *            cycle independently and undo the determinism.
 *   WRROBIN  the same pointer, spent in proportion to the declared weights.
 *
 * The destination NODE is what the rule fixes; the entry PHASE is drawn after,
 * from `pie` among that node's candidates, since taking the first match would
 * bias the service time.
 */
template <class T>
std::size_t NrmEngine<T>::resolve_state_dependent_dest(std::size_t kfire,
                                                       const std::vector<double>& X) {
    Rx& x = rx_[kfire];
    const std::size_t nd = x.sd_node.size();
    const std::size_t cursor_key = x.node * K_ + x.cls;
    std::size_t pick = 0;
    if (x.sd_strategy == lang::RoutingStrategy::JSQ) {
        std::vector<double> load(nd, 0.0);
        double best = 0.0;
        for (std::size_t d = 0; d < nd; ++d) {
            const std::vector<double> cc = class_counts(X, x.sd_node[d]);
            for (double v : cc) load[d] += v;
            if (d == 0 || load[d] < best) best = load[d];
        }
        std::vector<std::size_t> amins;
        for (std::size_t d = 0; d < nd; ++d)
            if (load[d] == best) amins.push_back(d);
        pick = amins[std::min(amins.size() - 1,
                              std::size_t(rng_.uniform() * double(amins.size())))];
    } else if (x.sd_strategy == lang::RoutingStrategy::SQ) {
        const std::size_t dsample = std::min<std::size_t>(2, nd);
        double best = 0.0;
        for (std::size_t t = 0; t < dsample; ++t) {
            const std::size_t cand =
                std::min(nd - 1, std::size_t(rng_.uniform() * double(nd)));
            const std::vector<double> cc = class_counts(X, x.sd_node[cand]);
            double load = 0.0;
            for (double v : cc) load += v;
            if (t == 0 || load < best) {
                best = load;
                pick = cand;
            }
        }
    } else if (x.sd_strategy == lang::RoutingStrategy::WRROBIN) {
        // The pointer is spent in proportion to the declared weights: candidate
        // d owns the slice [sum w_{<d}, sum w_{<=d}) of one full turn.
        const std::size_t c = rr_cursor_[cursor_key];
        const double u = double(c % nd) / double(nd);
        double acc = 0.0;
        pick = nd - 1;
        for (std::size_t d = 0; d < nd; ++d) {
            acc += x.sd_weight[d];
            if (u < acc - 1e-12) {
                pick = d;
                break;
            }
        }
        rr_cursor_[cursor_key] = c + 1;
    } else {
        pick = rr_cursor_[cursor_key] % nd;
        rr_cursor_[cursor_key] = rr_cursor_[cursor_key] + 1;
    }
    // the entry phase, drawn within the chosen destination
    const std::vector<std::size_t>& slots = x.sd_slots[pick];
    if (slots.size() == 1) return slots[0];
    const double u = rng_.uniform();
    double acc = 0.0;
    for (std::size_t k = 0; k + 1 < slots.size(); ++k) {
        acc += x.sd_pie[pick][k];
        if (u < acc) return slots[k];
    }
    return slots.back();
}
 
/**
 * The dependency sets, `solver_ssa_nrm.m` lines 1003-1050.
 *
 * A reaction's firing changes the slots its stoichiometry touches; every rate
 * law reads its node's WHOLE class-count vector (and, for the per-phase share,
 * the sibling phases of its own class), so the affected set is every slot of
 * every touched node, and the reactions to refresh are those that consume from
 * one of those slots.
 */
template <class T>
void NrmEngine<T>::build_dependencies() {
    const std::size_t nrx = rx_.size();
    D_.assign(nrx, std::vector<std::size_t>());
    for (std::size_t k = 0; k < nrx; ++k) {
        std::vector<bool> touched_node(I_, false);
        for (std::size_t i = 0; i < NS_; ++i)
            if (S_(i, k) != 0.0) touched_node[slot_node_[i]] = true;
        std::vector<bool> hit(nrx, false);
        for (std::size_t ind = 0; ind < I_; ++ind) {
            if (!touched_node[ind]) continue;
            for (std::size_t r = 0; r < K_; ++r)
                for (std::size_t p = 0; p < nph_[ind][r]; ++p) {
                    const std::size_t slot = phoff_[ind][r] + p;
                    for (std::size_t j = 0; j < nrx; ++j)
                        if (S_(slot, j) < 0.0) hit[j] = true;
                }
        }
        for (std::size_t j = 0; j < nrx; ++j)
            if (hit[j]) D_[k].push_back(j);
    }
}
 
/**
 * The initial state.
 *
 * THIS PORT HAS NO State PACKAGE, so the reference's `sn.state` cannot be read;
 * the same position `solver_fluid` takes applies here and for the same reason.
 * A closed class starts entirely at its reference station and an open class
 * puts the single fictitious token at the Source that the EXT rate law needs
 * (the reference substitutes exactly that for the infinite marginal). The chain
 * is ergodic, so the steady-state averages this engine reports do not depend on
 * the choice; a TRANSIENT would, which is why `getTranAvg` is not offered.
 */
template <class T>
void NrmEngine<T>::build_initial_state() {
    nvec0_.assign(NS_, 0.0);
    // The cache contents: the declared warm start when there is one (its row
    // after the per-class counts), else slot i holding item i, as the reference
    // NRM seeds it; with a retrieval system blocks A and B follow, empty.
    cache_var0_.assign(I_, std::vector<T>());
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (!is_cache_node_[ind]) continue;
        const qn::CacheParam<T>& cp = sn_.nodeparam.find(ind + 1)->second;
        std::size_t tcc = 0;
        for (int c : cp.itemcap)
            if (c > 0) tcc += static_cast<std::size_t>(c);
        std::vector<std::size_t> rcl, rci, rco;
        const bool retr = cp.retrieval_capacity > 0 && !cp.retrieval_classes.empty();
        if (retr) qn::cache_retrieval_class_map(cp, rcl, rci, rco);
        std::vector<T>& v = cache_var0_[ind];
        if (cp.initstate.size() > K_) {
            v.assign(cp.initstate.begin() + static_cast<std::ptrdiff_t>(K_), cp.initstate.end());
        } else {
            for (std::size_t c = 0; c < tcc; ++c)
                v.push_back(num_traits<T>::from_int(static_cast<long>(c + 1)));
        }
        const std::size_t want = tcc + (retr ? cp.nitems + rcl.size() : 0);
        v.resize(want, num_traits<T>::from_int(0));
    }
    buffers0_.assign(I_, std::vector<std::size_t>());
    svcph0_.assign(I_, Matrix<double>());
    for (std::size_t i = 0; i < I_; ++i)
        if (buf_ph_node_[i]) svcph0_[i] = Matrix<double>(K_, maxnph_, 0.0);
 
    std::vector<std::vector<double>> nir(I_, std::vector<double>(K_, 0.0));
    for (std::size_t r = 0; r < K_; ++r) {
        const double pop = sn_.classes[r].population;
        if (std::isinf(pop)) {
            if (sn_.sourceIdx == 0)
                throw InputError(
                    "solver_ssa_nrm: the open class '" + sn_.classes[r].name +
                    "' has no Source; an open model must carry one for the arrival reaction");
            nir[sn_.station_to_node[sn_.sourceIdx - 1] - 1][r] = 1.0;
        } else if (pop > 0.0) {
            const std::size_t rs = sn_.classes[r].refstat;
            if (rs < 1 || rs > M_)
                throw InputError("solver_ssa_nrm: class '" + sn_.classes[r].name +
                                 "' has no reference station");
            nir[sn_.station_to_node[rs - 1] - 1][r] = pop;
        }
    }
 
    for (std::size_t ind = 0; ind < I_; ++ind) {
        for (std::size_t r = 0; r < K_; ++r) {
            const double n = nir[ind][r];
            if (n <= 0.0) continue;
            if (nph_[ind][r] <= 1 || buf_ph_class_[ind][r]) {
                nvec0_[phoff_[ind][r]] = n;
            } else {
                double left = n;
                for (std::size_t ke = 0; ke < nph_[ind][r]; ++ke) {
                    const double take = (ke + 1 == nph_[ind][r])
                                            ? left
                                            : std::min(left, std::round(n * pie_[ind][r][ke]));
                    nvec0_[phoff_[ind][r] + ke] = take;
                    left -= take;
                }
            }
        }
        // A PAS / OI list holds EVERY job, oldest first, in class order.
        if (is_station_[ind] && detail::sched_is_list(sched_[to_station_[ind]]))
            for (std::size_t r = 0; r < K_; ++r)
                for (std::size_t c = 0; c < static_cast<std::size_t>(nir[ind][r]); ++c)
                    buffers0_[ind].push_back(r);
        // The buffer holds exactly the waiting jobs, `numel(buf) == max(0,
        // total - mi)`; which classes they are is immaterial to the steady
        // state, so they are taken in class order.
        if (is_station_[ind] && detail::sched_is_buffered(sched_[to_station_[ind]])) {
            double total = 0.0;
            for (std::size_t r = 0; r < K_; ++r) total += nir[ind][r];
            double waiting = std::max(0.0, total - mi_[ind]);
            for (std::size_t r = 0; r < K_ && waiting > 0.0; ++r) {
                double take = std::min(waiting, nir[ind][r]);
                for (std::size_t c = 0; c < static_cast<std::size_t>(take); ++c)
                    buffers0_[ind].push_back(r);
                waiting -= take;
            }
        }
        if (!buf_ph_node_[ind]) continue;
        for (std::size_t r = 0; r < K_; ++r) {
            double waiting_r = 0.0;
            for (std::size_t c : buffers0_[ind])
                if (c == r) waiting_r += 1.0;
            const double insvc = std::max(0.0, nir[ind][r] - waiting_r);
            if (nph_[ind][r] <= 1) {
                svcph0_[ind](r, 0) = insvc;
            } else {
                double left = insvc;
                for (std::size_t ke = 0; ke < nph_[ind][r]; ++ke) {
                    const double take = (ke + 1 == nph_[ind][r])
                                            ? left
                                            : std::min(left, std::round(insvc * pie_[ind][r][ke]));
                    svcph0_[ind](r, ke) = take;
                    left -= take;
                }
            }
        }
    }
}
 
// ---------------------------------------------------------------------------
// The rate laws
// ---------------------------------------------------------------------------
 
/**
 * The propensity of reaction `j`, `solver_ssa_nrm.m` lines 795-948.
 *
 * The reference builds one closure per reaction; this is the same switch
 * evaluated on demand, which avoids an indirect call per reaction per step and
 * cannot drift from the utilization accumulators that reuse the same sharing
 * factors.
 */
template <class T>
double NrmEngine<T>::propensity(std::size_t j, const std::vector<double>& X,
                                const std::vector<std::vector<std::size_t>>& bufs,
                                const std::vector<Matrix<double>>& svc) const {
    const Rx& x = rx_[j];
    const std::size_t ind = x.node, r = x.cls;
    if (!is_station_[ind])
        return x.rate * kir_frac(X, x.from, ind, r) * std::min(1.0, class_pop(X, ind, r));
 
    const std::size_t ist = to_station_[ind];
    const double eps = lang::GlobalConstants::Zero;
 
    // `pollSwRate`: the total leaving rate -D0(swk,swk) of the switchover phase
    // the walking server occupies; which transition it takes is drawn at firing.
    if (x.is_sw) {
        const PollCtrl& c = poll_ctrl_[ind];
        if (c.mode != 2) return 0.0;
        return -num_traits<T>::to_double(poll_info_[ind].sw_d0[c.pos - 1](c.swk, c.swk));
    }
 
    // Buffered phase-type service: the rate reads the in-service multiset, not
    // the class total, which also counts the jobs still waiting.
    if (buf_ph_class_[ind][r]) {
        const std::size_t kk = x.is_phase ? x.phase_from : x.dep_phase;
        const std::vector<double> n = class_counts(X, ind);
        double tot = 0.0;
        for (double v : n) tot += v;
        return x.rate * svc[ind](r, kk) * lldfac(ist, tot) * cdfac(ist, n, r);
    }
 
    const double kf = kir_frac(X, x.from, ind, r);
    switch (sched_[ist]) {
        case SchedStrategy::EXT:
            // A Source fires at a constant arrival rate: its token is
            // fictitious, so applying the population share would silence it and
            // deadlock every open model.
            return x.rate;
        case SchedStrategy::INF:
            return x.rate * kf * class_pop(X, ind, r) * cdfac(ist, class_counts(X, ind), r);
        case SchedStrategy::PS:
        case SchedStrategy::LPS: {
            if (K_ == 1)
                return x.rate * kf * std::min(mi_[ind], class_pop(X, ind, r)) *
                       lldfac(ist, class_pop(X, ind, r)) * cdfac(ist, class_counts(X, ind), r);
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * kf * (n[r] / (eps + tot)) * std::min(mi_[ind], eps + tot) *
                   lldfac(ist, tot) * cdfac(ist, n, r);
        }
        case SchedStrategy::DPS: {
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * kf * dpsshare(wnorm_[ist], n, r) * lldfac(ist, tot) *
                   cdfac(ist, n, r);
        }
        case SchedStrategy::GPS: {
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * kf * gpsshare(wnorm_[ist], n, r) * lldfac(ist, tot) *
                   cdfac(ist, n, r);
        }
        // The three PRIORITY sharing disciplines. Below the server count each
        // is its non-priority twin; above it only the most urgent non-empty
        // group shares the servers. Note that PSPRIO reads the lld factor at
        // the FULL population and the other two at the priority-restricted one
        // -- `prio_pop` carries that asymmetry, which comes from the reference.
        case SchedStrategy::PSPRIO: {
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * kf * psprioshare(n, r, mi_[ind]) * lldfac(ist, tot) *
                   cdfac(ist, n, r);
        }
        case SchedStrategy::DPSPRIO: {
            const std::vector<double> n = class_counts(X, ind);
            return x.rate * kf * dpsprioshare(wnorm_[ist], n, r, mi_[ind]) *
                   lldfac(ist, prio_pop(n, r, mi_[ind])) *
                   cdfac(ist, prio_vec(n, r, mi_[ind]), r);
        }
        case SchedStrategy::GPSPRIO: {
            const std::vector<double> n = class_counts(X, ind);
            return x.rate * kf * gpsprioshare(wnorm_[ist], n, r, mi_[ind]) *
                   lldfac(ist, prio_pop(n, r, mi_[ind])) *
                   cdfac(ist, prio_vec(n, r, mi_[ind]), r);
        }
        case SchedStrategy::FCFS:
        case SchedStrategy::LCFS:
        case SchedStrategy::SIRO:
        case SchedStrategy::HOL:
        case SchedStrategy::SEPT:
        case SchedStrategy::LEPT:
        case SchedStrategy::LCFSPR: {
            // Invariant: numel(buf) == max(0, total - mi), so the jobs in
            // service are the class population minus its share of the buffer.
            double waiting = 0.0;
            for (std::size_t c : bufs[ind])
                if (c == r) waiting += 1.0;
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * kf * std::max(0.0, n[r] - waiting) * lldfac(ist, tot) *
                   cdfac(ist, n, r);
        }
        case SchedStrategy::POLLING: {
            // One server serving one job, of the class its controller attends:
            // the class-r departure fires only while the controller SERVES r, at
            // the plain rate, never scaled by the class population.
            const PollCtrl& c = poll_ctrl_[ind];
            if (c.mode != 1 || c.pos != r + 1) return 0.0;
            const std::vector<double> n = class_counts(X, ind);
            double tot = 0.0;
            for (double v : n) tot += v;
            return x.rate * lldfac(ist, tot) * cdfac(ist, n, r);
        }
        case SchedStrategy::PAS:
        case SchedStrategy::OI:
            // Position p is served at Delta_mu(c1..cp) and pass-and-swap
            // decides which class that completion ejects, so the class-r rate
            // is the total increment over the positions that eject class r.
            return pas_rate(ist, bufs[ind], r);
        default:
            throw UnsupportedError(
                "solver_ssa_nrm: the scheduling policy '" +
                std::string(lang::sched_to_text(sched_[ist])) + "' at station '" +
                sn_.stations[ist].name + "' has no NRM rate law in this port");
    }
}
 
// ---------------------------------------------------------------------------
// Buffer maintenance
// ---------------------------------------------------------------------------
 
/** `pickFromBuffer`: the waiting job the discipline promotes, newest-first buffer. */
template <class T>
std::size_t NrmEngine<T>::pick_from_buffer(const std::vector<std::size_t>& buf, std::size_t ist) {
    switch (sched_[ist]) {
        case SchedStrategy::FCFS:
            return buf.size() - 1;  // oldest
        case SchedStrategy::LCFS:
        case SchedStrategy::LCFSPR:
            return 0;  // newest, or the most recently preempted
        case SchedStrategy::SIRO:
            return rng_.index(buf.size());
        case SchedStrategy::HOL: {
            // Highest priority (lowest value), FCFS within the group, so the
            // oldest is the LAST matching position.
            double best = std::numeric_limits<double>::infinity();
            for (std::size_t c : buf) best = std::min(best, classprio_[c]);
            for (std::size_t p = buf.size(); p-- > 0;)
                if (classprio_[buf[p]] == best) return p;
            return buf.size() - 1;
        }
        case SchedStrategy::SEPT:
        case SchedStrategy::LEPT: {
            // schedparam carries the rank of the class's mean service time,
            // ascending for SEPT and descending for LEPT, so both promote the
            // waiting class of least rank.
            double best = std::numeric_limits<double>::infinity();
            for (std::size_t c : buf)
                best = std::min(best, num_traits<T>::to_double(sn_.stations[ist].schedparam[c]));
            for (std::size_t p = buf.size(); p-- > 0;)
                if (num_traits<T>::to_double(sn_.stations[ist].schedparam[buf[p]]) == best)
                    return p;
            return buf.size() - 1;
        }
        default:
            throw UnsupportedError("solver_ssa_nrm: '" +
                                   std::string(lang::sched_to_text(sched_[ist])) +
                                   "' is not a buffered policy this port promotes from");
    }
}
 
/** `drawEntryPhase`: the phase a job starts service in, drawn from pie. */
template <class T>
std::size_t NrmEngine<T>::draw_entry_phase(std::size_t ind, std::size_t r) {
    if (nph_[ind][r] <= 1) return 0;
    return rng_.draw(pie_[ind][r]);
}
 
/**
 * `capacityLoss`: an OPEN arrival at a full physically-capped station is lost.
 *
 * A refused CLOSED job must block rather than vanish from the conserved
 * population, so it is not dropped here; `State.arrivalIsLost` draws that line
 * on the class type and this reproduces it. `isPhysicalCapacity` is the drop
 * rule test -- a bound that exists only as a state-space cutoff keeps the WAITQ
 * default and must not become a loss.
 */
template <class T>
bool NrmEngine<T>::capacity_loss(const std::vector<double>& X, std::size_t slot) const {
    const std::size_t jnd = slot_node_[slot], dst = slot_class_[slot];
    if (!is_station_[jnd]) return false;
    const std::size_t ist = to_station_[jnd];
    const qn::DropStrategy dr = sn_.droprule[ist][dst];
    if (dr == qn::DropStrategy::WAITQ) return false;  // cutoff-only, or unset
    if (dr == qn::DropStrategy::BAS || dr == qn::DropStrategy::BBS ||
        dr == qn::DropStrategy::RSRD)
        throw UnsupportedError("solver_ssa_nrm: station '" + sn_.stations[ist].name +
                               "' declares the blocking drop rule for class '" +
                               sn_.classes[dst].name +
                               "'; blocking-after-service is not ported to the C++ NRM");
    if (!std::isinf(sn_.classes[dst].population)) return false;  // closed: block, never lose
    const std::vector<double> cc = class_counts(X, jnd);
    double tot = 0.0;
    for (double v : cc) tot += v;
    if (!std::isinf(sn_.cap[ist]) && tot >= sn_.cap[ist]) return true;
    const double ccap = sn_.classcap[ist][dst];
    return ccap > 0.0 && !std::isinf(ccap) && cc[dst] >= ccap;
}
 
/**
 * `capacityBlock`: a refused arrival that may NOT be dropped cancels the firing.
 *
 * `capacity_loss` above answers the OPEN half of the same question and returns
 * false for a closed class precisely because a closed network's population is an
 * invariant. Nothing then stopped the reaction, so the NRM fired into the full
 * station anyway and reported the UNCONSTRAINED answer (BUG-81): on a closed
 * 3-queue tandem, N=6, Q2 capped at 2, QLen came back [1.96 2.07 1.97] against
 * the exact [3.609 0.971 1.420] -- a mean of 2.07 at a station that holds 2 --
 * while the serial engine, whose producer already implements the contract, gave
 * the exact answer.
 *
 * Blocking is the THIRD outcome of a firing, next to moving the job and losing
 * it. The departure does not occur, the source is not decremented, no buffer
 * moves; only the reaction's own clock is redrawn, which is exact by
 * memorylessness -- the residual of an exponential, or of the current PH phase,
 * is that same exponential. It is what `solver_ctmc` does when the target state
 * is absent from the enumerated space and the arc is dropped, which is why the
 * two agree.
 *
 * The population read is the PRE-arrival one, minus the departing job when it
 * currently sits at the destination node: a self-loop or a feedback arc at a
 * station already at cap would otherwise block itself forever, while the
 * reference producer sees the state AFTER the departure half.
 */
template <class T>
bool NrmEngine<T>::capacity_block(const std::vector<double>& X, std::size_t slot,
                                  std::size_t src_slot) const {
    if (!blk_on_) return false;
    const std::size_t jnd = slot_node_[slot], dst = slot_class_[slot];
    if (!is_station_[jnd]) return false;
    const std::size_t ist = to_station_[jnd];
    if (!blk_can_[ist][dst]) return false;
    const bool same_node = src_slot != npos && slot_node_[src_slot] == jnd;
    const std::vector<double> cc = class_counts(X, jnd);
    if (!std::isinf(sn_.cap[ist])) {
        double tot = 0.0;
        for (double v : cc) tot += v;
        if (same_node) tot -= 1.0;
        if (tot >= sn_.cap[ist]) return true;
    }
    const double ccap = sn_.classcap[ist][dst];
    if (ccap > 0.0 && !std::isinf(ccap)) {
        double pop = cc[dst];
        if (same_node && slot_class_[src_slot] == dst) pop -= 1.0;
        if (pop >= ccap) return true;
    }
    return false;
}
 
/** `pickPreempted`: the incumbent an LCFSPR arrival displaces, by server occupancy. */
template <class T>
std::size_t NrmEngine<T>::pick_preempted(const std::vector<double>& X,
                                         const std::vector<std::size_t>& buf, std::size_t jnd,
                                         std::size_t arr_class) {
    std::vector<double> insvc = class_counts(X, jnd);
    for (std::size_t r = 0; r < K_; ++r) {
        for (std::size_t c : buf)
            if (c == r) insvc[r] -= 1.0;
        if (r == arr_class) insvc[r] -= 1.0;  // the job that just arrived
        if (insvc[r] < 0.0) insvc[r] = 0.0;
    }
    double tot = 0.0;
    for (double v : insvc) tot += v;
    if (tot <= 0.0) return npos;
    const double u = rng_.uniform() * tot;
    double acc = 0.0;
    std::size_t last = npos;
    for (std::size_t r = 0; r < K_; ++r) {
        if (!(insvc[r] > 0.0)) continue;
        last = r;
        acc += insvc[r];
        if (u < acc) return r;
    }
    return last;
}
 
/** `applyArrivalBuffer`: join a just-arrived class-s job to node jnd's buffer. */
template <class T>
void NrmEngine<T>::apply_arrival_buffer(std::size_t jnd, std::size_t s,
                                        const std::vector<double>& X,
                                        std::vector<std::vector<std::size_t>>& bufs,
                                        std::vector<Matrix<double>>& svc, bool& svc_changed) {
    if (is_station_[jnd] && detail::sched_is_list(sched_[to_station_[jnd]])) {
        // PAS / OI: the arrival joins the BACK of the ordered list; there is no
        // server/buffer split. An arrival past the station's cap was already
        // lost or blocked by the capacity gates before it got here.
        const std::vector<std::size_t> old = bufs[jnd];
        bufs[jnd].push_back(s);
        pas_tag_started(jnd, old, bufs[jnd]);
        return;
    }
    // Polling: the controller, not the arrival, opens a service (`poll_land`).
    if (!is_poll_.empty() && is_poll_[jnd]) return;
    if (!is_station_[jnd] || !detail::sched_is_buffered(sched_[to_station_[jnd]])) {
        // A station that holds no waiting line serves every job it admits at
        // once (INF, PS, LPS, DPS, GPS), so the arrival IS the service start.
        if (is_station_[jnd]) tag(jnd, s, false);
        return;
    }
    const std::vector<double> cc = class_counts(X, jnd);
    double total = 0.0;
    for (double v : cc) total += v;
    bool entered = false;
    if (total > mi_[jnd]) {
        if (sched_[to_station_[jnd]] == SchedStrategy::LCFSPR) {
            // Preempt-resume: the arrival seizes a server and the incumbent it
            // displaces is the one that joins the buffer, which is what leaves
            // the new job in service (in-service is population minus buffer).
            const std::size_t c = pick_preempted(X, bufs[jnd], jnd, s);
            if (c != npos) {
                bufs[jnd].insert(bufs[jnd].begin(), c);
                tag(jnd, c, true);  // the incumbent is pushed back into the buffer
            }
            entered = true;
        } else {
            bufs[jnd].insert(bufs[jnd].begin(), s);  // it waits: it starts nothing
        }
    } else {
        entered = true;
    }
    if (entered) tag(jnd, s, false);  // the arrival seized a server
    if (entered && buf_ph_node_[jnd]) {
        svc[jnd](s, draw_entry_phase(jnd, s)) += 1.0;
        svc_changed = true;
    }
}
 
/**
 * `cacheAccess`: one READ at a cache, drawn from `after_event_cache`.
 *
 * The reference NRM re-implements State.afterEventCache's simulation branch;
 * this port draws from the enumerating port itself instead, so the NRM, the
 * serial engine and SolverCTMC read one set of replacement rules. The reading
 * job is presented ALONE (the enumerator requires exactly one job at the cache,
 * and other jobs at an immediate node are immaterial to a read), and the one
 * successor is drawn in proportion to its rate. Its per-class block, against the
 * one-hot input, says what happened:
 *
 *   all zero             the request MERGED onto an in-flight fetch (block B);
 *   a retrieval class    a fetch BEGINS: the in-flight bit is set, as the
 *                        departure that carries the job out would set it, and
 *                        the job is placed straight at the retrieval station,
 *                        as `cacheRetrDest` does, so that a job at
 *                        [cache, retrieval class] is always an inbound one;
 *   otherwise            a hit or a completed miss, plus the delayed hits the
 *                        completion releases, all at the cache in their classes.
 */
template <class T>
void NrmEngine<T>::cache_fire(const Rx& x, std::vector<double>& X, std::size_t& dest_pos,
                              bool& have_dest) {
    const std::size_t ind = x.node, r = x.cls;
    const qn::CacheParam<T>& cp = sn_.nodeparam.find(ind + 1)->second;
    std::vector<T> row(K_, num_traits<T>::from_int(0));
    row[r] = num_traits<T>::from_int(1);
    row.insert(row.end(), cache_var_[ind].begin(), cache_var_[ind].end());
    const qn::EventOutcome<T> o =
        qn::after_event_cache(sn_, ind + 1, row, lang::EventType::READ, r + 1);
    std::vector<double> w(o.space.size(), 0.0);
    double tot = 0.0;
    for (std::size_t j = 0; j < o.space.size(); ++j) {
        w[j] = num_traits<T>::to_double(o.rate[j]) * num_traits<T>::to_double(o.prob[j]);
        if (!(w[j] > 0.0)) w[j] = 0.0;
        tot += w[j];
    }
    if (!(tot > 0.0))
        throw NumericError("solver_ssa_nrm: the READ of class '" + sn_.classes[r].name +
                           "' at Cache '" + sn_.nodes[ind].name + "' has no enabled outcome");
    const std::vector<T>& nx = o.space[rng_.draw(w)];
    X[x.from] -= 1.0;
    cache_var_[ind].assign(nx.begin() + static_cast<std::ptrdiff_t>(K_), nx.end());
 
    std::vector<std::size_t> rcl, rci, rco;
    if (cp.retrieval_capacity > 0) qn::cache_retrieval_class_map(cp, rcl, rci, rco);
    const bool from_retr = std::find(rcl.begin(), rcl.end(), r + 1) != rcl.end();
    double produced = 0.0;
    for (std::size_t c = 0; c < K_; ++c) produced += num_traits<T>::to_double(nx[c]);
    if (produced == 0.0) {
        cache_dly_(ind, r) += 1.0;  // merged: held in block B until its fetch completes
        return;
    }
    for (std::size_t c = 0; c < K_; ++c) {
        const double d = num_traits<T>::to_double(nx[c]);
        if (!(d > 0.0)) continue;
        const bool begins = !from_retr && std::find(rcl.begin(), rcl.end(), c + 1) != rcl.end();
        if (!begins) {
            X[phoff_[ind][c]] += d;
            cache_prod_(ind, c) += d;
            if (!have_dest) {
                dest_pos = phoff_[ind][c];
                have_dest = true;
            }
            continue;
        }
        // BEGIN a fetch. The in-flight bit of the item whose retrieval class c is.
        std::size_t tcc = 0;
        for (int m : cp.itemcap)
            if (m > 0) tcc += static_cast<std::size_t>(m);
        for (std::size_t it = 0; it < cp.retrieval_classes.size(); ++it)
            if (r < cp.retrieval_classes[it].size() && cp.retrieval_classes[it][r] == c + 1 &&
                tcc + it < cache_var_[ind].size())
                cache_var_[ind][tcc + it] = num_traits<T>::from_int(1);
        // The retrieval station the class routes to, drawn over the routing row
        // and the entry phases there.
        std::vector<std::size_t> slots;
        std::vector<double> sw;
        for (std::size_t jnd = 0; jnd < I_; ++jnd)
            for (std::size_t s = 0; s < K_; ++s) {
                const double p =
                    num_traits<T>::to_double(sn_.rtnodes(ind * K_ + c, jnd * K_ + s));
                if (!(p > 0.0)) continue;
                if (buf_ph_class_[jnd][s]) {
                    slots.push_back(phoff_[jnd][s]);
                    sw.push_back(p);
                    continue;
                }
                for (std::size_t ke = 0; ke < nph_[jnd][s]; ++ke)
                    if (pie_[jnd][s][ke] > 0.0) {
                        slots.push_back(phoff_[jnd][s] + ke);
                        sw.push_back(p * pie_[jnd][s][ke]);
                    }
            }
        if (slots.empty())
            throw InputError("solver_ssa_nrm: retrieval class '" + sn_.classes[c].name +
                             "' at Cache '" + sn_.nodes[ind].name + "' routes nowhere");
        const std::size_t slot = slots[rng_.draw(sw)];
        X[slot] += 1.0;
        dest_pos = slot;
        have_dest = true;
    }
}
 
/**
 * The realized hit, miss and delayed-hit shares, as the serial engine reports
 * them: every read leaves as exactly one of the three, the released delayed
 * hits are carved out of the hit class's production, and the retrieval
 * classes, which have no hit class, get no share of their own.
 */
template <class T>
void NrmEngine<T>::cache_write_back() {
    cache_out_.clear();
    const double nan = std::numeric_limits<double>::quiet_NaN();
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (!is_cache_node_[ind]) continue;
        const qn::CacheParam<T>& cp = sn_.nodeparam.find(ind + 1)->second;
        SsaCacheRatio cr;
        cr.node = ind + 1;
        cr.hitprob.assign(K_, nan);
        cr.missprob.assign(K_, nan);
        cr.residt.assign(K_, nan);
        std::vector<double> dly(K_, 0.0);
        bool any_delayed = false;
        for (std::size_t k = 0; k < K_; ++k) {
            if (k >= cp.hitclass.size() || k >= cp.missclass.size()) continue;
            const std::size_t h = cp.hitclass[k], mi = cp.missclass[k];
            if (h == 0 || mi == 0 || h > K_ || mi > K_) continue;
            const double th = cache_prod_(ind, h - 1), tm = cache_prod_(ind, mi - 1);
            const double td = cache_dly_(ind, k);
            if (th + tm > 0.0) {
                cr.hitprob[k] = std::max(th - td, 0.0) / (th + tm);
                cr.missprob[k] = tm / (th + tm);
                dly[k] = td / (th + tm);
                if (td > 0.0) any_delayed = true;
            }
        }
        if (any_delayed) cr.delayedprob = dly;
        cache_out_.push_back(cr);
    }
}
 
/** `updateBuffers`: promote on a departure, join on an arrival. */
template <class T>
void NrmEngine<T>::update_buffers(std::size_t kfire, const std::vector<double>& X,
                                  std::vector<std::vector<std::size_t>>& bufs,
                                  std::size_t dest_pos, bool have_dest,
                                  std::vector<Matrix<double>>& svc, bool& svc_changed) {
    const Rx& x = rx_[kfire];
    const std::size_t ind = x.node;
    // The completing job leaves the phase it occupied; the promotion below
    // refills the freed server at a fresh entry phase.
    if (buf_ph_node_[ind] && x.is_buf_svc && !x.is_phase) {
        svc[ind](x.cls, x.dep_phase) -= 1.0;
        svc_changed = true;
    }
    // A PAS / OI departure is not a promotion but a pass-and-swap rewrite of the
    // whole list. Unlike the reference, which returns here, the ARRIVAL half
    // below still runs: a PAS station feeding a buffered one must join the job
    // to that station's buffer, and returning early would leave it counted in
    // the population but absent from the buffer.
    if (is_station_[ind] && detail::sched_is_list(sched_[to_station_[ind]]) &&
        !bufs[ind].empty() && !x.is_phase) {
        const std::vector<std::size_t> old = bufs[ind];
        pas_depart(ind, bufs[ind], x.cls);
        pas_tag_started(ind, old, bufs[ind]);
    } else if (is_station_[ind] && detail::sched_is_buffered(sched_[to_station_[ind]]) &&
        !bufs[ind].empty() && !x.is_phase) {
        const std::size_t pos = pick_from_buffer(bufs[ind], to_station_[ind]);
        const std::size_t promoted = bufs[ind][pos];
        bufs[ind].erase(bufs[ind].begin() + static_cast<std::ptrdiff_t>(pos));
        tag(ind, promoted, false);  // it takes the server the completion freed
        if (buf_ph_node_[ind]) {
            svc[ind](promoted, draw_entry_phase(ind, promoted)) += 1.0;
            svc_changed = true;
        }
    }
    if (have_dest)
        apply_arrival_buffer(slot_node_[dest_pos], slot_class_[dest_pos], X, bufs, svc,
                             svc_changed);
}
 
// ---------------------------------------------------------------------------
// The sample path
// ---------------------------------------------------------------------------
 
/** `next_reaction_method_direct`: the Gibson and Bruck loop with metrics in line. */
template <class T>
SsaSolution NrmEngine<T>::run() {
    const std::size_t nrx = rx_.size();
    SsaSolution out;
    out.QN = Matrix<double>(M_, K_, 0.0);
    out.UN = Matrix<double>(M_, K_, 0.0);
    out.RN = Matrix<double>(M_, K_, 0.0);
    out.TN = Matrix<double>(M_, K_, 0.0);
    out.CN.assign(K_, 0.0);
    out.XN.assign(K_, 0.0);
    out.StartN = Matrix<double>(M_, K_, 0.0);
    out.PreemptN = Matrix<double>(M_, K_, 0.0);
    out.method = "nrm";
    start_cnt_ = Matrix<double>(M_, K_, 0.0);
    preempt_cnt_ = Matrix<double>(M_, K_, 0.0);
    block_cnt_ = Matrix<double>(M_, K_, 0.0);
    if (nrx == 0) return out;
 
    std::vector<double> nvec = nvec0_;
    std::vector<std::vector<std::size_t>> bufs = buffers0_;
    std::vector<Matrix<double>> svc = svcph0_;
    // The WAITQ FIFO of each region, tokens `node * K + class` of the parked jobs.
    std::vector<std::vector<std::size_t>> fcr_buf(fcr_.size());
    cache_var_ = cache_var0_;
    cache_prod_ = Matrix<double>(I_, K_, 0.0);
    cache_dly_ = Matrix<double>(I_, K_, 0.0);
    // Seed each polling controller as `State.pollingInit` does: walk from buffer
    // 1 and settle on the first tangible state, so the seed is reachable.
    poll_ctrl_.assign(I_, PollCtrl());
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (!is_poll_[ind]) continue;
        std::size_t q = 1;
        int mode = 0;
        long budget = 0;
        qn::polling_next(poll_info_[ind], 1, poll_nbuf(nvec, ind), K_, true, q, mode, budget);
        poll_ctrl_[ind] = poll_land(ind, q, mode, budget, false);
    }
 
    std::vector<double> Ak(nrx, 0.0), Pk(nrx, 0.0), Tk(nrx, 0.0), tau(nrx, 0.0);
    for (std::size_t k = 0; k < nrx; ++k) {
        Ak[k] = propensity(k, nvec, bufs, svc);
        Pk[k] = -std::log(rng_.uniform());
        tau[k] = Ak[k] > 0.0 ? (Pk[k] - Tk[k]) / Ak[k] : std::numeric_limits<double>::infinity();
    }
 
    double total_time = 0.0;
    std::size_t n = 0;
    line::util::LineConsole::loop("drawing the sample path: %zu samples requested",
                                  static_cast<std::size_t>(opt_.samples));
    const std::size_t console_every = std::max<std::size_t>(1, opt_.samples / 20);
    for (; n < opt_.samples; ++n) {
        if ((n + 1) % console_every == 0)
            line::util::LineConsole::iter(
                static_cast<long>((n + 1) / console_every),
                "simulated %zu of %zu samples (%.0f%%), simulated time %.4g", n + 1,
                static_cast<std::size_t>(opt_.samples),
                100.0 * static_cast<double>(n + 1) / static_cast<double>(opt_.samples),
                total_time);
        std::size_t kfire = 0;
        double dt = std::numeric_limits<double>::infinity();
        for (std::size_t k = 0; k < nrx; ++k)
            if (tau[k] < dt) {
                dt = tau[k];
                kfire = k;
            }
        if (std::isinf(dt))
            throw NumericError(
                "solver_ssa_nrm: deadlock -- every reaction has propensity zero, so the sample "
                "path cannot advance");
        total_time += dt;
 
        // Time-integrate the metrics over this interval.
        for (std::size_t ist = 0; ist < M_; ++ist) {
            const std::size_t ind = sn_.station_to_node[ist] - 1;
            const std::vector<double> npop = class_counts(nvec, ind);
            double totpop = 0.0;
            for (double v : npop) totpop += v;
            // A Source's state slot is a FICTITIOUS TOKEN, not a queue. It is
            // seeded at 1, decremented by every arrival and replenished only
            // when a job reaches the Sink and is forwarded back, so its time
            // average is 1 - E[jobs in system] -- negative on any model with
            // more than one job in system, which is what it was reporting. A
            // Source holds no jobs, so QLen, Util and hence RespT/ResidT are
            // ZERO BY DEFINITION there and only its throughput (the arrival
            // rate) is a quantity; MVA, NC and MAM print exactly that row, and
            // the reference SSA applies the same rule downstream in
            // `solver_ssa.py` / `@@NetworkSolver/getAvg.m`. This is a
            // definition, not a clamp: it is keyed on the station's discipline
            // and never on the sign, so a negative anywhere else stays visible.
            const bool is_source = sched_[ist] == SchedStrategy::EXT;
            for (std::size_t k = 0; k < K_; ++k) {
                double dep = 0.0;
                for (std::size_t jd : dep_rx_[ist][k]) dep += Ak[jd];
                out.TN(ist, k) += dep * dt;
                if (is_source) continue;
                out.QN(ist, k) += npop[k] * dt;
                switch (sched_[ist]) {
                    case SchedStrategy::INF:
                        out.UN(ist, k) += npop[k] * dt;
                        break;
                    case SchedStrategy::PS:
                    case SchedStrategy::LPS:
                        if (totpop > 0.0)
                            out.UN(ist, k) += (npop[k] / totpop) *
                                              std::min(nservers_[ist], totpop) / nservers_[ist] *
                                              dt;
                        break;
                    case SchedStrategy::DPS:
                        out.UN(ist, k) += dpsshare(wnorm_[ist], npop, k) / nservers_[ist] * dt;
                        break;
                    case SchedStrategy::GPS:
                        out.UN(ist, k) += gpsshare(wnorm_[ist], npop, k) / nservers_[ist] * dt;
                        break;
                    // The utilization of a PRIORITY discipline is the SAME
                    // share the rate law uses, which is the invariant this
                    // accumulator exists to preserve: reusing the sharing
                    // factor is what keeps U = T E[S] from drifting away from
                    // the sample path that produced T.
                    case SchedStrategy::PSPRIO:
                        out.UN(ist, k) += psprioshare(npop, k, nservers_[ist]) / nservers_[ist] * dt;
                        break;
                    case SchedStrategy::DPSPRIO:
                        out.UN(ist, k) +=
                            dpsprioshare(wnorm_[ist], npop, k, nservers_[ist]) / nservers_[ist] * dt;
                        break;
                    case SchedStrategy::GPSPRIO:
                        out.UN(ist, k) +=
                            gpsprioshare(wnorm_[ist], npop, k, nservers_[ist]) / nservers_[ist] * dt;
                        break;
                    // Busy on one class-k job exactly while the controller serves k.
                    case SchedStrategy::POLLING:
                        if (poll_ctrl_[ind].mode == 1 && poll_ctrl_[ind].pos == k + 1)
                            out.UN(ist, k) += dt / nservers_[ist];
                        break;
                    // The positions with a positive Delta_mu over the servers, so a
                    // job served by several server types counts once.
                    case SchedStrategy::PAS:
                    case SchedStrategy::OI:
                        out.UN(ist, k) += pas_in_service(ist, bufs[ind], k) / nservers_[ist] * dt;
                        break;
                    case SchedStrategy::FCFS:
                    case SchedStrategy::LCFS:
                    case SchedStrategy::SIRO:
                    case SchedStrategy::HOL:
                    case SchedStrategy::SEPT:
                    case SchedStrategy::LEPT:
                    case SchedStrategy::LCFSPR: {
                        if (sn_.disabled[ist][k]) break;
                        double waiting = 0.0;
                        for (std::size_t c : bufs[ind])
                            if (c == k) waiting += 1.0;
                        out.UN(ist, k) += ((npop[k] - waiting) / nservers_[ist]) * dt;
                        break;
                    }
                    default:
                        break;
                }
            }
        }
 
        // Apply the firing.
        const Rx& x = rx_[kfire];
        std::size_t dest_pos = npos;
        bool have_dest = false;
        // A firing has THREE outcomes, not two: it moves the job, it loses it
        // (the source departs either way), or it is BLOCKED -- cancelled with
        // the source keeping the job and no slot changing. See `capacity_block`.
        bool blocked = false;
        if (x.is_cache) {
            cache_fire(x, nvec, dest_pos, have_dest);
        } else if (x.nnzP > 1 || !x.sd_node.empty()) {
            std::size_t slot;
            if (!x.sd_node.empty()) {
                // the dispatcher decides, from the state, not from a draw
                slot = resolve_state_dependent_dest(kfire, nvec);
            } else {
                const double u = rng_.uniform();
                std::size_t sel = x.cdf.size() - 1;
                for (std::size_t i = 0; i < x.cdf.size(); ++i)
                    if (x.cdf[i] > u) {
                        sel = i;
                        break;
                    }
                slot = x.to_slots[sel];
            }
            // A closed job that finds no room BLOCKS: the firing is cancelled
            // outright and nothing below runs, because a job that cannot leave
            // its station never departs at all. See `capacity_block`.
            if (capacity_block(nvec, slot, x.from)) {
                blocked = true;
            } else {
            // The loss is decided on the PRE-arrival population, so it is drawn
            // before the state is updated: the source still releases the job and
            // the destination never receives it.
            bool lost = capacity_loss(nvec, slot);
            // A region refuses on the same pre-arrival population; see
            // `fcr_outcome` for what each rule does with the refused job.
            if (!lost && !fcr_.empty())
                fcr_outcome(nvec, x, slot_node_[slot], slot_class_[slot], fcr_buf, lost, blocked);
            if (!blocked) {
            for (std::size_t s : x.from_slots) nvec[s] -= 1.0;
            if (!lost) {
                nvec[slot] += 1.0;
                dest_pos = slot;
                have_dest = true;
            }
            }
            }
        } else if (x.det_dest != npos && !x.is_phase &&
                   capacity_block(nvec, x.det_dest, x.from)) {
            blocked = true;   // same closed-class block; a phase change is never one
        } else {
            bool lost = x.det_dest != npos && !x.is_phase && capacity_loss(nvec, x.det_dest);
            // Single-destination departures cross region boundaries too; a phase
            // change is not an arrival and is never gated.
            if (!lost && x.det_dest != npos && !x.is_phase && !fcr_.empty())
                fcr_outcome(nvec, x, slot_node_[x.det_dest], slot_class_[x.det_dest], fcr_buf, lost,
                            blocked);
            if (blocked) {
                // censored: nothing moves, as for `capacity_block`
            } else if (lost) {
                nvec[x.from] -= 1.0;
            } else {
                for (std::size_t i = 0; i < NS_; ++i)
                    if (S_(i, kfire) != 0.0) nvec[i] += S_(i, kfire);
                if (x.det_dest != npos) {
                    dest_pos = x.det_dest;
                    have_dest = true;
                }
            }
        }
 
        bool svc_changed = false;
        if (blocked) {
            // nothing moved, so no buffer and no service phase may change
            if (is_station_[x.node]) block_cnt_(to_station_[x.node], x.cls) += 1.0;
        } else if (x.is_phase && buf_ph_node_[x.node]) {
            // A buffered-PH phase change frees no server and adds no arrival:
            // only the in-service multiset moves.
            svc[x.node](x.cls, x.phase_from) -= 1.0;
            svc[x.node](x.cls, x.phase_to) += 1.0;
            svc_changed = true;
        } else if (!x.is_sw) {
            update_buffers(kfire, nvec, bufs, dest_pos, have_dest, svc, svc_changed);
        }
        // A read rewrites the contents and places jobs outside the stoichiometry.
        if (x.is_cache) svc_changed = true;
 
        // Polling controller advance, `solver_ssa_nrm.m` lines 1710-1774: a
        // switchover step, a completion that ends or continues the visit, and an
        // arrival that wakes a parked server. Each can flip a service gate or the
        // switchover rate outside the static dependency set, so any move forces a
        // full refresh below.
        if (!is_poll_.empty() && !blocked) {
            if (x.is_sw) {
                const std::size_t ind = x.node;
                PollCtrl& c = poll_ctrl_[ind];
                const Matrix<T>& D0 = poll_info_[ind].sw_d0[c.pos - 1];
                const std::size_t Ksw = D0.rows();
                std::vector<double> w(Ksw + 1, 0.0);
                double rowsum = 0.0;
                for (std::size_t kd = 0; kd < Ksw; ++kd) {
                    const double v = num_traits<T>::to_double(D0(c.swk, kd));
                    rowsum += v;
                    if (kd != c.swk && v > 0.0) w[kd] = v;
                }
                w[Ksw] = std::max(0.0, -rowsum);  // absorption: the D1 row sum
                const std::size_t pick = rng_.draw(w);
                if (pick < Ksw && pick != c.swk) {
                    c.swk = pick;  // an internal phase advance of the walk
                } else {
                    std::size_t q = c.pos;
                    int mode = 0;
                    long budget = 0;
                    qn::polling_next(poll_info_[ind], c.pos, poll_nbuf(nvec, ind), K_, true, q,
                                     mode, budget);
                    c = poll_land(ind, q, mode, budget);
                }
                svc_changed = true;
            } else if (is_poll_[x.node] && !x.is_phase) {
                const std::size_t ind = x.node;
                PollCtrl& c = poll_ctrl_[ind];
                const qn::PollingInfo<T>& pi = poll_info_[ind];
                const std::vector<long> nb = poll_nbuf(nvec, ind);  // after the departure
                long ctrnext = 0;
                bool goon = false;
                switch (pi.ptype) {
                    case lang::PollingType::EXHAUSTIVE:
                        goon = nb[c.pos - 1] > 0;
                        break;
                    case lang::PollingType::GATED:
                        ctrnext = c.ctr - 1;
                        goon = ctrnext > 0;
                        break;
                    case lang::PollingType::KLIMITED:
                        ctrnext = c.ctr - 1;
                        goon = ctrnext > 0 && nb[c.pos - 1] > 0;
                        break;
                    case lang::PollingType::DECREMENTING:
                        ctrnext = c.ctr;
                        goon = nb[c.pos - 1] > c.ctr;
                        break;
                    default:
                        throw UnsupportedError("solver_ssa_nrm: unsupported polling type");
                }
                if (goon) {
                    c = poll_land(ind, c.pos, 1, ctrnext);
                } else {
                    std::size_t q = c.pos;
                    int mode = 0;
                    long budget = 0;
                    qn::polling_next(pi, c.pos, nb, K_, false, q, mode, budget);
                    c = poll_land(ind, q, mode, budget);
                }
                svc_changed = true;
            }
            if (have_dest && is_poll_[slot_node_[dest_pos]]) {
                const std::size_t jnd = slot_node_[dest_pos];
                PollCtrl& c = poll_ctrl_[jnd];
                if (c.mode == 0) {
                    std::size_t q = c.pos;
                    int mode = 0;
                    long budget = 0;
                    qn::polling_next(poll_info_[jnd], c.pos, poll_nbuf(nvec, jnd), K_, true, q,
                                     mode, budget);
                    c = poll_land(jnd, q, mode, budget);
                    svc_changed = true;
                }
            }
        }
 
        // WAITQ: admit parked jobs whose regions this firing may have relieved.
        // A release changes populations at arbitrary nodes, so it forces a full
        // propensity refresh below, as a service-phase move does.
        if (fcr_any_waitq_ && !blocked &&
            fcr_release_cascade(nvec, bufs, fcr_buf, svc, svc_changed) > 0)
            svc_changed = true;
 
        for (std::size_t k = 0; k < nrx; ++k) Tk[k] += Ak[k] * dt;
 
        // svcph is not part of the stoichiometry, so the static dependency set
        // cannot see a change to it; a move there forces a full refresh.
        if (svc_changed) {
            for (std::size_t k = 0; k < nrx; ++k) Ak[k] = propensity(k, nvec, bufs, svc);
        } else {
            for (std::size_t k : D_[kfire]) Ak[k] = propensity(k, nvec, bufs, svc);
        }
 
        Pk[kfire] -= std::log(rng_.uniform());
        for (std::size_t k = 0; k < nrx; ++k)
            tau[k] =
                Ak[k] > 0.0 ? (Pk[k] - Tk[k]) / Ak[k] : std::numeric_limits<double>::infinity();
    }
 
    if (total_time > 0.0)
        for (std::size_t ist = 0; ist < M_; ++ist)
            for (std::size_t k = 0; k < K_; ++k) {
                out.QN(ist, k) /= total_time;
                out.UN(ist, k) /= total_time;
                // net the blocked firings out of the departure-rate integral
                out.TN(ist, k) = (out.TN(ist, k) - block_cnt_(ist, k)) / total_time;
                // COUNTS over the simulated time, where TN above integrates a
                // rate: the two estimators of the same quantity agree in the
                // limit and differ by simulation error at any finite budget.
                out.StartN(ist, k) = start_cnt_(ist, k) / total_time;
                out.PreemptN(ist, k) = preempt_cnt_(ist, k) / total_time;
            }
 
    // A class- or joint-dependent station reports U = T*S/peak against the
    // product of its declared cd and jd peaks, the T*S/c convention of the
    // analytic solvers and serial SSA; busy time divided by the server count is
    // not that quantity once beta(n) != 1.
    for (std::size_t ist = 0; ist < M_; ++ist) {
        const auto& st = sn_.stations[ist];
        if (!st.cdscaling && !st.jdscaling) continue;
        for (std::size_t k = 0; k < K_; ++k) {
            double peak = 1.0;
            if (st.cdscaling) peak *= num_traits<T>::to_double(st.cdscalingpeak[k]);
            if (st.jdscaling) peak *= num_traits<T>::to_double(st.jdscalingpeak[k]);
            const double rate = num_traits<T>::to_double(sn_.rates(ist, k));
            out.UN(ist, k) = (std::isfinite(rate) && rate > 0.0 && peak > 0.0)
                                 ? out.TN(ist, k) / rate / peak
                                 : 0.0;
        }
    }
 
    // Every LOAD-DEPENDENT station reports the work-based T*S/peak instead,
    // peak = max(c, max(alpha)). The loop above integrates BUSY TIME, which is a
    // different quantity once alpha(n) != 1: a server running alpha(n) times
    // faster does the same work in less time, so busy time reads it as no busier
    // than one at its nominal rate. That put NRM at 0.9587 on a 4-job closed
    // model with alpha = [1 1.5 2 2.5] where CTMC, MVA, NC and serial SSA all
    // report 0.6612, and left the two SSA engines disagreeing with each other.
    // INF and EXT keep U = Q, as everywhere else.
    for (std::size_t ist = 0; ist < M_; ++ist) {
        if (lld_[ist].empty()) continue;
        double peak = nservers_[ist];
        bool non_unit = false;
        for (double a : lld_[ist]) {
            if (a != 1.0) non_unit = true;
            if (a > peak) peak = a;
        }
        if (!non_unit) continue;
        if (sched_[ist] == SchedStrategy::INF || sched_[ist] == SchedStrategy::EXT)
            continue;
        for (std::size_t k = 0; k < K_; ++k) {
            const double rate = num_traits<T>::to_double(sn_.rates(ist, k));
            out.UN(ist, k) = (std::isfinite(rate) && rate > 0.0 && peak > 0.0)
                                 ? out.TN(ist, k) / rate / peak
                                 : 0.0;
        }
    }
 
    for (std::size_t k = 0; k < K_; ++k) {
        out.XN[k] = out.TN(sn_.classes[k].refstat - 1, k);
        for (std::size_t ist = 0; ist < M_; ++ist)
            out.RN(ist, k) = out.TN(ist, k) > 0.0 ? out.QN(ist, k) / out.TN(ist, k) : 0.0;
        if (out.XN[k] > 0.0) out.CN[k] = sn_.classes[k].population / out.XN[k];
    }
    out.simulated_time = total_time;
    out.samples = n;
    cache_write_back();
    return out;
}
 
}  // namespace ssa
}  // namespace line
 
#endif  // LINE_SOLVERS_SSA_SOLVER_SSA_NRM_H