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/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;
}
 
/**
 * 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
 
    /** `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();
        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(); }
 
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_;
 
    // ---- construction ------------------------------------------------------
    void build_layout();
    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;
    }
 
    // ---- 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;
    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.
 */
template <class T>
void NrmEngine<T>::build_reactions() {
    rx_.clear();
    // 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];
                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);
                }
        }
    }
 
    // 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_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) 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);
    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;
                }
            }
        }
        // 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;
 
    // 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);
    }
 
    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);
        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));
            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);
        }
        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);
        }
        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);
        }
        // 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);
        }
        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]));
        }
        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]));
        }
        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);
        }
        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_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;
    }
}
 
/** `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;
    }
    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_;
 
    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;
                    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.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.
            const bool lost = capacity_loss(nvec, slot);
            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 {
            const bool lost =
                x.det_dest != npos && !x.is_phase && capacity_loss(nvec, x.det_dest);
            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 {
            update_buffers(kfire, nvec, bufs, dest_pos, have_dest, svc, svc_changed);
        }
 
        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;
            }
 
    // 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;
    return out;
}
 
}  // namespace ssa
}  // namespace line
 
#endif  // LINE_SOLVERS_SSA_SOLVER_SSA_NRM_H