solver_ssa_nrm_spn.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_SPN_H
#define LINE_SOLVERS_SSA_SOLVER_SSA_NRM_SPN_H
 
/**
 * @file
 * @ingroup line_solvers
 * SolverSSA on a STOCHASTIC PETRI NET: the port of `solver_ssa_nrm_spn`, the
 * sub-engine `solver_ssa_nrm.m:174` hands a net carrying Transition nodes to.
 *
 * WHY IT IS A SEPARATE ENGINE. The queueing reaction builder of
 * `solver_ssa_nrm.h` reads a state vector of JOBS PER (node, class, phase) and
 * a departure is the absorption of a service process. A Petri net has neither:
 * a Place holds TOKENS, a Transition mode is a reaction whose stoichiometry is
 * the arc incidence, and nothing in the net has a service phase. The reference
 * therefore branches on `any(sn.nodetype == NodeType.Transition)` before it
 * builds a single reaction, and so does this port.
 *
 * THE MAPPING IS EXACT, not an approximation. A Place holds a per-class token
 * count (one slot of the state vector), and a timed Transition mode is a
 * reaction: input (enabling) arcs consume, output (firing) arcs produce.
 * Enabling is a propensity gate -- every input place at or above its arc
 * weight, every inhibitor place strictly below its threshold -- and a
 * single-server mode fires at its exponential rate while an infinite or
 * k-server mode fires at that rate times its enabling degree. One firing
 * applies the stoichiometry once, which is the atomic GSPN firing the exact
 * CTMC, JMT and the standard GSPN tools all take.
 *
 * IMMEDIATE MODES ARE NOT REACTIONS. They fire in zero time, so they are
 * resolved by VANISHING-MARKING ELIMINATION: after every timed firing (and once
 * on the initial marking) every enabled immediate mode is fired -- highest
 * firing priority first and, among equal priority, drawn in proportion to
 * firing weight -- until the marking is tangible. The timed race therefore only
 * ever samples from tangible markings and the immediate modes consume no
 * simulated time.
 *
 * A SOURCE IS NOT A TRANSITION and needs a reaction of its own, or the Place it
 * feeds stays empty and the net deadlocks on the first draw. Splitting a
 * Poisson stream by independent routing probabilities yields independent
 * Poisson streams, so the edge of probability p carries rate lambda*p exactly;
 * the reaction has an EMPTY enabling set (a constant propensity) and deposits
 * one token into the routed Place slot.
 *
 * WHAT IT REFUSES, and each refusal is the reference's: a non-exponential timed
 * firing (representing an in-flight firing's phase would need per-mode phase
 * state the reaction network does not carry), a non-exponential Source arrival
 * into a Place, a Source that reaches no Place, an infinite initial marking, and
 * a marking-dependent firing rate -- the last under EVERY SSA method, since no
 * SSA engine applies the g(marking) multiplier and answering with the nominal
 * rate would be silently wrong. `spn_nrm_supported` is the predicate; the
 * `default` dispatch reads it to decide whether to prefer this engine, and the
 * explicit `nrm` arm reads it to refuse by name.
 *
 * THE SLOT LAYOUT IS FLAT, one slot per (node, class), where the queueing
 * engine's is per (node, class, PHASE). Nothing here has a phase: a Place holds
 * tokens rather than jobs in service, and the reference's own SPN path only
 * ever addresses `phOff(place, class) + 1`, the first phase slot of the pair.
 * The two layouts therefore agree on every slot this engine touches.
 *
 * NOT A SAMPLE-PATH TWIN OF THE REFERENCE. The random stream is this port's
 * MT19937 (`SsaRng`) and the draw order is the algorithm's, so a seeded run
 * matches the reference STATISTICALLY and never bit for bit, exactly as the
 * queueing NRM does.
 *
 * A FLUID NET (FSPN, `Place::set_fluid`) runs a PIECEWISE-DETERMINISTIC loop instead, the reference's
 * `spnFluidRun`. Between two firings the marking is constant, so every fluid level moves LINEARLY at the slope
 * its running flows give it: a flow of a timed mode counts once per running server, min(enabling degree,
 * servers), and only while every guard of the mode holds. Three things end a segment and the first one wins:
 * a FIRING (a Gibson-Bruck clock runs out; a mode whose rate depends on a level has a hazard that varies along
 * the segment, so its clock is integrated by adaptive quadrature and the epoch is a root of the cumulative
 * hazard), a GUARD crossing (nothing fires; the guards are re-read and the immediate modes the crossing enabled
 * fire in zero time), or a BOUND hit at 0 or the capacity (nothing fires; the slope that would carry the level
 * out is cut from then on). A firing applies the token stoichiometry and then the mode's IMPULSES, each clipped
 * at the bounds. A guard is inclusive and read with a LOOKAHEAD at the threshold: it holds there only if the
 * slope keeps it holding, and a guarded flow that pushes its own level back across its own guard (a sliding
 * mode) is refused by name. A fluid place reports its time-average level as QN and everything that leaves it
 * as TN. `ssa_fluid_refusal` is the gate: a fluid net has this one engine, so `serial` is refused.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/sn/sn_fluid_refusal.h"
#include "line/lang/qn/network_struct.h"
#include "line/solvers/ssa/ssa_types.h"
#include "line/util/error.h"
#include "line/util/line_console.h"
#include "line/util/matrix.h"
#include "line/util/rootfind.h"
 
namespace line {
namespace ssa {
 
namespace detail {
 
/**
 * One reaction of the net: a timed Transition mode, an immediate one, or a
 * Source arrival. The reference's `spnEmptyRx` record, field for field.
 */
struct SpnReaction {
    std::size_t node = 0;  ///< 1-based node the reaction belongs to
    std::size_t mode = 0;  ///< 1-based mode; 0 marks a Source arrival
    /** Stoichiometry over the (node, class) slots: consume negative, produce positive. */
    std::vector<double> S;
    std::vector<std::size_t> en_slot;   ///< input arcs: the slots read
    std::vector<double> en_w;           ///< and the tokens each one needs
    std::vector<std::size_t> inh_slot;  ///< inhibitor arcs: the slots read
    std::vector<double> inh_thr;        ///< and the count at which each blocks
    double base_rate = 0.0;             ///< exponential firing rate, sum(D1)
    double nservers = 1.0;              ///< concurrent firings the mode may run
    double weight = 1.0;                ///< immediate: weight among equal priorities
    double prio = 1.0;                  ///< immediate: firing priority
    std::vector<std::size_t> dep_slots;  ///< slots this reaction DEPOSITS into
    // FSPN arcs, each fluid place named by its fluid index q (0-based)
    std::vector<std::size_t> flow_q;    ///< flows: the level each one moves
    std::vector<double> flow_rate;      ///< and its signed rate per running server
    std::vector<std::size_t> guard_q;   ///< guards: the level each one reads
    std::vector<int> guard_sgn;         ///< +1 for '>=', -1 for '<='
    std::vector<double> guard_thr;      ///< and the threshold
    std::vector<std::size_t> imp_q;     ///< impulses: the level a firing moves
    std::vector<double> imp_d;          ///< and by how much
    std::vector<std::size_t> fdep_q;    ///< level-dependent rate: the level read
    std::vector<std::function<double(double)>> fdep_fun;  ///< and the multiplier it gives
    /** Marking-dependent flow multiplier per flow (aligned with flow_q; empty: none), on the flat I*K marking. */
    std::vector<std::function<double(const std::vector<double>&)>> flow_fun;
    std::vector<std::vector<std::size_t>> flow_read;  ///< and the 0-based token places each one reads
    /** Marking-dependent firing multiplier g on the flat (node, class) marking of `gdep_read`, as `state_events.h`. */
    std::function<double(const std::vector<double>&)> gdep;
    std::vector<std::size_t> gdep_read;  ///< 0-based input places g reads
};
 
/**
 * Enabling degree: how many concurrent firings the marking supports, the
 * minimum over the input arcs of floor(tokens / weight), zeroed by any active
 * inhibitor. A mode with no input arc is single-degree -- which is what makes a
 * Source arrival a constant-propensity reaction.
 */
inline double spn_en_degree(const std::vector<double>& n, const SpnReaction& rx) {
    for (std::size_t i = 0; i < rx.inh_slot.size(); ++i)
        if (n[rx.inh_slot[i]] >= rx.inh_thr[i]) return 0.0;
    if (rx.en_slot.empty()) return 1.0;
    double d = std::numeric_limits<double>::infinity();
    for (std::size_t i = 0; i < rx.en_slot.size(); ++i)
        d = std::min(d, std::floor(n[rx.en_slot[i]] / rx.en_w[i]));
    return d;
}
 
/**
 * Propensity of a timed mode: the exponential rate times the effective server
 * count, min(enabling degree, mode servers). A single-server mode therefore
 * fires at its rate whenever enabled and an infinite-server one at the rate
 * scaled by the enabling degree.
 */
inline double spn_propensity(const std::vector<double>& n, const SpnReaction& rx) {
    const double eff = std::min(spn_en_degree(n, rx), rx.nservers);
    return eff <= 0.0 ? 0.0 : rx.base_rate * eff;
}
 
/**
 * `ssa_firingdep_refusal.m`: the NRM cannot answer a net whose firing rates
 * depend on the marking, and says so by name.
 *
 * The NRM builds one CONSTANT-propensity reaction per timed mode, so the
 * g(marking) multiplier has no place in it and answering with the nominal rate
 * would be a wrong number under the caller's own method name.
 *
 * THE SAME SCOPE AS THE REFERENCE'S. The refusal binds the NRM only, because
 * the SERIAL engine reaches the multiplier for real -- the FIRE arm of
 * `after_global_event` (`state_events.h`) scales the mode's firing rate by
 * `firingdep[m](marking)` on the state it is leaving, which is exactly the
 * semantics SolverCTMC gives it. The gate therefore reads it with `raise=false`,
 * so `default` FALLS BACK to the serial engine, and with `raise=true` only when
 * the caller named `nrm`. MATLAB, the JAR and Python refused it under every SSA
 * method until 2026-10-18; they now share this rule.
 *
 * A FLUID NET IS EXEMPT. Its runner (`run_fluid`) recomputes every propensity at each segment, and the marking is
 * constant within a segment, so it applies g(marking) exactly, on the input-place marking the serial engine reads.
 */
template <class T>
bool ssa_check_firingdep(const qn::NetworkStruct<T>& sn, bool raise = true) {
    if (sn.has_fluid_places()) return true;
    for (typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
             sn.transparam.begin();
         it != sn.transparam.end(); ++it) {
        const qn::TransitionParam<T>& tp = it->second;
        for (std::size_t m = 0; m < tp.nmodes && m < tp.firingdep.size(); ++m) {
            if (!tp.firingdep[m]) continue;
            if (!raise) return false;
            const std::size_t ind = it->first;
            throw UnsupportedError(ErrorKind{"MarkingDependence"},
                "SolverSSA: transition '" +
                (ind >= 1 && ind <= sn.nodes.size() ? sn.nodes[ind - 1].name : std::string("?")) +
                "' mode " + std::to_string(m + 1) +
                " uses a marking-dependent firing rate (setFiringRateDependence), which the "
                "NRM does not support: it builds one constant-propensity reaction per timed "
                "mode and never evaluates the handle. Ask for method='serial', which applies "
                "it, or use SolverCTMC or SolverLDES");
        }
    }
    return true;
}
 
/**
 * `ssa_nrm_guards.m`'s `g.spn`: the shapes the SPN reaction builder can read.
 *
 * FOUR CONDITIONS, and all four are the reference's. Each one the reference
 * RAISES from inside the builder, and each is a condition the serial engine
 * serves, so `default` falls back rather than refusing -- which is why this is
 * asked with `raise = false` by the eligibility test and with `raise = true`
 * only when the caller named `nrm` itself.
 */
template <class T>
bool spn_nrm_supported(const qn::NetworkStruct<T>& sn, bool raise = true) {
    const std::size_t I = sn.nodes.size();
    const std::size_t K = sn.nclasses;
    bool any_transition = false;
    for (std::size_t i = 0; i < I; ++i)
        if (sn.nodes[i].nodetype == qn::NodeType::Transition) any_transition = true;
    if (!any_transition) return true;
    if (!ssa_check_firingdep(sn, raise)) return false;
 
    for (std::size_t i = 0; i < I; ++i) {
        if (sn.nodes[i].nodetype == qn::NodeType::Transition) {
            const typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
                sn.transparam.find(i + 1);
            if (it == sn.transparam.end()) continue;
            const qn::TransitionParam<T>& tp = it->second;
            for (std::size_t m = 0; m < tp.nmodes; ++m) {
                if (m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE)
                    continue;
                const bool one_phase = m < tp.firingphases.size() && tp.firingphases[m] == 1;
                const bool has_proc = m < tp.firingproc.size() && tp.firingproc[m].D1.rows() > 0;
                if (one_phase && has_proc) continue;
                if (!raise) return false;
                throw UnsupportedError(ErrorKind{"NonExponential"},
                    "SolverSSA(method='nrm'): transition '" + sn.nodes[i].name + "' mode " +
                    std::to_string(m + 1) +
                    " has a non-exponential timed firing. The reaction network carries one "
                    "constant-rate reaction per mode and no per-mode phase, so an in-flight "
                    "firing has nowhere to keep its phase; use method='serial'");
            }
        } else if (sn.nodes[i].nodetype == qn::NodeType::Source) {
            const std::size_t ist = sn.nodes[i].station;
            if (ist == 0) continue;
            for (std::size_t r = 0; r < K; ++r) {
                const double lambda = num_traits<T>::to_double(sn.rates(ist - 1, r));
                if (std::isnan(lambda) || lambda <= 0.0) continue;
                if (sn.service[ist - 1][r].type != lang::ProcessType::EXP) {
                    if (!raise) return false;
                    throw UnsupportedError(ErrorKind{"NonExponential"},
                        "SolverSSA(method='nrm'): source '" + sn.nodes[i].name + "' class '" +
                        sn.classes[r].name +
                        "' has a non-exponential arrival. An arrival into a Place is one "
                        "constant-propensity reaction here, which only a Poisson stream is; "
                        "use method='serial'");
                }
                bool feeds = false;
                for (std::size_t j = 0; j < I && !feeds; ++j) {
                    if (sn.nodes[j].nodetype != qn::NodeType::Place) continue;
                    for (std::size_t s = 0; s < K && !feeds; ++s)
                        if (num_traits<T>::to_double(sn.rtnodes(i * K + r, j * K + s)) > 0.0)
                            feeds = true;
                }
                if (!feeds) {
                    if (!raise) return false;
                    throw UnsupportedError(ErrorKind{"Topology"},
                        "SolverSSA(method='nrm'): source '" + sn.nodes[i].name + "' class '" +
                        sn.classes[r].name +
                        "' routes to no Place. The SPN path deposits an arrival into a Place "
                        "slot and has nowhere else to put one; use method='serial'");
                }
            }
        } else if (sn.nodes[i].nodetype == qn::NodeType::Place) {
            const typename std::map<std::size_t, std::vector<T>>::const_iterator im =
                sn.initmarking.find(i + 1);
            if (im == sn.initmarking.end()) continue;
            for (std::size_t r = 0; r < im->second.size(); ++r) {
                if (std::isfinite(num_traits<T>::to_double(im->second[r]))) continue;
                if (!raise) return false;
                throw UnsupportedError(ErrorKind{"Place"}, "SolverSSA(method='nrm'): place '" + sn.nodes[i].name +
                                       "' declares an infinite initial marking, which the "
                                       "reaction network cannot count; use method='serial'");
            }
        }
    }
    return true;
}
 
/**
 * `ssa_fluid_refusal.m`: whether SolverSSA can simulate this Fluid Stochastic Petri Net under `method`.
 *
 * A FLUID NET RUNS ON THE NRM PETRI RUNNER ONLY. Its levels move continuously between firings, so the
 * simulator is a piecewise-deterministic Markov process whose levels, guard crossings and bound hits ride
 * beside the NRM clocks as scheduled non-firing events. The serial engine has no level to move, so `serial`
 * (and its alias `ssa`) is refused by name and `default`/`parallel` run the NRM. A net the reaction builder
 * cannot read (a non-exponential firing) has no engine left and is refused with the builder's reason. A
 * marking-dependent firing or flow rate is read by the runner, segment by segment. Always OK for a net with no fluid.
 *
 * @param sn     the network structure
 * @param method the SSA method, without an `ssa.` prefix
 */
template <class T>
api::SnFluidVerdict ssa_fluid_refusal(const qn::NetworkStruct<T>& sn, const std::string& method) {
    api::SnFluidVerdict v;
    if (!sn.has_fluid_places()) return v;
    v = api::sn_fluid_refusal(sn);
    if (!v.ok) {
        v.reason = "SolverSSA: " + v.reason;
        return v;
    }
    if (method == "serial" || method == "ssa") {
        v.ok = false;
        v.reason = "SolverSSA: the model is a fluid stochastic Petri net (Place.setFluid), whose levels only the "
                   "NRM Petri runner moves. Use options.method='nrm' or 'default'.";
        return v;
    }
    for (const qn::NodeDef& nd : sn.nodes)
        if (nd.nodetype != qn::NodeType::Place && nd.nodetype != qn::NodeType::Transition &&
            nd.nodetype != qn::NodeType::Source && nd.nodetype != qn::NodeType::Sink) {
            v.ok = false;
            v.reason = "SolverSSA: a fluid stochastic Petri net may hold only Place, Transition, Source and Sink "
                       "nodes; the NRM Petri runner has no hand-off between a level and a queueing station.";
            return v;
        }
    if (!spn_nrm_supported(sn, false)) {
        v.ok = false;
        v.reason = "SolverSSA: the fluid net has a non-exponential firing or Source arrival, or an infinite "
                   "marking, which the NRM Petri runner (the only engine that moves fluid levels) has no form for.";
    }
    return v;
}
 
/**
 * `spnGuardHolds`: every fluid guard of a mode, inclusive, with the lookahead of the slope `v` at the
 * threshold itself, so a crossing switches the guard at the instant of the crossing.
 */
inline bool spn_fluid_guard_holds(const SpnReaction& rec, const std::vector<double>& f,
                                  const std::vector<double>& v) {
    for (std::size_t g = 0; g < rec.guard_q.size(); ++g) {
        const std::size_t q = rec.guard_q[g];
        const double sg = static_cast<double>(rec.guard_sgn[g]);
        const double thr = rec.guard_thr[g];
        const double tolf = 1e-12 * std::max(1.0, std::fabs(thr));
        const double d = sg * (f[q] - thr);
        if (d > tolf) continue;
        if (d >= -tolf && sg * v[q] >= 0.0) continue;
        return false;
    }
    return true;
}
 
/** `spnFluidDep`: the level-dependent multiplier of a mode at time `t` into the segment. */
inline double spn_fluid_dep(const SpnReaction& rec, const std::vector<double>& f, const std::vector<double>& v,
                            double t) {
    double h = 1.0;
    for (std::size_t i = 0; i < rec.fdep_q.size(); ++i) {
        const std::size_t q = rec.fdep_q[i];
        h *= rec.fdep_fun[i](std::max(0.0, f[q] + v[q] * t));
    }
    return h;
}
 
/** `spnFluidVarying`: does the hazard of a mode change along the segment? */
inline bool spn_fluid_varying(const SpnReaction& rec, const std::vector<double>& v) {
    for (std::size_t i = 0; i < rec.fdep_q.size(); ++i)
        if (v[rec.fdep_q[i]] != 0.0) return true;
    return false;
}
 
/** One adaptive Simpson panel: recurse until the Richardson estimate meets `tol` or the depth runs out. */
template <class F>
double spn_simpson_panel(const F& h, double a, double b, double fa, double fm, double fb, double whole,
                         double tol, int depth) {
    const double m = 0.5 * (a + b);
    const double lm = 0.5 * (a + m), rm = 0.5 * (m + b);
    const double flm = h(lm), frm = h(rm);
    const double left = (m - a) / 6.0 * (fa + 4.0 * flm + fm);
    const double right = (b - m) / 6.0 * (fm + 4.0 * frm + fb);
    const double err = left + right - whole;
    if (depth <= 0 || std::fabs(err) <= 15.0 * tol) return left + right + err / 15.0;
    return spn_simpson_panel(h, a, m, fa, flm, fm, left, 0.5 * tol, depth - 1) +
           spn_simpson_panel(h, m, b, fm, frm, fb, right, 0.5 * tol, depth - 1);
}
 
/**
 * `spnFluidHazardIntegral`: the cumulative hazard of a mode over [0, dt] of the segment. Closed form when the
 * hazard is constant; otherwise adaptive Simpson at the reference's `integral` tolerances (AbsTol 1e-12,
 * RelTol 1e-10).
 */
inline double spn_fluid_hazard(const SpnReaction& rec, double a0, const std::vector<double>& f,
                               const std::vector<double>& v, double dt) {
    if (dt <= 0.0) return 0.0;
    if (!spn_fluid_varying(rec, v)) return a0 * spn_fluid_dep(rec, f, v, 0.0) * dt;
    const auto h = [&](double t) { return spn_fluid_dep(rec, f, v, t); };
    const double fa = h(0.0), fm = h(0.5 * dt), fb = h(dt);
    const double whole = dt / 6.0 * (fa + 4.0 * fm + fb);
    const double tol = std::max(1e-12, 1e-10 * std::fabs(whole));
    return a0 * spn_simpson_panel(h, 0.0, dt, fa, fm, fb, whole, tol, 50);
}
 
/**
 * `spnFluidFireTime`: when the cumulative hazard of a mode reaches `need`, or infinity if not before
 * `horizon`. A constant hazard is closed form; a varying one is bracketed (doubling when the horizon is
 * infinite, at most 200 times) and the epoch found by Brent's method, the reference's `fzero`.
 */
inline double spn_fluid_fire_time(const SpnReaction& rec, double a0, const std::vector<double>& f,
                                  const std::vector<double>& v, double need, double horizon) {
    const double inf = std::numeric_limits<double>::infinity();
    if (!spn_fluid_varying(rec, v)) {
        const double a = a0 * spn_fluid_dep(rec, f, v, 0.0);
        return a <= 0.0 ? inf : need / a;
    }
    if (need <= 0.0) return 0.0;
    double hi;
    if (std::isfinite(horizon)) {
        hi = horizon;
        if (spn_fluid_hazard(rec, a0, f, v, hi) < need) return inf;
    } else {
        const double a = std::max(a0 * spn_fluid_dep(rec, f, v, 0.0), std::numeric_limits<double>::min());
        hi = need / a;
        int grow = 0;
        while (spn_fluid_hazard(rec, a0, f, v, hi) < need) {
            hi *= 2.0;
            if (++grow > 200) return inf;
        }
    }
    const auto g = [&](double t) { return spn_fluid_hazard(rec, a0, f, v, t) - need; };
    return root_brent<double>(g, 0.0, hi, 1e-14 * std::max(1.0, hi)).root;
}
 
}  // namespace detail
 
/**
 * The SPN Next-Reaction-Method engine, `solver_ssa_nrm_spn` of the reference.
 *
 * DOUBLE ONLY, for the reason `NrmEngine` is: the sample path is generated from
 * exponential clocks, which are logarithms of uniform draws, so there is no
 * exact value a wider arithmetic could carry. The analyzer refuses a non-double
 * backend by name before this class is instantiated.
 */
template <class T>
class NrmSpnEngine {
public:
    NrmSpnEngine(const qn::NetworkStruct<T>& sn, const SsaOptions& opt)
        : sn_(sn), opt_(opt), rng_(opt.seed) {
        build();
    }
 
    /** Run `opt.samples` firings and return the time-averaged metrics. */
    SsaSolution run();
 
    /** The reaction count, for the tests that assert the builder. */
    std::size_t nreactions() const { return rx_.size(); }
    /** The immediate-mode count, likewise. */
    std::size_t nimmediate() const { return imm_.size(); }
 
private:
    using Rx = detail::SpnReaction;
 
    const qn::NetworkStruct<T>& sn_;
    SsaOptions opt_;
    SsaRng rng_;
 
    std::size_t I_ = 0, K_ = 0, M_ = 0, NS_ = 0;
 
    std::vector<Rx> rx_;   ///< the timed reactions, the ones that race
    std::vector<Rx> imm_;  ///< the immediate modes, resolved between races
    /** Timed reactions consuming from each (0-based node, class): the Place throughput. */
    std::vector<std::vector<std::vector<std::size_t>>> consumers_;
    /** Arrival reactions injecting at each (0-based Source node, class). */
    std::vector<std::vector<std::vector<std::size_t>>> producers_;
 
    std::vector<double> nvec0_;  ///< the initial marking, before the collapse
    std::vector<double> pcap_slot_;  ///< per-(place, class) capacity, infinite when none
    /** Per-place total capacities: the bound and the slots it is taken over. */
    std::vector<std::pair<double, std::vector<std::size_t>>> place_total_caps_;
    bool has_caps_ = false;
 
    // FSPN fluid places (`fl` of the reference), indexed by the fluid index q
    std::vector<long> fq_;               ///< per 0-based node: its fluid index, -1 when not fluid
    std::vector<std::size_t> fl_node_;   ///< 0-based node of each fluid place
    std::vector<std::size_t> fl_ist_;    ///< 0-based station of each fluid place
    std::vector<double> fl_cap_;         ///< level ceiling, infinite when unbounded
    std::vector<double> fl_f0_;          ///< level at time zero
    /** Per timed reaction, the marking multiplier of each of its flows for the current segment (empty: all 1). */
    std::vector<std::vector<double>> fm_;
 
    /** The livelock guard of the vanishing-marking collapse, the reference's 1e5. */
    static const std::size_t kMaxImmSteps = 100000;
 
    std::size_t slot(std::size_t node0, std::size_t cls) const { return node0 * K_ + cls; }
 
    void build();
    Rx build_mode(std::size_t ind0, std::size_t m) const;
    void apply_caps(std::vector<double>& n, const std::vector<std::size_t>& deposited) const;
    void collapse(std::vector<double>& n);
 
    // the piecewise-deterministic loop of a fluid net, `spnFluidRun` and its helpers
    SsaSolution run_fluid();
    std::size_t fluid_index(std::size_t node1) const;
    void fluid_mark_mul(const std::vector<double>& n, std::vector<double>& gk);
    void fluid_resolve(const std::vector<double>& n, const std::vector<double>& f, std::vector<double>& eff,
                       std::vector<double>& v, std::vector<double>& vin, std::vector<double>& vout) const;
    void fluid_slope(const std::vector<double>& f, const std::vector<double>& eff, std::vector<double>& v,
                     std::vector<double>& vin, std::vector<double>& vout) const;
    void fluid_next_det(const std::vector<double>& f, const std::vector<double>& v, double& tdet,
                        std::size_t& qdet, double& fdet) const;
    void fluid_impulse(std::vector<double>& f, const Rx& rec, std::vector<double>& OF) const;
    void collapse_fluid(std::vector<double>& n, std::vector<double>& f, const std::vector<double>& v,
                        std::vector<double>& OF);
};
 
/**
 * The reaction list, the initial marking and the capacity tables.
 *
 * ONE PASS PER NODE KIND, in the reference's order: the Transition modes become
 * reactions (timed) or immediate records, then the Source arrivals become
 * constant-propensity reactions, then the Places contribute the marking and
 * their capacities.
 */
template <class T>
void NrmSpnEngine<T>::build() {
    I_ = sn_.nodes.size();
    K_ = sn_.nclasses;
    M_ = sn_.nstations;
    NS_ = I_ * K_;
    consumers_.assign(I_, std::vector<std::vector<std::size_t>>(K_));
    producers_.assign(I_, std::vector<std::vector<std::size_t>>(K_));
 
    // Refused here as well as in the gate, as the reference raises it from
    // inside its own builder too (`solver_ssa_nrm.m:3426`): that is the path a
    // caller reaches with the checks disabled, and this engine must not run a
    // net whose rates it would silently ignore however it was entered.
    detail::ssa_check_firingdep(sn_, true);
 
    // Fluid places (FSPN): one real level each, indexed q = fq_[node]. The structural rules are re-read here
    // for the same reason as the firingdep refusal above: a builder entered with the gate skipped must not run.
    fq_.assign(I_, -1);
    if (sn_.has_fluid_places()) {
        const api::SnFluidVerdict fv = api::sn_fluid_refusal(sn_);
        if (!fv.ok) throw UnsupportedError(ErrorKind{"Structure"}, "solver_ssa_nrm_spn: " + fv.reason);
        for (std::size_t ind = 0; ind < I_; ++ind) {
            if (sn_.nodes[ind].nodetype != qn::NodeType::Place || !sn_.is_fluid_node(ind + 1)) continue;
            const std::size_t ist = sn_.nodes[ind].station;
            const qn::Station<T>& st = sn_.stations[ist - 1];
            fq_[ind] = static_cast<long>(fl_node_.size());
            fl_node_.push_back(ind);
            fl_ist_.push_back(ist - 1);
            fl_cap_.push_back(st.fluidcap);
            fl_f0_.push_back(st.fluidinit);
        }
    }
 
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (sn_.nodes[ind].nodetype != qn::NodeType::Transition) continue;
        const typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
            sn_.transparam.find(ind + 1);
        if (it == sn_.transparam.end()) continue;
        const qn::TransitionParam<T>& tp = it->second;
        for (std::size_t m = 0; m < tp.nmodes; ++m) {
            const Rx rec = build_mode(ind, m);
            if (m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE) {
                imm_.push_back(rec);
                continue;
            }
            rx_.push_back(rec);
            const std::size_t ridx = rx_.size() - 1;
            for (std::size_t a = 0; a < rec.en_slot.size(); ++a)
                consumers_[rec.en_slot[a] / K_][rec.en_slot[a] % K_].push_back(ridx);
        }
    }
 
    // Source arrivals. Thinning a Poisson stream by the independent routing
    // probabilities leaves independent Poisson streams, so one reaction per
    // routed (Source, class) -> (Place, class) edge carries rate lambda*p
    // exactly. `producers_` is what makes the Source station report that rate as
    // its throughput, which is the open class's reference-station throughput.
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (sn_.nodes[ind].nodetype != qn::NodeType::Source) continue;
        const std::size_t ist = sn_.nodes[ind].station;
        if (ist == 0) continue;
        for (std::size_t r = 0; r < K_; ++r) {
            const double lambda = num_traits<T>::to_double(sn_.rates(ist - 1, r));
            if (std::isnan(lambda) || lambda <= 0.0) continue;
            if (sn_.service[ist - 1][r].type != lang::ProcessType::EXP)
                throw UnsupportedError(ErrorKind{"NonExponential"},
                    "solver_ssa_nrm_spn: source '" + sn_.nodes[ind].name + "' class '" +
                    sn_.classes[r].name +
                    "' has a non-exponential arrival, which the SPN reaction network cannot "
                    "express; use method='serial' or SolverJMT");
            bool found_place = false;
            for (std::size_t jnd = 0; jnd < I_; ++jnd) {
                if (sn_.nodes[jnd].nodetype != qn::NodeType::Place) continue;
                for (std::size_t s = 0; s < K_; ++s) {
                    const double p =
                        num_traits<T>::to_double(sn_.rtnodes(ind * K_ + r, jnd * K_ + s));
                    if (p <= 0.0) continue;
                    found_place = true;
                    Rx rec;
                    rec.node = ind + 1;
                    rec.mode = 0;
                    rec.S.assign(NS_, 0.0);
                    rec.S[slot(jnd, s)] += 1.0;
                    rec.base_rate = lambda * p;
                    rec.nservers = 1.0;
                    rec.dep_slots.push_back(slot(jnd, s));
                    rx_.push_back(rec);
                    producers_[ind][r].push_back(rx_.size() - 1);
                }
            }
            if (!found_place)
                throw UnsupportedError(ErrorKind{"Topology"}, "solver_ssa_nrm_spn: source '" + sn_.nodes[ind].name +
                                       "' class '" + sn_.classes[r].name +
                                       "' does not route to any Place; the SPN path needs a "
                                       "Source->Place arc");
        }
    }
 
    if (rx_.empty())
        throw InputError(
            "solver_ssa_nrm_spn: the stochastic Petri net has no timed reaction; there is "
            "nothing to simulate");
 
    // THE INITIAL MARKING, in the reference's two steps. MATLAB reads it off
    // `sn.state`, which `initDefault` fills by putting each closed class's whole
    // population at its reference station and which `Place.setState` then
    // overrides. This port has no State package -- the same position
    // `NrmEngine::build_initial_state` takes, and for the same reason -- so the
    // two steps are taken directly: a closed class whose reference station IS a
    // Place seeds that Place with its population, and a DECLARED `initmarking`
    // then replaces that Place's whole row, because an explicit marking is a
    // statement about every class of the place and not an addition to one.
    //
    // WITHOUT THE FIRST STEP a closed Petri net starts empty and deadlocks on the
    // first draw, since `initmarking` is only written by an explicit
    // `set_initial_marking` and a net that declares its tokens through a
    // ClosedClass population carries none.
    nvec0_.assign(NS_, 0.0);
    for (std::size_t r = 0; r < K_; ++r) {
        const double pop = sn_.classes[r].population;
        if (!std::isfinite(pop) || pop <= 0.0) continue;
        const std::size_t rs = sn_.classes[r].refstat;
        if (rs < 1 || rs > M_) continue;
        const std::size_t ind = sn_.station_to_node[rs - 1] - 1;
        if (sn_.nodes[ind].nodetype != qn::NodeType::Place) continue;
        nvec0_[slot(ind, r)] = pop;
    }
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (sn_.nodes[ind].nodetype != qn::NodeType::Place) continue;
        const typename std::map<std::size_t, std::vector<T>>::const_iterator im =
            sn_.initmarking.find(ind + 1);
        if (im == sn_.initmarking.end()) continue;
        for (std::size_t r = 0; r < K_; ++r) {
            const double v =
                r < im->second.size() ? num_traits<T>::to_double(im->second[r]) : 0.0;
            if (!std::isfinite(v))
                throw UnsupportedError(ErrorKind{"Place"}, "solver_ssa_nrm_spn: place '" + sn_.nodes[ind].name +
                                       "' declares an infinite initial marking, which the "
                                       "reaction network cannot count");
            nvec0_[slot(ind, r)] = v > 0.0 ? v : 0.0;
        }
    }
 
    // FINITE-CAPACITY PLACES DROP. A Place with a finite per-class (`classcap`)
    // or total (`cap`) capacity loses any arriving token that would exceed it,
    // which is the loss semantics JMT and the exact CTMC both give it: an
    // M/M/1/1 Place at rho = 0.5 holds a mean 1/3, not the unbounded 1. Without
    // the clamp the deposit accumulates past the bound.
    const double inf = std::numeric_limits<double>::infinity();
    pcap_slot_.assign(NS_, inf);
    for (std::size_t ind = 0; ind < I_; ++ind) {
        if (sn_.nodes[ind].nodetype != qn::NodeType::Place) continue;
        const std::size_t ist = sn_.nodes[ind].station;
        if (ist == 0) continue;
        std::vector<std::size_t> slots_here;
        for (std::size_t r = 0; r < K_; ++r) {
            slots_here.push_back(slot(ind, r));
            if (ist - 1 < sn_.classcap.size() && r < sn_.classcap[ist - 1].size()) {
                const double cc = sn_.classcap[ist - 1][r];
                if (std::isfinite(cc)) pcap_slot_[slot(ind, r)] = cc;
            }
        }
        if (ist - 1 < sn_.cap.size() && std::isfinite(sn_.cap[ist - 1]))
            place_total_caps_.push_back(std::make_pair(sn_.cap[ist - 1], slots_here));
    }
    has_caps_ = !place_total_caps_.empty();
    for (std::size_t j = 0; j < NS_ && !has_caps_; ++j)
        if (std::isfinite(pcap_slot_[j])) has_caps_ = true;
}
 
/**
 * `spnBuildMode`: the reaction record of transition `ind0` mode `m`.
 *
 * The enabling, firing and inhibiting matrices are (nnodes x nclasses), so a
 * nonzero entry (p, c) is an arc on class c of place p and its slot is the
 * (p, c) pair. An inhibiting entry is an arc only when it is FINITE: infinity
 * is how "this place never blocks the mode" is written, here as in the
 * reference.
 */
template <class T>
typename NrmSpnEngine<T>::Rx NrmSpnEngine<T>::build_mode(std::size_t ind0, std::size_t m) const {
    const qn::TransitionParam<T>& tp = sn_.transparam.at(ind0 + 1);
    Rx rec;
    rec.node = ind0 + 1;
    rec.mode = m + 1;
    rec.S.assign(NS_, 0.0);
 
    if (m < tp.enabling.size()) {
        const Matrix<T>& en = tp.enabling[m];
        for (std::size_t p = 0; p < en.rows() && p < I_; ++p)
            for (std::size_t c = 0; c < en.cols() && c < K_; ++c) {
                const double w = num_traits<T>::to_double(en(p, c));
                if (w == 0.0) continue;
                rec.en_slot.push_back(slot(p, c));
                rec.en_w.push_back(w);
                rec.S[slot(p, c)] -= w;
            }
    }
    if (m < tp.firing.size()) {
        const Matrix<T>& fir = tp.firing[m];
        for (std::size_t p = 0; p < fir.rows() && p < I_; ++p)
            for (std::size_t c = 0; c < fir.cols() && c < K_; ++c) {
                const double w = num_traits<T>::to_double(fir(p, c));
                if (w == 0.0) continue;
                rec.S[slot(p, c)] += w;
            }
    }
    if (m < tp.inhibiting.size()) {
        const Matrix<T>& inh = tp.inhibiting[m];
        for (std::size_t p = 0; p < inh.rows() && p < I_; ++p)
            for (std::size_t c = 0; c < inh.cols() && c < K_; ++c) {
                const double thr = num_traits<T>::to_double(inh(p, c));
                if (!std::isfinite(thr)) continue;
                rec.inh_slot.push_back(slot(p, c));
                rec.inh_thr.push_back(thr);
            }
    }
    for (std::size_t j = 0; j < NS_; ++j)
        if (rec.S[j] > 0.0) rec.dep_slots.push_back(j);
 
    // The exponential firing rate is the single-phase completion rate sum(D1).
    // A non-exponential firing is refused by `spn_nrm_supported` upstream and
    // again here, so no route reaches the run loop with a rate it invented.
    if (!(m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE)) {
        const bool one_phase = m < tp.firingphases.size() && tp.firingphases[m] == 1;
        if (!one_phase || m >= tp.firingproc.size() || tp.firingproc[m].D1.rows() == 0)
            throw UnsupportedError(ErrorKind{"NonExponential"}, "solver_ssa_nrm_spn: transition '" + sn_.nodes[ind0].name +
                                   "' mode " + std::to_string(m + 1) +
                                   " has a non-exponential firing, which the SPN reaction "
                                   "network cannot express; use method='serial'");
        const Matrix<T>& d1 = tp.firingproc[m].D1;
        double s = 0.0;
        for (std::size_t a = 0; a < d1.rows(); ++a)
            for (std::size_t b = 0; b < d1.cols(); ++b) s += num_traits<T>::to_double(d1(a, b));
        rec.base_rate = s;
    }
    rec.nservers = m < tp.nmodeservers.size() ? tp.nmodeservers[m] : 1.0;
    if (std::isinf(rec.nservers)) rec.nservers = lang::GlobalConstants::MaxInt;
    rec.weight = m < tp.fireweight.size() ? num_traits<T>::to_double(tp.fireweight[m]) : 1.0;
    rec.prio = m < tp.firingprio.size() ? tp.firingprio[m] : 1.0;
    // Fluid arcs, with each place named by its fluid index q.
    for (const auto& a : tp.fluid_rates(m)) {
        rec.flow_q.push_back(fluid_index(a.node));
        rec.flow_rate.push_back(a.value);
    }
    for (const auto& g : tp.fluid_guards(m)) {
        rec.guard_q.push_back(fluid_index(g.node));
        rec.guard_sgn.push_back(g.sgn);
        rec.guard_thr.push_back(g.thr);
    }
    for (const auto& a : tp.fluid_impulses(m)) {
        rec.imp_q.push_back(fluid_index(a.node));
        rec.imp_d.push_back(a.value);
    }
    for (const auto& d : tp.fluid_deps(m)) {
        rec.fdep_q.push_back(fluid_index(d.node));
        rec.fdep_fun.push_back(d.fun);
    }
    // marking-dependent flow multipliers, aligned with flow_q
    const std::vector<typename qn::TransitionParam<T>::FluidArc>& fr = tp.fluid_rates(m);
    rec.flow_fun.assign(fr.size(), std::function<double(const std::vector<double>&)>());
    rec.flow_read.assign(fr.size(), std::vector<std::size_t>());
    for (const auto& fd : tp.fluid_flow_deps(m))
        for (std::size_t t = 0; t < fr.size(); ++t) {
            if (fr[t].node != fd.node) continue;
            rec.flow_fun[t] = fd.fun;
            for (std::size_t q : fd.read) rec.flow_read[t].push_back(q - 1);
            break;
        }
    // the marking-dependent firing multiplier, read on the input places as the serial engine reads it
    // (`after_global_event`'s ep); only a fluid net reaches here with one, `ssa_check_firingdep` refusing the rest
    const bool timed = !(m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE);
    if (timed && m < tp.firingdep.size() && tp.firingdep[m]) {
        const std::function<T(const std::vector<T>&)> g = tp.firingdep[m];
        rec.gdep = [g](const std::vector<double>& mk) -> double {
            std::vector<T> mt(mk.size());
            for (std::size_t i = 0; i < mk.size(); ++i) mt[i] = num_traits<T>::from_double(mk[i]);
            return num_traits<T>::to_double(g(mt));
        };
        if (m < tp.enabling.size()) {
            const Matrix<T>& en = tp.enabling[m];
            for (std::size_t p = 0; p < en.rows() && p < I_; ++p)
                for (std::size_t c = 0; c < en.cols(); ++c)
                    if (num_traits<T>::to_double(en(p, c)) > 0.0) {
                        rec.gdep_read.push_back(p);
                        break;
                    }
        }
    }
    return rec;
}
 
/**
 * `spnFluidMarkMul`: the marking-dependent multipliers of marking `n`, `gk[k]` the firing multiplier g of timed
 * reaction k and `fm_[k]` the multiplier of each of its flows (empty when none carries one). Each function reads
 * the marking restricted to its declared token places, as the flat (node, class) marking, whose slot layout is
 * the engine's own and the node-major one of `state_events.h`.
 */
template <class T>
void NrmSpnEngine<T>::fluid_mark_mul(const std::vector<double>& n, std::vector<double>& gk) {
    const std::size_t nrx = rx_.size();
    gk.assign(nrx, 1.0);
    fm_.assign(nrx, std::vector<double>());
    for (std::size_t k = 0; k < nrx; ++k) {
        const Rx& rec = rx_[k];
        if (rec.gdep) {
            std::vector<double> mk(NS_, 0.0);
            for (std::size_t p : rec.gdep_read)
                for (std::size_t c = 0; c < K_; ++c) mk[slot(p, c)] = n[slot(p, c)];
            gk[k] = rec.gdep(mk);
        }
        for (std::size_t t = 0; t < rec.flow_fun.size(); ++t) {
            if (!rec.flow_fun[t]) continue;
            if (fm_[k].empty()) fm_[k].assign(rec.flow_q.size(), 1.0);
            std::vector<double> mk(NS_, 0.0);
            for (std::size_t p : rec.flow_read[t])
                for (std::size_t c = 0; c < K_; ++c) mk[slot(p, c)] = n[slot(p, c)];
            fm_[k][t] = rec.flow_fun[t](mk);
        }
    }
}
 
/** The fluid index of 1-based node `node1`; `sn_fluid_refusal` has already refused a fluid arc on a token place. */
template <class T>
std::size_t NrmSpnEngine<T>::fluid_index(std::size_t node1) const {
    if (node1 == 0 || node1 > I_ || fq_[node1 - 1] < 0)
        throw UnsupportedError(ErrorKind{"Structure"},
            "solver_ssa_nrm_spn: a fluid arc names node " + std::to_string(node1) + ", which is not a fluid place");
    return static_cast<std::size_t>(fq_[node1 - 1]);
}
 
/**
 * `applyPlaceCaps`: drop the tokens a firing pushed above a Place's bound.
 *
 * Only the just-deposited slots can overflow, so the clamp is local to them --
 * which also keeps it from disturbing a marking that was already over a bound
 * when the engine started.
 */
template <class T>
void NrmSpnEngine<T>::apply_caps(std::vector<double>& n,
                                 const std::vector<std::size_t>& deposited) const {
    for (std::size_t a = 0; a < deposited.size(); ++a) {
        const std::size_t j = deposited[a];
        if (n[j] > pcap_slot_[j]) n[j] = pcap_slot_[j];
    }
    for (std::size_t p = 0; p < place_total_caps_.size(); ++p) {
        const double tcap = place_total_caps_[p].first;
        const std::vector<std::size_t>& slots = place_total_caps_[p].second;
        double total = 0.0;
        for (std::size_t a = 0; a < slots.size(); ++a) total += n[slots[a]];
        double excess = total - tcap;
        for (std::size_t a = 0; a < deposited.size() && excess > 0.0; ++a) {
            const std::size_t j = deposited[a];
            if (std::find(slots.begin(), slots.end(), j) == slots.end()) continue;
            if (n[j] <= 0.0) continue;
            const double d = std::min(excess, n[j]);
            n[j] -= d;
            excess -= d;
        }
    }
}
 
/**
 * `spnCollapse`: vanishing-marking elimination.
 *
 * Fire enabled immediate modes until the marking is tangible, highest firing
 * priority first and ties drawn in proportion to firing weight. These firings
 * take zero time and advance no clock, so the timed race only ever resumes from
 * a tangible marking.
 */
template <class T>
void NrmSpnEngine<T>::collapse(std::vector<double>& n) {
    if (imm_.empty()) return;
    std::size_t steps = 0;
    while (true) {
        std::vector<std::size_t> enabled;
        for (std::size_t m = 0; m < imm_.size(); ++m)
            if (detail::spn_en_degree(n, imm_[m]) >= 1.0) enabled.push_back(m);
        if (enabled.empty()) return;
        double top_prio = -std::numeric_limits<double>::infinity();
        for (std::size_t i = 0; i < enabled.size(); ++i)
            top_prio = std::max(top_prio, imm_[enabled[i]].prio);
        std::vector<std::size_t> top;
        std::vector<double> w;
        for (std::size_t i = 0; i < enabled.size(); ++i)
            if (imm_[enabled[i]].prio == top_prio) {
                top.push_back(enabled[i]);
                w.push_back(imm_[enabled[i]].weight);
            }
        const std::size_t pick = top.size() == 1 ? top[0] : top[rng_.draw(w)];
        for (std::size_t j = 0; j < NS_; ++j) n[j] += imm_[pick].S[j];
        if (++steps > kMaxImmSteps)
            throw NumericError(ErrorKind{"NonConvergence"},
                "solver_ssa_nrm_spn: immediate-transition livelock -- the vanishing-marking "
                "collapse did not reach a tangible marking");
    }
}
 
/**
 * The Next-Reaction-Method run loop over the tangible markings.
 *
 * EVERY PROPENSITY IS REFRESHED after a firing rather than a dependency subset:
 * a firing plus its immediate cascade can change any place, and the reaction
 * count of a Petri net is small, so refreshing all of them removes any
 * dependency-graph blind spot at no cost worth measuring.
 *
 * A PLACE IS AN INF STATION, so its utilization IS its mean token count -- the
 * SPN convention the CTMC analyzer also reports -- and its throughput is the
 * summed firing rate of the modes consuming from it, counted once per firing.
 */
template <class T>
SsaSolution NrmSpnEngine<T>::run() {
    // A fluid net runs the piecewise-deterministic loop instead (`spnFluidRun`).
    if (!fl_node_.empty()) return run_fluid();
    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";
 
    std::vector<double> nvec = nvec0_;
    collapse(nvec);
    if (has_caps_) {
        std::vector<std::size_t> all(NS_);
        for (std::size_t j = 0; j < NS_; ++j) all[j] = j;
        apply_caps(nvec, all);
    }
 
    std::vector<double> Ak(nrx, 0.0), Pk(nrx, 0.0), Tk(nrx, 0.0), tau(nrx, 0.0);
    const double inf = std::numeric_limits<double>::infinity();
    for (std::size_t k = 0; k < nrx; ++k) {
        Ak[k] = detail::spn_propensity(nvec, rx_[k]);
        Pk[k] = -std::log(rng_.uniform());
        tau[k] = Ak[k] > 0.0 ? (Pk[k] - Tk[k]) / Ak[k] : inf;
    }
 
    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);
    // Warmup discard, as in the queueing-network NRM loop (`ssa_resolve_warmupfrac`).
    const std::size_t n_warm = ssa_warmup_firings(opt_.warmupfrac, opt_.samples);
    for (; n < opt_.samples; ++n) {
        if (n == n_warm && n_warm > 0) {
            out.QN = Matrix<double>(M_, K_, 0.0);
            out.UN = Matrix<double>(M_, K_, 0.0);
            out.TN = Matrix<double>(M_, K_, 0.0);
            total_time = 0.0;
        }
        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 = inf;
        for (std::size_t k = 0; k < nrx; ++k)
            if (tau[k] < dt) {
                dt = tau[k];
                kfire = k;
            }
        if (std::isinf(dt))
            throw NumericError(ErrorKind{"Deadlock"},
                "solver_ssa_nrm_spn: deadlock -- no transition is enabled, so the sample path "
                "cannot advance");
        total_time += dt;
 
        for (std::size_t ist = 0; ist < M_; ++ist) {
            const std::size_t ind = sn_.station_to_node[ist] - 1;
            for (std::size_t c = 0; c < K_; ++c) {
                const double tokens = nvec[slot(ind, c)];
                out.QN(ist, c) += tokens * dt;
                out.UN(ist, c) += tokens * dt;
                double depr = 0.0;
                const std::vector<std::size_t>& cons = consumers_[ind][c];
                for (std::size_t a = 0; a < cons.size(); ++a) depr += Ak[cons[a]];
                // A Source has no consuming transition: its throughput is the
                // aggregate arrival rate it injects, which is what makes the
                // reference station report the open class's arrival rate.
                const std::vector<std::size_t>& prod = producers_[ind][c];
                for (std::size_t a = 0; a < prod.size(); ++a) depr += Ak[prod[a]];
                out.TN(ist, c) += depr * dt;
            }
        }
 
        // One atomic firing, then the immediate cascade the new marking enabled.
        // A capacity-bound Place loses the tokens pushed above its bound BEFORE
        // the cascade sees the marking.
        for (std::size_t j = 0; j < NS_; ++j) nvec[j] += rx_[kfire].S[j];
        if (has_caps_) apply_caps(nvec, rx_[kfire].dep_slots);
        collapse(nvec);
 
        // Advance the Gibson and Bruck clocks with the PRE-firing propensities,
        // then refresh every propensity from the new marking.
        for (std::size_t k = 0; k < nrx; ++k) Tk[k] += Ak[k] * dt;
        for (std::size_t k = 0; k < nrx; ++k) Ak[k] = detail::spn_propensity(nvec, rx_[k]);
        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] : inf;
    }
 
    if (total_time > 0.0)
        for (std::size_t ist = 0; ist < M_; ++ist)
            for (std::size_t c = 0; c < K_; ++c) {
                out.QN(ist, c) /= total_time;
                out.UN(ist, c) /= total_time;
                out.TN(ist, c) /= total_time;
            }
    for (std::size_t c = 0; c < K_; ++c) {
        out.XN[c] = out.TN(sn_.classes[c].refstat - 1, c);
        for (std::size_t ist = 0; ist < M_; ++ist)
            out.RN(ist, c) = out.TN(ist, c) > 0.0 ? out.QN(ist, c) / out.TN(ist, c) : 0.0;
        if (out.XN[c] > 0.0) out.CN[c] = sn_.classes[c].population / out.XN[c];
    }
    out.simulated_time = total_time;
    out.samples = n;
    return out;
}
 
/**
 * `spnFluidRun`: the piecewise-deterministic loop of a fluid stochastic Petri net.
 *
 * Each pass resolves the running servers and level slopes, schedules the first guard or bound event, races
 * it against the firing epochs of the timed clocks and advances the levels linearly to whichever comes first.
 * Only a FIRING counts as a sample; a guard or bound event re-reads the guards and collapses the immediate
 * modes it enabled. The level reported is its time average, and the outflow (TN) is everything that leaves the
 * place: running outflows, cut at the empty level to the inflow, plus the part of each negative impulse that
 * was actually there to remove.
 */
template <class T>
SsaSolution NrmSpnEngine<T>::run_fluid() {
    const std::size_t nrx = rx_.size();
    const std::size_t nf = fl_node_.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";
    std::vector<bool> fluid_st(M_, false);
    for (std::size_t q = 0; q < nf; ++q) fluid_st[fl_ist_[q]] = true;
 
    std::vector<double> f = fl_f0_;
    std::vector<double> nvec = nvec0_;
    std::vector<double> QF(nf, 0.0), OF(nf, 0.0);
    {
        // the outflow of the initial collapse precedes the observation window and is discarded
        std::vector<double> of0(nf, 0.0);
        collapse_fluid(nvec, f, std::vector<double>(nf, 0.0), of0);
    }
    if (has_caps_) {
        std::vector<std::size_t> all(NS_);
        for (std::size_t j = 0; j < NS_; ++j) all[j] = j;
        apply_caps(nvec, all);
    }
    std::vector<double> Pk(nrx, 0.0), Tk(nrx, 0.0), a0(nrx, 0.0), Ik(nrx, 0.0), tau(nrx, 0.0);
    for (std::size_t k = 0; k < nrx; ++k) Pk[k] = -std::log(rng_.uniform());
 
    const double inf = std::numeric_limits<double>::infinity();
    double total_time = 0.0;
    std::size_t n = 0, zero_run = 0, ndet = 0;
    const std::size_t n_warm = ssa_warmup_firings(opt_.warmupfrac, opt_.samples);
    bool warm_pending = n_warm > 0;
    std::vector<double> eff, v, vin, vout;
    std::vector<double> gk(nrx, 1.0);
    fm_.assign(nrx, std::vector<double>());
    bool has_mark_dep = false;
    for (std::size_t k = 0; k < nrx; ++k) {
        has_mark_dep = has_mark_dep || static_cast<bool>(rx_[k].gdep);
        for (const auto& ff : rx_[k].flow_fun) has_mark_dep = has_mark_dep || static_cast<bool>(ff);
    }
    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);
    while (n < opt_.samples) {
        if (warm_pending && n == n_warm) {
            out.QN = Matrix<double>(M_, K_, 0.0);
            out.UN = Matrix<double>(M_, K_, 0.0);
            out.TN = Matrix<double>(M_, K_, 0.0);
            std::fill(QF.begin(), QF.end(), 0.0);
            std::fill(OF.begin(), OF.end(), 0.0);
            total_time = 0.0;
            warm_pending = false;
        }
        // the marking multipliers hold for the whole segment: only a firing or an immediate collapse changes the
        // marking, and either ends the segment
        if (has_mark_dep) fluid_mark_mul(nvec, gk);
        fluid_resolve(nvec, f, eff, v, vin, vout);
        for (std::size_t k = 0; k < nrx; ++k) a0[k] = rx_[k].base_rate * eff[k] * gk[k];
 
        // the next scheduled non-firing event
        double tdet = inf, fdet = 0.0;
        std::size_t qdet = 0;
        fluid_next_det(f, v, tdet, qdet, fdet);
 
        // the next firing, from each clock's remaining unit-rate time
        double dt_r = inf;
        std::size_t kfire = 0;
        for (std::size_t k = 0; k < nrx; ++k) {
            tau[k] = a0[k] > 0.0 ? detail::spn_fluid_fire_time(rx_[k], a0[k], f, v, Pk[k] - Tk[k], tdet) : inf;
            if (tau[k] < dt_r) {
                dt_r = tau[k];
                kfire = k;
            }
        }
        const bool fire = dt_r <= tdet;
        const double dt = fire ? dt_r : tdet;
        if (std::isinf(dt))
            throw NumericError(ErrorKind{"Deadlock"},
                "solver_ssa_nrm_spn: deadlock -- no transition is enabled and no fluid level is moving, so the "
                "sample path cannot advance");
 
        // tallies over the segment
        for (std::size_t k = 0; k < nrx; ++k)
            Ik[k] = a0[k] > 0.0 ? detail::spn_fluid_hazard(rx_[k], a0[k], f, v, dt) : 0.0;
        total_time += dt;
        for (std::size_t ist = 0; ist < M_; ++ist) {
            if (fluid_st[ist]) continue;
            const std::size_t ind = sn_.station_to_node[ist] - 1;
            for (std::size_t c = 0; c < K_; ++c) {
                const double tokens = nvec[slot(ind, c)];
                out.QN(ist, c) += tokens * dt;
                out.UN(ist, c) += tokens * dt;
                double fired = 0.0;
                for (std::size_t k : consumers_[ind][c]) fired += Ik[k];
                for (std::size_t k : producers_[ind][c]) fired += Ik[k];
                out.TN(ist, c) += fired;
            }
        }
        for (std::size_t q = 0; q < nf; ++q) {
            QF[q] += f[q] * dt + 0.5 * v[q] * dt * dt;
            // an empty place cut to its inflow passes on exactly what it receives
            if (v[q] == 0.0 && vin[q] < vout[q] && f[q] <= 0.0)
                OF[q] += vin[q] * dt;
            else
                OF[q] += vout[q] * dt;
        }
        for (std::size_t k = 0; k < nrx; ++k) Tk[k] += Ik[k];
        for (std::size_t q = 0; q < nf; ++q) f[q] = std::min(std::max(f[q] + v[q] * dt, 0.0), fl_cap_[q]);
 
        if (!fire) {
            f[qdet] = fdet;
            ++ndet;
            if (dt <= 0.0) {
                if (++zero_run > 1000)
                    throw NumericError(ErrorKind{"Deadlock"},
                        "solver_ssa_nrm_spn: the fluid level of place " + sn_.nodes[fl_node_[qdet]].name +
                        " keeps re-crossing a guard at the same instant: a guarded flow drives its own level back "
                        "across its own guard (a sliding mode), which a sample path has no value for");
            } else {
                zero_run = 0;
            }
            collapse_fluid(nvec, f, v, OF);
            if (ndet > 1000 * std::max<std::size_t>(opt_.samples, 1))
                throw NumericError(ErrorKind{"Deadlock"},
                    "solver_ssa_nrm_spn: more than 1000 guard or bound events per firing: the fluid levels "
                    "chatter. Check the guards against the flows they gate");
            continue;
        }
        zero_run = 0;
 
        // One atomic firing, its impulses, then the immediate cascade the new marking enabled.
        for (std::size_t j = 0; j < NS_; ++j) nvec[j] += rx_[kfire].S[j];
        if (has_caps_) apply_caps(nvec, rx_[kfire].dep_slots);
        fluid_impulse(f, rx_[kfire], OF);
        collapse_fluid(nvec, f, v, OF);
        Pk[kfire] -= std::log(rng_.uniform());
 
        ++n;
        if (n % console_every == 0)
            line::util::LineConsole::iter(
                static_cast<long>(n / console_every),
                "simulated %zu of %zu samples (%.0f%%), simulated time %.4g", n,
                static_cast<std::size_t>(opt_.samples),
                100.0 * static_cast<double>(n) / static_cast<double>(opt_.samples), total_time);
    }
 
    if (total_time > 0.0) {
        for (std::size_t ist = 0; ist < M_; ++ist)
            for (std::size_t c = 0; c < K_; ++c) {
                out.QN(ist, c) /= total_time;
                out.UN(ist, c) /= total_time;
                out.TN(ist, c) /= total_time;
            }
        for (std::size_t q = 0; q < nf; ++q) {
            QF[q] /= total_time;
            OF[q] /= total_time;
        }
    }
    // A fluid place reports its level in the first class column: QN the mean level, UN the level over the
    // capacity (the level itself when unbounded), TN the outflow.
    for (std::size_t q = 0; q < nf; ++q) {
        const std::size_t ist = fl_ist_[q];
        for (std::size_t c = 0; c < K_; ++c) out.QN(ist, c) = out.UN(ist, c) = out.TN(ist, c) = 0.0;
        if (K_ == 0) continue;
        out.QN(ist, 0) = QF[q];
        out.UN(ist, 0) = std::isfinite(fl_cap_[q]) ? QF[q] / fl_cap_[q] : QF[q];
        out.TN(ist, 0) = OF[q];
    }
    for (std::size_t c = 0; c < K_; ++c) {
        out.XN[c] = out.TN(sn_.classes[c].refstat - 1, c);
        for (std::size_t ist = 0; ist < M_; ++ist)
            out.RN(ist, c) = out.TN(ist, c) > 0.0 ? out.QN(ist, c) / out.TN(ist, c) : 0.0;
        if (out.XN[c] > 0.0) out.CN[c] = sn_.classes[c].population / out.XN[c];
    }
    out.simulated_time = total_time;
    out.samples = n;
    return out;
}
 
/**
 * `spnFluidResolve`: the running servers of every timed mode and the slope of every level, solved together.
 * A guard is read with the lookahead of the slope and the slope sums the flows of the modes whose guards hold;
 * starting from the inclusive reading the two are iterated to agreement, and a pair that never agrees within
 * 20 rounds is a sliding mode.
 */
template <class T>
void NrmSpnEngine<T>::fluid_resolve(const std::vector<double>& n, const std::vector<double>& f,
                                    std::vector<double>& eff, std::vector<double>& v, std::vector<double>& vin,
                                    std::vector<double>& vout) const {
    const std::size_t nrx = rx_.size();
    const std::vector<double> still(f.size(), 0.0);
    std::vector<double> deg(nrx, 0.0);
    std::vector<bool> gok(nrx, true), gnew(nrx, true);
    for (std::size_t k = 0; k < nrx; ++k) {
        deg[k] = std::min(detail::spn_en_degree(n, rx_[k]), rx_[k].nservers);
        gok[k] = detail::spn_fluid_guard_holds(rx_[k], f, still);
    }
    eff.assign(nrx, 0.0);
    for (int it = 0; it < 20; ++it) {
        for (std::size_t k = 0; k < nrx; ++k) eff[k] = gok[k] ? deg[k] : 0.0;
        fluid_slope(f, eff, v, vin, vout);
        for (std::size_t k = 0; k < nrx; ++k) gnew[k] = detail::spn_fluid_guard_holds(rx_[k], f, v);
        if (gnew == gok) return;
        gok = gnew;
    }
    throw NumericError(ErrorKind{"Deadlock"},
        "solver_ssa_nrm_spn: a guarded fluid flow switches its own guard on and off at the threshold (a sliding "
        "mode), which a sample path has no value for");
}
 
/**
 * `spnFluidSlope`: running inflows minus running outflows, each scaled by its marking multiplier (`fm_`), cut to
 * zero where a bound stops the level.
 */
template <class T>
void NrmSpnEngine<T>::fluid_slope(const std::vector<double>& f, const std::vector<double>& eff,
                                  std::vector<double>& v, std::vector<double>& vin,
                                  std::vector<double>& vout) const {
    const std::size_t nf = fl_node_.size();
    vin.assign(nf, 0.0);
    vout.assign(nf, 0.0);
    for (std::size_t k = 0; k < rx_.size(); ++k) {
        if (eff[k] <= 0.0) continue;
        for (std::size_t t = 0; t < rx_[k].flow_q.size(); ++t) {
            const std::size_t q = rx_[k].flow_q[t];
            double rt = rx_[k].flow_rate[t] * eff[k];
            if (k < fm_.size() && !fm_[k].empty()) rt *= fm_[k][t];
            if (rt > 0.0)
                vin[q] += rt;
            else
                vout[q] -= rt;
        }
    }
    v.assign(nf, 0.0);
    for (std::size_t q = 0; q < nf; ++q) {
        v[q] = vin[q] - vout[q];
        if ((f[q] <= 0.0 && v[q] < 0.0) || (f[q] >= fl_cap_[q] && v[q] > 0.0)) v[q] = 0.0;
    }
}
 
/** `spnFluidNextDet`: the first guard threshold or bound a moving level reaches, when, and at which level. */
template <class T>
void NrmSpnEngine<T>::fluid_next_det(const std::vector<double>& f, const std::vector<double>& v, double& tdet,
                                     std::size_t& qdet, double& fdet) const {
    tdet = std::numeric_limits<double>::infinity();
    qdet = 0;
    fdet = 0.0;
    for (std::size_t q = 0; q < fl_node_.size(); ++q) {
        if (v[q] > 0.0 && std::isfinite(fl_cap_[q])) {
            const double t = (fl_cap_[q] - f[q]) / v[q];
            if (t < tdet) {
                tdet = t;
                qdet = q;
                fdet = fl_cap_[q];
            }
        } else if (v[q] < 0.0) {
            const double t = f[q] / (-v[q]);
            if (t < tdet) {
                tdet = t;
                qdet = q;
                fdet = 0.0;
            }
        }
    }
    const std::vector<Rx>* sets[2] = {&rx_, &imm_};
    for (int s = 0; s < 2; ++s)
        for (const Rx& rec : *sets[s])
            for (std::size_t g = 0; g < rec.guard_q.size(); ++g) {
                const std::size_t q = rec.guard_q[g];
                const double thr = rec.guard_thr[g];
                if (v[q] == 0.0 || std::fabs(f[q] - thr) <= 1e-12 * std::max(1.0, std::fabs(thr))) continue;
                const double t = (thr - f[q]) / v[q];
                if (t > 0.0 && t < tdet) {
                    tdet = t;
                    qdet = q;
                    fdet = thr;
                }
            }
    tdet = std::max(tdet, 0.0);
}
 
/**
 * `spnFluidImpulse`: the impulses of a firing, each clipped at the bounds of its place; what an impulse would
 * carry past them is lost, and the part of a negative impulse that was there to remove is outflow.
 */
template <class T>
void NrmSpnEngine<T>::fluid_impulse(std::vector<double>& f, const Rx& rec, std::vector<double>& OF) const {
    for (std::size_t t = 0; t < rec.imp_q.size(); ++t) {
        const std::size_t q = rec.imp_q[t];
        const double fnew = std::min(std::max(f[q] + rec.imp_d[t], 0.0), fl_cap_[q]);
        if (fnew < f[q]) OF[q] += f[q] - fnew;
        f[q] = fnew;
    }
}
 
/**
 * `spnCollapseFluid`: `collapse` with the guards and impulses of the immediate modes. A guard is read with the
 * lookahead of the slope `v` the level had coming in.
 */
template <class T>
void NrmSpnEngine<T>::collapse_fluid(std::vector<double>& n, std::vector<double>& f, const std::vector<double>& v,
                                     std::vector<double>& OF) {
    if (imm_.empty()) return;
    std::size_t steps = 0;
    while (true) {
        std::vector<std::size_t> enabled;
        for (std::size_t m = 0; m < imm_.size(); ++m)
            if (detail::spn_en_degree(n, imm_[m]) >= 1.0 && detail::spn_fluid_guard_holds(imm_[m], f, v))
                enabled.push_back(m);
        if (enabled.empty()) return;
        double top_prio = -std::numeric_limits<double>::infinity();
        for (std::size_t i = 0; i < enabled.size(); ++i) top_prio = std::max(top_prio, imm_[enabled[i]].prio);
        std::vector<std::size_t> top;
        std::vector<double> w;
        for (std::size_t i = 0; i < enabled.size(); ++i)
            if (imm_[enabled[i]].prio == top_prio) {
                top.push_back(enabled[i]);
                w.push_back(imm_[enabled[i]].weight);
            }
        const std::size_t pick = top.size() == 1 ? top[0] : top[rng_.draw(w)];
        for (std::size_t j = 0; j < NS_; ++j) n[j] += imm_[pick].S[j];
        fluid_impulse(f, imm_[pick], OF);
        if (++steps > kMaxImmSteps)
            throw NumericError(ErrorKind{"NonConvergence"},
                "solver_ssa_nrm_spn: immediate-transition livelock -- the vanishing-marking collapse did not reach "
                "a tangible marking");
    }
}
 
}  // namespace ssa
}  // namespace line
 
#endif  // LINE_SOLVERS_SSA_SOLVER_SSA_NRM_SPN_H