ldes_engine.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_LDES_LDES_ENGINE_H
#define LINE_SOLVERS_LDES_LDES_ENGINE_H
 
/**
 * @file
 * @ingroup line_solvers
 * The NATIVE LDES discrete-event engine.
 *
 * WHAT THIS IS. `solvers/wrappers/ldes/solver_ldes.h` is a CLIENT: it writes a
 * model.json, runs the Java engine as a subprocess and parses the result back.
 * This header is the engine itself, in process, and is the C++ twin of
 * `jline/solvers/ldes/handlers/Solver_ssj.java`. Do not confuse the two
 * directories -- `wrappers/ldes/` shells out, `ldes/` simulates.
 *
 * WHAT IT COVERS TODAY. Networks over Source, Queue, Delay, Router,
 * ClassSwitch, Logger and Sink; every renewal, Markovian and trace-driven
 * service law the language carries (see `ldes_sampler.h`); the buffered
 * disciplines FCFS, LCFS, SIRO, HOL, the LCFS priority variants, SJF, LJF,
 * SEPT, LEPT, EDD, EDF, INF; the preemptive disciplines FCFSPR/PI,
 * LCFSPR/PI, their PRIO variants, SRPT, SRPTPRIO, PSJF, FB, LRPT, SETF and
 * FSP, in both resume and restart; the sharing disciplines PS, DPS, GPS,
 * their priority variants and LPS; multi-server stations; class priorities;
 * finite station and per-class capacity with DROP; balking and reneging;
 * load dependence (`sn.lldscaling`) and class dependence (`sn.cdscaling`);
 * open and closed classes; class switching; G-network NEGATIVE and CATASTROPHE
 * signals, with their batch-size law and removal policy; and confidence
 * intervals by overlapping or non-overlapping batch means.
 *
 * WHAT IT REFUSES, by name, in `ldes_engine_reject` (this list was rewritten on
 * 2026-08-01 against the code; the tranche that followed the first draft added
 * most of what the old text said was missing): the EXT discipline; a REPLY
 * signal, whose synchronous call holds a server this engine does not model; a
 * blocking rule that is neither DROP, BAS nor BBS; a continuously shared server
 * or an NHPP under SLOTTED mode, neither of which lands on the slot lattice; and any
 * node kind outside Source, Sink, Queue, Delay, Router, ClassSwitch, Logger,
 * Cache, Place and Transition. POLLING, PAS, FSP, SETF, EDD and EDF ARE
 * simulated, as are caches, Petri nets, fork-join, finite capacity regions,
 * retrial orbits, setup/delay-off, MAPt/PHt, slotted time, transient runs, busy
 * periods and replications. SERVER BREAKDOWNS are simulated as of 2026-08-15:
 * `NetworkStruct::breakdownparam` carries the two clocks and the degraded
 * `downServiceRates`, and the engine runs the outage as a state-dependent rate
 * change -- a zero multiplier stalls the job in service with its residual
 * intact, a positive one is the degraded server -- so a job is interrupted and
 * never evicted. They are refused under SLOTTED mode, where the failure and
 * repair epochs would have to land on the lattice too. Each refusal is a
 * definite feature with a definite semantics, and serving one as though it were
 * something else produces a run that finishes and reports numbers no downstream
 * test can distinguish from correct ones.
 *
 * WHAT IS TRANSCRIBED FROM THE JAVA ENGINE, because these choices decide the
 * numbers and not just the speed:
 *
 *  - THE BUDGET IS SERVICE COMPLETIONS, not simulated time and not arrivals.
 *  - WARMUP IS MSER-5 over EVENT-SPACED observations, interval
 *    `maxEvents / 1000` completions, truncation chosen on the AGGREGATE queue
 *    length with a per-series fallback (the closed-network case, whose
 *    constant total carries no transient signal).
 *  - THE ESTIMATORS ARE TIME-WEIGHTED AND TRUNCATION-DIFFERENCED: QLen and
 *    Util are (integral at the end - integral at the truncation) / elapsed,
 *    never an unweighted mean of the per-interval averages, which overweights
 *    congested epochs because they contain more events per unit time.
 *  - RESPONSE TIME IS THE FULL TALLY. `getAvgResponseTime` averages every
 *    observation including the warmup ones, so Little's law holds only up to
 *    the warmup mass. That is the reference's behaviour, mirrored rather than
 *    corrected, or the two engines would disagree on a short run.
 *  - THE SOURCE ROW OF TN IS THE NOMINAL ARRIVAL RATE, not a measured one.
 *
 * RANDOMNESS, AND HOW FAR IT NOW AGREES WITH THE JAVA ENGINE. The generators
 * are no longer Mersenne Twisters: `ldes_sampler.h` draws every renewal family
 * from `rng::Mrg32k3a`, a bit-exact reproduction of SSJ's MRG32k3a, through the
 * inverse-CDF quantiles of `ldes_ssj_variates.h`, and every Markovian family
 * from `rng::JavaRandom` through the transcribed `Map_sample.MapSampler` walk,
 * because the reference sends those two families to those two generators. What
 * is NOT yet aligned is stream OWNERSHIP: `Solver_ssj` seeds one stream per
 * (node, class) from `seed + offset`, while this engine still carries three
 * streams for the whole run (arrivals, services, routing). Until that is
 * matched a fixed seed does not reproduce the Java sample path, so a
 * cross-engine test must still compare means against a tolerance derived from
 * the standard error rather than element by element. The seed bands are kept separate
 * per role (arrivals, services, routing) for the reason the Java engine
 * separates them: one shared stream makes the routing draw the very uniform
 * that generated the interarrival time, which correlates a probabilistic split
 * with the gap length and biases the destinations in opposite directions.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <functional>
#include <limits>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <string>
#include <vector>
 
#include "line/lang/qn/network_struct.h"
#include "line/solvers/ldes/ldes_busyperiod.h"
#include "line/solvers/ldes/ldes_cache.h"
#include "line/solvers/ldes/ldes_region.h"
#include "line/solvers/ldes/ldes_spn.h"
#include "line/solvers/ldes/ldes_sampler.h"
#include "line/solvers/ldes/ldes_station.h"
#include "line/solvers/ldes/ldes_stats.h"
#include "line/solvers/wrappers/ldes/ldes_options.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace ldes {
namespace engine {
 
/** Station roles this engine simulates. */
/**
 * Station roles this engine simulates.
 *
 * A Fork and a Join appear in `sn.stations` because the refresh registers
 * them, but they SERVE NOTHING: they are handled by the routing layer, which
 * replicates at one and synchronizes at the other. Giving them a role of their
 * own keeps them out of the service-station setup, which would otherwise
 * demand a service law they do not have.
 */
enum class Role { Source, Queue, Delay, Synchronization };
 
/**
 * Event tags, ordered so that a service completion precedes an arrival
 * scheduled for the same instant.
 *
 * The order matters on any model whose rates produce exact ties (a Det
 * service, a closed class placed at t=0), and it is the completions-then-
 * arrivals phase order of the Java engine's `slotPhase`. Ties inside one tag
 * break on the insertion sequence, so the event list is a total order and the
 * path is reproducible from the seed.
 */
enum EventKind {
    EV_DEPARTURE = 0,
    EV_FIRING = 6,
    EV_SWITCHOVER = 1,
    EV_SETUP = 2,
    EV_RETRIAL = 3,
    EV_RENEGE = 4,
    EV_ARRIVAL = 5,
    EV_BREAKDOWN = 7
};
 
/** One scheduled event. */
struct Event {
    double t = 0.0;
    int kind = EV_ARRIVAL;
    std::size_t station = 0;
    std::size_t cls = 0;
    std::size_t slot = 0;   ///< server slot at a buffered station
    /**
     * Identifies the job the event was scheduled for.
     *
     * A preempted or resharing job leaves a departure event behind that the
     * event list gives no way to cancel, so every completion carries the tag
     * its job held when it was scheduled and is DISCARDED if the two no longer
     * agree. Without it a preempted job would depart twice: once on its stale
     * event and once on the one its resumption schedules.
     */
    std::uint64_t tag = 0;
    Job job;                ///< Delay and PS carry the job in the event
    std::uint64_t seq = 0;
};
 
/** Later-is-greater, so the priority queue pops the earliest event. */
struct EventLater {
    bool operator()(const Event& a, const Event& b) const {
        if (a.t != b.t) return a.t > b.t;
        if (a.kind != b.kind) return a.kind > b.kind;
        return a.seq > b.seq;
    }
};
 
/** One routing alternative out of a (node, class) pair. */
struct RouteEntry {
    std::size_t node = 0;     ///< 1-based destination NODE
    std::size_t station = 0;  ///< its station index, or M when it is not one
    std::size_t cls = 0;
    bool sink = false;
    double cum = 0.0;
};
 
/**
 * One destination NODE of a (node, class) pair, with the class it switches to.
 *
 * THE STRATEGY CHOOSES THE NODE AND THE SWITCH ROW CHOOSES THE CLASS, in that
 * order, which is `selectDestinationWithClassSwitch` and the reason the two are
 * separated at all. A strategy applied to (node, class) PAIRS would resample
 * the class as well -- uniformly under RAND, in turn under RROBIN -- and
 * destroy the class-switch probabilities the pair was carrying.
 *
 * The list is built in ASCENDING NODE INDEX because a round-robin pointer walks
 * it in order: any other order is a different, and equally arbitrary, cycle.
 */
struct RouteDest {
    std::size_t node = 0;     ///< 1-based destination NODE
    std::size_t station = 0;  ///< its station index, or M when it is not one
    bool sink = false;
    double mass = 0.0;        ///< total probability into this node, over classes
    std::vector<std::pair<std::size_t, double>> cls;  ///< (class, probability)
};
 
/**
 * Refuse, by name, every model outside this engine's scope.
 *
 * The refusal is the point. A PS station served FCFS, a preemptive-resume
 * station served non-preemptively, or a blocking model run without blocking
 * all finish and report plausible numbers; nothing downstream can tell them
 * from correct ones, so the check has to happen here, before the run.
 */
template <class T>
void ldes_engine_reject(const qn::NetworkStruct<T>& sn, const LdesOptions& o) {
    using lang::NodeType;
    using lang::SchedStrategy;
 
    for (const auto& nd : sn.nodes) {
        switch (nd.nodetype) {
            case NodeType::Source:
            case NodeType::Sink:
            case NodeType::Queue:
            case NodeType::Delay:
            case NodeType::Router:
            case NodeType::ClassSwitch:
            case NodeType::Logger:
            case NodeType::Fork:
            case NodeType::Join:
            case NodeType::Cache:
            case NodeType::Place:
            case NodeType::Transition:
                break;
            default:
                throw UnsupportedError(
                    std::string("SolverLDES (native engine): node kind ") +
                    lang::node_type_to_text(nd.nodetype) +
                    " is not ported yet; use the subprocess client");
        }
    }
    // HETEROGENEOUS SERVER POOLS are simulated here as of 2026-09-04: each job
    // occupies one server of one compatible pool and is served at that pool's
    // own law, which is what makes a pooled model a different system from the
    // multiserver it would otherwise look like. What is NOT ported is the pools
    // under a SHARING discipline, where the reference has no per-pool treatment
    // either -- `calculatePSRates` never consults the server types, so the JAR
    // silently flattens them. Rather than reproduce that, the sharing case is
    // refused by name; the max-min fair allocation it would need is the layered
    // engine's `SnCompatShare`, and neither engine carries it for a flat model.
    //
    // A FIFO depository is still understood by no engine at all.
    for (std::size_t si = 0; si < sn.stations.size(); ++si) {
        const auto& st = sn.stations[si];
        if (!st.server_types.empty()) {
            const SchedStrategy sc = st.sched;
            const bool sharing = (sc == SchedStrategy::PS || sc == SchedStrategy::DPS ||
                                  sc == SchedStrategy::GPS || sc == SchedStrategy::LPS ||
                                  sc == SchedStrategy::PSPRIO || sc == SchedStrategy::DPSPRIO ||
                                  sc == SchedStrategy::GPSPRIO);
            if (sharing)
                throw UnsupportedError(
                    std::string("SolverLDES (native engine): station '") + st.name +
                    "' declares heterogeneous server pools under a sharing discipline, which "
                    "no engine allocates per pool; use a queueing discipline, or the "
                    "subprocess client, which flattens them");
            if (sc == SchedStrategy::POLLING || sc == SchedStrategy::PAS)
                throw UnsupportedError(
                    std::string("SolverLDES (native engine): station '") + st.name +
                    "' declares heterogeneous server pools under " +
                    std::string(lang::sched_to_text(sc)) +
                    ", whose controller picks the server itself and has no pool to pick from");
        }
        for (lang::DepartureDiscipline d : st.departure_discipline)
            if (d != lang::DepartureDiscipline::NORMAL)
                throw UnsupportedError(
                    std::string("SolverLDES (native engine): station '") + st.name +
                    "' sets a non-NORMAL depository departure discipline, which no "
                    "engine implements");
    }
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        const auto& st = sn.stations[i];
        if (st.nodetype == NodeType::Source) continue;
        switch (st.sched) {
            case SchedStrategy::INF:
            case SchedStrategy::FCFS:
            case SchedStrategy::LCFS:
            case SchedStrategy::SIRO:
            case SchedStrategy::HOL:
            case SchedStrategy::LCFSPRIO:
            case SchedStrategy::SJF:
            case SchedStrategy::LJF:
            case SchedStrategy::SEPT:
            case SchedStrategy::LEPT:
            case SchedStrategy::PS:
            case SchedStrategy::DPS:
            case SchedStrategy::GPS:
            case SchedStrategy::PSPRIO:
            case SchedStrategy::DPSPRIO:
            case SchedStrategy::GPSPRIO:
            case SchedStrategy::LPS:
            case SchedStrategy::LCFSPR:
            case SchedStrategy::LCFSPI:
            case SchedStrategy::LCFSPRPRIO:
            case SchedStrategy::LCFSPIPRIO:
            case SchedStrategy::FCFSPR:
            case SchedStrategy::FCFSPI:
            case SchedStrategy::FCFSPRPRIO:
            case SchedStrategy::FCFSPIPRIO:
            case SchedStrategy::SRPT:
            case SchedStrategy::SRPTPRIO:
            case SchedStrategy::PSJF:
            case SchedStrategy::FB:
            case SchedStrategy::LRPT:
            case SchedStrategy::SETF:
            case SchedStrategy::FSP:
            case SchedStrategy::EDD:
            case SchedStrategy::EDF:
            case SchedStrategy::POLLING:
            case SchedStrategy::PAS:
            case SchedStrategy::OI:
                break;
            default:
                // EXT IS NOT A SERVICE DISCIPLINE AND IS REFUSED ON PURPOSE.
                // `Source.m` and `Sink.m` set it as their own marker, so it
                // never reaches a served station by any supported route; a
                // Queue explicitly built with it is refused by the reference
                // too, at its model layer -- `Queue.<init>` throws "Routing
                // Strategy is not supported in JLINE" before the Java engine
                // sees the model. Serving it as INF, which is how `npfqn_sqd`
                // and `sn_has_product_form` group it, would be a guess at a
                // semantics neither reference implements.
                throw UnsupportedError(
                    "SolverLDES (native engine): station '" + st.name +
                    "' uses a scheduling this engine does not simulate (EXT is the "
                    "Source and Sink marker, not a queueing discipline)");
        }
    }
    // BAS and BBS are simulated; RSRD and the retrial rules are not. Only a
    // rule on a CAPPED station can ever fire, and the refresh gives every
    // station a default rule including the Source, so an uncapped station is
    // not consulted -- rejecting on the default would refuse every model.
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        if (sn.stations[i].nodetype == NodeType::Source) continue;
        if (!std::isfinite(sn.stations[i].cap)) continue;
        for (std::size_t r = 0; r < sn.nclasses; ++r) {
            if (i >= sn.droprule.size() || r >= sn.droprule[i].size()) continue;
            const lang::DropStrategy d = sn.droprule[i][r];
            if (d != lang::DropStrategy::DROP && d != lang::DropStrategy::BAS &&
                d != lang::DropStrategy::BBS && d != lang::DropStrategy::WAITQ)
                throw UnsupportedError(
                    "SolverLDES (native engine): station '" + sn.stations[i].name +
                    "' declares a blocking rule that is not DROP, BAS or BBS");
        }
    }
    // AN ARRIVAL BATCH IS STORED BUT NOT SIMULATED HERE. `Station::arrival_batch`
    // is read by the reader and written by the writer, and BatchArrival is in
    // this engine's declared feature set, but nothing below draws a batch size:
    // every epoch releases one job. That answered a DIFFERENT model in silence
    // -- an M[3]/M/2 queue came back as the M/M/2 queue at a third of the task
    // rate -- so it is refused by name instead. The Java engine implements it,
    // and the wrappers try the runners in order, so a refusal here routes the
    // model to the jar rather than losing the feature.
    for (std::size_t i = 0; i < sn.nstations; ++i) {
        const auto& stb = sn.stations[i];
        for (std::size_t r = 0; r < stb.arrival_batch.size(); ++r)
            if (!stb.arrival_batch[r].disabled)
                throw UnsupportedError(
                    "SolverLDES (native engine): station '" + stb.name +
                    "' declares an arrival batch, which this engine does not draw yet");
    }
    // A REPLY SIGNAL completes a SYNCHRONOUS CALL: the caller keeps its server
    // while the callee runs and gets it back only when the reply arrives. That
    // is simulated -- see the `pending_reply` map below -- but only where a
    // server EXISTS to hold. A sharing station divides one server continuously
    // among everyone present and has no slot to park, so a caller there would
    // be served as though its call were asynchronous.
    for (std::size_t r = 0; r < sn.nclasses; ++r) {
        if (r >= sn.syncreply.size() || sn.syncreply[r] == 0) continue;
        for (std::size_t i = 0; i < sn.nstations; ++i) {
            if (sn.stations[i].nodetype == NodeType::Source) continue;
            if (i < sn.disabled.size() && r < sn.disabled[i].size() && sn.disabled[i][r]) continue;
            if (is_ps_family(sn.stations[i].sched) || sn.stations[i].sched == SchedStrategy::INF)
                throw UnsupportedError(
                    "SolverLDES (native engine): class '" + sn.classes[r].name +
                    "' makes a synchronous call at station '" + sn.stations[i].name +
                    "', which has no server slot to hold across it");
        }
    }
    // SLOTTED TIME EXCLUDES TWO FAMILIES, and refuses rather than approximates.
    // NHPP integrates a cumulative intensity over CONTINUOUS time, so its
    // inverse transform does not land on the lattice; the PS family divides the
    // server continuously among the jobs in service, which has no discrete-time
    // counterpart at all. Serving either on a lattice would answer a different
    // model.
    if (o.slotted) {
        for (std::size_t i = 0; i < sn.nstations; ++i) {
            if (sn.stations[i].nodetype == NodeType::Source) continue;
            if (is_ps_family(sn.stations[i].sched))
                throw UnsupportedError(
                    "SolverLDES (native engine): station '" + sn.stations[i].name +
                    "' shares its server continuously, which slotted mode has no "
                    "counterpart for");
            for (std::size_t r = 0; r < sn.nclasses; ++r)
                if (!sn.disabled[i][r] && sn.procid(i + 1, r + 1) == lang::ProcessType::NHPP)
                    throw UnsupportedError(
                        "SolverLDES (native engine): station '" + sn.stations[i].name +
                        "' carries an NHPP, whose intensity is integrated over continuous "
                        "time and does not land on the slot lattice");
        }
        for (std::size_t i = 0; i < sn.nstations; ++i)
            if (sn.breakdownparam.find(i + 1) != sn.breakdownparam.end())
                throw UnsupportedError(
                    "SolverLDES (native engine): station '" + sn.stations[i].name +
                    "' declares a server breakdown, whose failure and repair epochs are drawn on "
                    "a continuous clock and do not land on the slot lattice");
    }
}
 
/** One station's mutable state. */
struct StationState {
    Role role = Role::Queue;
    lang::SchedStrategy sched = lang::SchedStrategy::FCFS;
    bool ps = false;
    bool preemptive = false;
    bool resume = false;   ///< PR (continue the residual) as against PI (redraw)
    std::size_t nservers = 1;
    std::size_t lps_limit = 0;  ///< LPS admission cap; 0 = none
    double cap = std::numeric_limits<double>::infinity();
    std::vector<double> classcap;
    /**
     * `sn.droprule`, per class, declared on the DESTINATION station.
     *
     * JMT's convention and the reference's: the policy belongs to the station
     * whose capacity is limited, not to the upstream one whose server ends up
     * blocked. A model that declared it on the source would be read as having
     * no policy at all.
     */
    std::vector<lang::DropStrategy> droprule;
    /**
     * Jobs held OUTSIDE this station's buffer but counted in its queue length:
     * a BAS or BBS job whose destination is full occupies the upstream server
     * and is charged to the destination, which is where it is queueing for.
     * Charging it upstream would let a full station report a queue below its
     * own capacity while jobs pile up in front of it.
     */
    std::vector<double> blocked_at;
    std::vector<double> weight;      ///< schedparam, the DPS/GPS weights
    std::vector<double> class_mean;  ///< mean service per class, for SEPT/LEPT
    /**
     * `sn.lldscaling`: the RATE multiplier when n jobs are present, indexed
     * n-1. It divides the sampled service time at an ordinary station and IS
     * the effective server count at a sharing one -- the same table read two
     * ways, because a share of c servers and a c-fold rate are the same thing
     * only for PS.
     */
    std::vector<double> lld;
    /** `sn.cdscaling`: beta_r(n), a per-class RATE at the station's population. */
    lang::CdScaling<double> cd;
    bool has_cd = false;
    /**
     * A PASS-AND-SWAP station serves an ORDERED LIST whose total rate is a
     * function of the whole class sequence (Dorsman and Gardner 2024).
     *
     * The rate of the token at position p is the INCREMENT
     * mu(c1..cp) - mu(c1..c_{p-1}), which is what makes the station
     * order-independent. On a completion at position p the departing job is not
     * the one at p: the SWAP GRAPH chains forward from p to the last position
     * reachable through it, and that job departs while the chain shifts into
     * the hole. An empty swap graph is the plain OI station, where p itself
     * departs.
     */
    bool pas = false;
    std::function<double(const std::vector<std::size_t>&)> pas_rate;
    std::vector<std::vector<bool>> pas_swap;
    std::vector<Job> pas_list;
    std::uint64_t pas_tag = 0;
 
    /**
     * A POLLING SERVER visits the per-class buffers in a fixed cyclic order,
     * paying a SWITCHOVER on each leg.
     *
     * It advances ONE STEP AT A TIME and pays that step's switchover whether or
     * not the buffer it reaches holds work. Skipping ahead to the next
     * non-empty buffer for a single switchover is a different and much faster
     * discipline: it shortens the cycle and understates the waiting time by
     * tens of percent against Takagi's exact formulas, which assume the walking
     * server. The one exception is a lap made entirely of zero-time legs, which
     * would spin forever at the same instant: there the server PARKS and the
     * next arrival restarts it.
     *
     * The switchover charged to a leg is the one stored against the buffer
     * being LEFT -- that is Takagi's r_i and what `Queue.setSwitchover` means --
     * not the one against the buffer being entered.
     */
    bool polling = false;
    std::size_t poll_at = 0;         ///< the buffer the server is serving or heading to
    bool poll_switching = false;
    bool poll_parked = false;
    std::size_t poll_budget = 0;     ///< remaining visit budget; 0 = exhausted
    lang::PollingType poll_type = lang::PollingType::EXHAUSTIVE;
    std::size_t poll_k = 1;
    std::vector<Sampler> switchover;
    std::vector<bool> has_switchover;
 
    /**
     * SETUP AND DELAY-OFF: a station that powers down when it falls idle.
     *
     * The server runs a four-state machine. It starts OFF; an arrival puts it
     * in SETUP for the setup time and only then does service begin; when it
     * empties it enters DELAYOFF and shuts down when that timer expires, so a
     * job arriving during the delay-off window pays NOTHING and one arriving
     * after it pays the full setup. Collapsing the two idle states makes every
     * arrival to an idle server pay the setup, which is the model of a server
     * with no delay-off at all.
     */
    bool has_setup = false;
    Sampler setup_time, delayoff_time;
    bool setup_on = false;      ///< the server has completed its setup
    bool setup_running = false; ///< a setup is in progress
    double delayoff_at = std::numeric_limits<double>::infinity();  ///< when it shuts down
 
    /**
     * RETRIAL: a job refused admission joins an ORBIT and tries again later,
     * instead of being dropped or queueing.
     *
     * The orbit is not a buffer. Its occupants are not at the station, do not
     * occupy its capacity and are not in its queue length; they are a separate
     * population that re-offers itself at the retrial rate. Treating a retrial
     * job as queued makes the station hold more than its own capacity, which is
     * exactly what the model says it cannot.
     */
    std::vector<Sampler> retrial;
    std::vector<bool> has_retrial;
    std::vector<int> max_attempts;
    std::vector<double> orbit_size;   ///< live, per class
    std::vector<double> tot_orbit;    ///< time integral of the orbit size
    double orbit_last = 0.0;
    std::vector<double> retried, retrial_lost;
 
    /** Reneging: the abandonment timer of a WAITING job, per class. */
    std::vector<Sampler> patience;
    std::vector<bool> has_patience;
    /** Balking: the (min, max, probability) triples an arrival consults. */
    struct BalkRule {
        double min_jobs = 0.0;
        double max_jobs = -1.0;  ///< -1 = unbounded above, the wire's spelling
        double probability = 0.0;
    };
    std::vector<std::vector<BalkRule>> balk;
    std::vector<Sampler> svc;
    std::vector<bool> off;
 
    /**
     * HETEROGENEOUS SERVER POOLS, empty on a station whose servers are alike.
     *
     * A pool is `count` servers that accept only `compatible` classes and serve
     * them at a law of their own, so a model that declares them is a DIFFERENT
     * system from the multiserver it would otherwise look like: flattening the
     * pools into `nservers` answers that other system, plausibly and silently.
     *
     * The slots are laid out POOL BY POOL in declaration order, so slot
     * `type_first[t] .. type_first[t] + type_count[t] - 1` belongs to pool `t`
     * and `server_type[slot]` is the inverse. Every other part of the engine
     * addresses a server by its slot and needs to know nothing about this.
     */
    bool has_pools = false;
    std::vector<std::size_t> server_type;   ///< per slot, its pool
    std::vector<std::size_t> type_first;    ///< per pool, its first slot
    std::vector<std::size_t> type_count;    ///< per pool, its server count
    std::vector<std::vector<bool>> type_compat;   ///< [pool][class]
    std::vector<std::vector<Sampler>> type_svc;   ///< [pool][class]
    std::vector<std::vector<bool>> type_has_svc;  ///< [pool][class], a law of its own
    std::vector<std::vector<double>> type_rate;   ///< [pool][class], 1/mean, for FSF
    lang::HeteroSchedPolicy hetero_policy = lang::HeteroSchedPolicy::ORDER;
    /** ALIS/FAIRNESS rotate this order; a used pool goes to the back. */
    std::vector<std::size_t> type_order;
    /** ALFS walks the pools least-flexible first: fewest compatible classes. */
    std::vector<std::size_t> alfs_order;
 
    /**
     * The waiting room, as an EXPLICIT HEAP rather than a `std::priority_queue`.
     *
     * FSP is why. Its order is the virtual finish time of a processor-sharing
     * schedule over the residual work PRESENT AT THE STATION, so it changes
     * whenever the station's composition does and cannot be stamped on a job
     * when it arrives. A container that hides its storage cannot be reordered
     * or rescanned; this one is scanned linearly for FSP and kept as a heap
     * for every other discipline, whose keys are per-job constants.
     */
    std::vector<Job> buffer;
    WaitCmp cmp;
    std::vector<Job> server;
    std::vector<bool> server_busy;
    std::vector<double> server_start;   ///< instant the held job entered service
    std::vector<std::uint64_t> server_tag;
    /**
     * A server holding a completed job it cannot hand on, because its
     * destination is full. It is NOT free: refusing to start a new service is
     * the whole content of blocking, and a blocked slot treated as idle turns
     * a blocking network into a lossless one with the same topology.
     */
    std::vector<bool> server_blocked;
    /**
     * A server HELD ACROSS A SYNCHRONOUS CALL. Its job has left for the callee
     * and the slot is NOT free: releasing it would let the caller serve someone
     * else while its own call is outstanding, which is the asynchronous model
     * under the synchronous one's name. Separate from `server_blocked`, which
     * holds a COMPLETED job whose destination is full.
     */
    std::vector<bool> server_held;
    std::vector<std::size_t> held_cls;
    std::vector<Job> blocked_job;
    std::vector<std::size_t> blocked_dest;
    std::vector<std::size_t> blocked_dest_cls;
 
    // sharing disciplines
    std::vector<PsJob> ps_jobs;
    double ps_last_update = 0.0;
 
    /**
     * True when the station's service RATE depends on its own population, i.e.
     * it declares a load- or class-dependence table.
     *
     * Such a station cannot schedule a departure once and forget it: the rate
     * in force changes every time the population does, so the residual work
     * must be integrated forward and the completion re-timed at each change --
     * the same discipline a sharing station needs, for a different reason.
     * Applying the scaling only at service start would time every completion
     * at the rate that happened to hold when it began.
     */
    bool state_dependent = false;
    /** max(c, max(alpha)): the peak capacity the utilization is reported against. */
    double util_peak = 1.0;
    /**
     * The class-dependent speed a SHARING station is running at right now, one
     * entry per class.
     *
     * It is cached rather than recomputed because `ps_advance` integrates the
     * interval that just ENDED while the caller has already written the new
     * population: reading the handle there would credit that interval at the
     * speed that only starts now. `ps_reschedule` refreshes it in the same
     * breath as it re-times the departures, which is the moment the new speed
     * actually takes effect -- the same flush-then-install order
     * `Accum::set_busy_scale` uses for load dependence.
     */
    std::vector<double> ps_cd;
    /**
     * Server breakdown: the server alternates up and down on two clocks.
     *
     * `down_scale[r]` is the RATE MULTIPLIER while down, `down_rate(r) *
     * E[S(i,r)]`, so the whole outage is expressed as a state-dependent rate
     * change and needs no separate machinery: a zero multiplier stalls the job
     * with its residual intact (nothing is rescheduled, so nothing departs
     * until repair) and a positive one is the degraded server. That is exactly
     * the reference's semantics -- a breakdown does not evict the job in
     * service, and `downServiceRates` degrades rather than stops.
     */
    bool has_breakdown = false;
    bool up = true;
    std::vector<double> down_scale;
    /** The declared absolute rate, kept so the scale can be formed once E[S] is known. */
    std::vector<double> down_rate_raw;
    Sampler failure_time, repair_time;
    double sd_last_update = 0.0;
};
 
}  // namespace engine
 
/**
 * Simulate `sn` in process and return the same record the subprocess client
 * parses out of an `ldes-result` document.
 *
 * A metric the run did not measure stays EMPTY rather than becoming a matrix
 * of zeros: zero is a measurement and absence is not.
 */
template <class T>
LdesResult ldes_engine_solve_one(const qn::NetworkStruct<T>& sn, const LdesOptions& o) {
    using namespace engine;
    using lang::NodeType;
    using lang::SchedStrategy;
 
    ldes_engine_reject(sn, o);
 
    const std::size_t M = sn.nstations, K = sn.nclasses;
    if (M == 0 || K == 0) throw InputError("SolverLDES (native engine): empty model");
 
    const std::uint64_t max_events =
        static_cast<std::uint64_t>(o.events > 0 ? o.events : o.samples);
    if (max_events == 0)
        throw InputError("SolverLDES (native engine): the completion budget is zero");
 
    // ---- random streams, one band per role ---------------------------------
    const std::uint64_t base =
        (o.seed >= 0) ? static_cast<std::uint64_t>(o.seed) : std::random_device{}();
    /**
     * Snap a sampled duration onto the slot lattice, or REFUSE it.
     *
     * A non-lattice sample is an ERROR, not something to round: rounding
     * silently changes the distribution, and a model whose service law is not
     * lattice-valued is not a discrete-time model. The reference refuses for
     * the same reason. Only Geometric, Det with an integral slot count, and
     * DMAP produce lattice values.
     */
    const bool slotted = o.slotted;
    const double slot_len = o.slot_length;
    auto slot_snap = [&](double v, const char* what) -> double {
        if (!slotted || v == 0.0) return v;
        const double slots = v / slot_len;
        const double rounded = std::floor(slots + 0.5);
        if (rounded < 1.0 || std::fabs(slots - rounded) > 1e-9 * std::max(1.0, slots))
            throw InputError(std::string("SolverLDES (native engine): slotted mode sampled a ") +
                             what + " of " + std::to_string(v) +
                             ", which is not a positive multiple of the slot length " +
                             std::to_string(slot_len) +
                             "; a discrete-time model needs lattice-valued interarrival and "
                             "service times, e.g. Geometric or Det on an integral slot count");
        return rounded * slot_len;
    };
 
    /**
     * ONE STREAM PER (node, class), with the reference's own offsets.
     *
     * `Solver_ssj` never shares a stream between two generators: it builds one
     * MRG32k3a and one java.util.Random per (node, class) and seeds them from
     * `seed + offset`, with
     *
     *     arrival (source s, class k)   stream (s*K + k)*10,  random the same + 2000
     *     service (svc j, class k)      both ((numSources + j)*K + k)*10 + 1000
     *     routing                       ROUTING_SEED_BASE = 900000
     *
     * where `numSources` counts SOURCE nodes and `j` indexes the SERVICE nodes
     * in order, neither of which is the station index. Sharing three run-wide
     * streams instead, as this engine did, changes the interleaving of every
     * draw: two stations that each consume their own stream in the reference
     * would take turns from one here, so the paths diverge from the first event
     * even with identical generators and identical quantiles.
     *
     * The station-level draws that are not per-class -- a polling switchover,
     * a setup or delay-off, an SPN transition -- sit on `g_aux[station]` in a
     * band of their own, mirroring the reference's own
     * `(numSources + numServiceNodes + j)` family without claiming to match it
     * draw for draw.
     */
    const long long seed_ll = static_cast<long long>(base);
    const long long Kll = static_cast<long long>(K);
    std::size_t num_sources = 0, num_service_nodes = 0;
    std::vector<std::size_t> svc_index(M, 0);
    for (std::size_t i = 0; i < M; ++i) {
        if (sn.stations[i].nodetype == NodeType::Source) {
            ++num_sources;
        } else {
            svc_index[i] = num_service_nodes++;
        }
    }
    const long long nsrc = static_cast<long long>(num_sources);
    const long long nsvc = static_cast<long long>(num_service_nodes);
 
    std::vector<Rng> g_arr;   ///< [class], at the (single) source
    g_arr.reserve(K);
    for (std::size_t k = 0; k < K; ++k) {
        const long long off = (static_cast<long long>(k)) * 10;  // srcIdx == 0
        g_arr.push_back(Rng(seed_ll, off, off + 2000));
    }
    std::vector<std::vector<Rng>> g_svc;  ///< [station][class]
    g_svc.reserve(M);
    for (std::size_t i = 0; i < M; ++i) {
        std::vector<Rng> row;
        row.reserve(K);
        for (std::size_t k = 0; k < K; ++k) {
            const long long off =
                ((nsrc + static_cast<long long>(svc_index[i])) * Kll + static_cast<long long>(k)) *
                    10 + 1000;
            row.push_back(Rng(seed_ll, off));
        }
        g_svc.push_back(row);
    }
    // [station][pool][class]: a heterogeneous pool's own service stream, banded
    // past `g_svc` so a pooled station's draws never interleave with the
    // class-level ones and a homogeneous model's sequence is untouched.
    std::vector<std::vector<std::vector<Rng>>> g_hsvc(M);
    for (std::size_t i = 0; i < M; ++i) {
        const std::size_t nT = sn.stations[i].server_types.size();
        g_hsvc[i].resize(nT);
        for (std::size_t t = 0; t < nT; ++t) {
            g_hsvc[i][t].reserve(K);
            for (std::size_t k = 0; k < K; ++k) {
                const long long off =
                    ((nsrc + static_cast<long long>(svc_index[i])) * Kll +
                     static_cast<long long>(k)) * 10 +
                    600000 + static_cast<long long>(t) * 137;
                g_hsvc[i][t].push_back(Rng(seed_ll, off));
            }
        }
    }
    std::vector<Rng> g_aux;  ///< [station], the draws that carry no class
    g_aux.reserve(M);
    for (std::size_t i = 0; i < M; ++i) {
        const long long off = ((nsrc + nsvc + static_cast<long long>(svc_index[i])) * Kll) * 10;
        g_aux.push_back(Rng(seed_ll, off));
    }
    Rng g_routing(seed_ll, 900000);  // ROUTING_SEED_BASE
    // The SPN transition clocks are a model-level draw, not a station's.
    // The reference bands the place generators at +5000; this is that band.
    Rng g_spn(seed_ll, 5000);
    // The fork degree and branch activation draws, on the reference's own
    // `new Random(seed + 88888)` band. DRAWN FROM ONLY BY A FORK THAT DECLARES
    // AN OVERRIDE: a plain fork must consume nothing here, or every seeded
    // golden and the bit-exact agreement with the Java engine would move.
    Rng g_fork(seed_ll, 88888);
 
    /**
     * True when a fork's parameters actually vary: some branch is uncertain,
     * some link carries a count other than the scalar, or some link draws its
     * count. A fork that declares the matrices but leaves them uniform stays on
     * the fixed path, so it keeps drawing nothing.
     */
    auto fork_is_variable = [&](std::size_t node) -> bool {
        const qn::ForkParam<double>* fp = sn.fork_param_of(node);
        if (fp == 0) return false;
        const double tpl = sn.nodes[node - 1].tasks_per_link;
        for (std::size_t k = 0; k < fp->fan_out_link.rows(); ++k)
            for (std::size_t r = 0; r < fp->fan_out_link.cols(); ++r) {
                const double p = fp->fan_out_prob(k, r);
                if (p == 0.0) continue;  // link not taken
                if (p != 1.0) return true;
                if (fp->fan_out_link(k, r) != tpl) return true;
                if (!fp->fan_out_dist[k][r].disabled) return true;
            }
        return false;
    };
 
    /** One draw from a jobs-per-link distribution, by inverse CDF on its pmf. */
    auto sample_fork_degree = [&](const lang::Distrib<double>& d) -> int {
        double tot = 0.0;
        for (std::size_t e = 0; e < d.params.size(); ++e) tot += d.params[e];
        const double u = g_fork.aux.next_double() * tot;
        double accp = 0.0;
        for (std::size_t e = 0; e < d.params.size(); ++e) {
            accp += d.params[e];
            if (u <= accp)
                return static_cast<int>((d.trace.empty() ? static_cast<double>(e + 1) : d.trace[e]) +
                                        0.5);
        }
        return static_cast<int>(
            (d.trace.empty() ? static_cast<double>(d.params.size()) : d.trace.back()) + 0.5);
    };
 
    // ---- stations ----------------------------------------------------------
    std::vector<StationState> S(M);
    std::size_t source_st = M;
    std::vector<int> classprio(K, 0);
    std::vector<double> classdeadline(K, std::numeric_limits<double>::infinity());
    for (std::size_t r = 0; r < K; ++r) {
        classprio[r] = sn.classes[r].prio;
        classdeadline[r] = sn.classes[r].deadline;
    }
 
    for (std::size_t i = 0; i < M; ++i) {
        const auto& st = sn.stations[i];
        StationState& s = S[i];
        s.sched = st.sched;
        s.off.assign(K, true);
        s.class_mean.assign(K, 0.0);
        s.weight.assign(K, 1.0);
        for (std::size_t r = 0; r < K; ++r)
            if (r < st.schedparam.size()) s.weight[r] = num_traits<T>::to_double(st.schedparam[r]);
 
        if (st.nodetype == NodeType::Source) {
            s.role = Role::Source;
            source_st = i;
        } else if (st.nodetype == NodeType::Fork || st.nodetype == NodeType::Join ||
                   st.nodetype == NodeType::Place || st.nodetype == NodeType::Transition) {
            // A Place holds tokens and a Transition moves them; neither serves,
            // so neither belongs in the service-station setup.
            s.role = Role::Synchronization;
            // A JOIN STILL REPORTS. It serves nothing, so it stays out of the
            // setup below, but it holds siblings and releases parents and the
            // result assembly must read those rows; leaving every class OFF is
            // what emitted it as an all-zero station. `off` gates the report
            // here, not admission: `admit` is never reached through a Join.
            if (st.nodetype == NodeType::Join) s.off.assign(K, false);
        } else if (st.nodetype == NodeType::Delay || st.sched == SchedStrategy::INF) {
            s.role = Role::Delay;
        } else {
            s.role = Role::Queue;
            s.ps = is_ps_family(st.sched);
            s.preemptive = is_preemptive(st.sched);
            s.resume = is_preemptive_resume(st.sched);
            const double c = st.nservers;
            if (!(c >= 1.0) || !std::isfinite(c))
                throw InputError("SolverLDES (native engine): station '" + st.name +
                                 "' has a server count that is neither finite nor at least one");
            s.nservers = static_cast<std::size_t>(c + 0.5);
            s.cap = st.cap;
            s.classcap.assign(K, std::numeric_limits<double>::infinity());
            s.blocked_at.assign(K, 0.0);
            s.droprule.assign(K, lang::DropStrategy::DROP);
            for (std::size_t r = 0; r < K; ++r)
                if (i < sn.droprule.size() && r < sn.droprule[i].size())
                    s.droprule[r] = sn.droprule[i][r];
            for (std::size_t r = 0; r < K; ++r)
                if (i < sn.classcap.size() && r < sn.classcap[i].size())
                    s.classcap[r] = sn.classcap[i][r];
            if (st.sched == SchedStrategy::LPS) s.lps_limit = s.nservers;
        }
 
        // One entry per class, neutral until a sharing station with a class
        // dependence refreshes it in `ps_reschedule`.
        s.ps_cd.assign(K, 1.0);
        for (const T& v : st.lldscaling) s.lld.push_back(num_traits<T>::to_double(v));
        // The peak capacity that normalizes this station's utilization is
        // max(c, max(alpha)) -- CTMC's own ceff, which is what makes the two
        // report the same number on a load-dependent station. Without load
        // dependence it is just c, so nothing else moves.
        s.util_peak = static_cast<double>(s.nservers);
        for (double a : s.lld) s.util_peak = std::max(s.util_peak, a);
        // CLASS DEPENDENCE beta_r(n) AND JOINT DEPENDENCE eta_i(n) ENTER THE
        // SAMPLE-PATH RATE THE SAME WAY, and a station declaring both is scaled
        // by their PRODUCT -- `initializeClassDependence` and `State.java` both
        // fold them that way. They stay separate on the struct because they
        // carry different modelling claims (beta keeps the product form, eta
        // does not), but a simulator has no use for that distinction: it needs
        // the rate in force, which is the product.
        if (st.cdscaling || st.jdscaling) {
            s.has_cd = true;
            // THE DECLARED PEAK IS THE NORMALIZER, not a peak swept out of the
            // handle. `solver_ssa_serial` and `solver_ssa_nrm_space` take the
            // PRODUCT of whichever declared peaks are present, and drop the
            // server count while doing so, because the rate in force is that
            // same product; this engine must answer with their number. A handle
            // with no declared peak is a model defect here exactly as it is
            // there -- an open class has no lattice to sweep, so a swept peak
            // would be an invented bound.
            double dep_peak = 1.0;
            const std::vector<T>* pks[2] = {&st.cdscalingpeak, &st.jdscalingpeak};
            const bool on[2] = {static_cast<bool>(st.cdscaling), static_cast<bool>(st.jdscaling)};
            const char* names[2] = {"setClassDependence", "setJointDependence"};
            for (std::size_t h = 0; h < 2; ++h) {
                if (!on[h]) continue;
                if (pks[h]->empty())
                    throw InputError(
                        std::string("SolverLDES: station '") + st.name +
                        "' declares a dependent scaling with no declared peak rate. Utilization "
                        "there is T*E[S]/peak, so pass the peak to " + names[h]);
                double m = 0.0;
                for (const T& v : *pks[h]) m = std::max(m, num_traits<T>::to_double(v));
                if (!(m > 0.0))
                    throw InputError(std::string("SolverLDES: station '") + st.name +
                                     "' declares a non-positive peak rate for " + names[h]);
                dep_peak *= m;
            }
            s.util_peak = dep_peak;
            const auto& beta = st.cdscaling;
            const auto& eta = st.jdscaling;
            s.cd = [beta, eta](const std::vector<double>& n) {
                std::vector<T> nt(n.size());
                for (std::size_t k = 0; k < n.size(); ++k) nt[k] = num_traits<T>::from_double(n[k]);
                std::vector<double> bd, ed;
                if (beta) {
                    const std::vector<T> b = beta(nt);
                    bd.resize(b.size());
                    for (std::size_t k = 0; k < b.size(); ++k) bd[k] = num_traits<T>::to_double(b[k]);
                }
                if (eta) {
                    const std::vector<T> e = eta(nt);
                    ed.resize(e.size());
                    for (std::size_t k = 0; k < e.size(); ++k) ed[k] = num_traits<T>::to_double(e[k]);
                }
                if (bd.empty()) return ed;
                if (ed.empty()) return bd;
                // Elementwise, with scalar broadcasting: either factor may be a
                // single value shared across the classes.
                const std::size_t n2 = std::max(bd.size(), ed.size());
                std::vector<double> out(n2, 1.0);
                for (std::size_t k = 0; k < n2; ++k)
                    out[k] = bd[bd.size() > 1 ? k : 0] * ed[ed.size() > 1 ? k : 0];
                return out;
            };
        }
        if (st.sched == SchedStrategy::PAS || st.sched == SchedStrategy::OI) {
            s.pas = true;
            auto pp = sn.pasparam.find(i + 1);
            if (pp != sn.pasparam.end()) {
                if (pp->second.svc_rate_fun) {
                    const auto f = pp->second.svc_rate_fun;
                    // mu(c) TAKES 1-BASED CLASS TAGS -- it is the state
                    // encoding's own convention, the one `state.h` and
                    // `after_event_station_pas` hand it and the one
                    // `oi_rate_from_json` decodes (it drops any entry below 1).
                    // This engine's `Job::cls` is 0-based, so shift on the way
                    // in: unshifted, class 0 is invisible to mu, contributes no
                    // rate increment at any position, and its jobs pile up at
                    // the station and never complete.
                    s.pas_rate = [f](const std::vector<std::size_t>& seq) {
                        std::vector<std::size_t> tags(seq.size());
                        for (std::size_t k = 0; k < seq.size(); ++k) tags[k] = seq[k] + 1;
                        return num_traits<T>::to_double(f(tags));
                    };
                }
                s.pas_swap = pp->second.swap_graph;
            }
            if (!s.pas_rate)
                throw InputError("SolverLDES (native engine): station '" + st.name +
                                 "' is a pass-and-swap station with no service rate function");
        }
        if (st.sched == SchedStrategy::POLLING) {
            s.polling = true;
            const auto pp = sn.effective_polling(i + 1);
            s.poll_type = pp.ptype;
            s.poll_k = (pp.pk >= 1) ? pp.pk : 1;
            s.poll_budget = (s.poll_type == lang::PollingType::EXHAUSTIVE)
                                ? std::numeric_limits<std::size_t>::max()
                                : ((s.poll_type == lang::PollingType::KLIMITED) ? s.poll_k : 1);
            s.poll_parked = true;
            s.switchover.resize(K);
            s.has_switchover.assign(K, false);
            for (std::size_t r = 0; r < K && r < pp.switchover.size(); ++r)
                if (!pp.switchover[r].disabled &&
                    num_traits<T>::to_double(pp.switchover[r].mean) > 0.0) {
                    s.switchover[r] = Sampler(pp.switchover[r],
                                              "the switchover of station '" + st.name +
                                                  "', class '" + sn.classes[r].name + "'");
                    s.has_switchover[r] = true;
                }
        }
        s.retrial.resize(K);
        s.has_retrial.assign(K, false);
        s.max_attempts.assign(K, 0);
        s.orbit_size.assign(K, 0.0);
        s.tot_orbit.assign(K, 0.0);
        s.retried.assign(K, 0.0);
        s.retrial_lost.assign(K, 0.0);
        {
            auto rp = sn.retrialparam.find(i + 1);
            if (rp != sn.retrialparam.end())
                for (std::size_t r = 0; r < K && r < rp->second.retrial_proc.size(); ++r)
                    if (!rp->second.retrial_proc[r].disabled) {
                        s.retrial[r] = Sampler(rp->second.retrial_proc[r],
                                               "the retrial process of station '" + st.name +
                                                   "', class '" + sn.classes[r].name + "'");
                        s.has_retrial[r] = true;
                        if (r < rp->second.max_attempts.size())
                            s.max_attempts[r] = rp->second.max_attempts[r];
                    }
        }
        s.patience.resize(K);
        s.has_patience.assign(K, false);
        s.balk.assign(K, std::vector<StationState::BalkRule>());
        for (std::size_t r = 0; r < K; ++r) {
            if (r < st.impatience.size() && st.impatience[r] == lang::ImpatienceType::RENEGING &&
                r < st.patience.size() && !st.patience[r].disabled) {
                s.patience[r] = Sampler(st.patience[r], "the patience of station '" + st.name +
                                                            "', class '" + sn.classes[r].name + "'");
                s.has_patience[r] = true;
            }
            if (r < st.balking.size() &&
                st.balking[r].strategy != lang::BalkingStrategy::NONE) {
                if (st.balking[r].strategy != lang::BalkingStrategy::QUEUE_LENGTH)
                    throw UnsupportedError(
                        "SolverLDES (native engine): station '" + st.name +
                        "' declares a balking rule that is not QUEUE_LENGTH; the "
                        "expected-wait rules are not ported yet");
                for (const auto& th : st.balking[r].thresholds) {
                    StationState::BalkRule br;
                    br.min_jobs = th.min_jobs;
                    br.max_jobs = th.max_jobs;
                    br.probability = num_traits<T>::to_double(th.probability);
                    s.balk[r].push_back(br);
                }
            }
        }
 
        {
            auto sp = sn.setupparam.find(i + 1);
            if (sp != sn.setupparam.end()) {
                lang::Distrib<T> su, doff;
                if (sp->second.last(su, doff) && !su.disabled) {
                    s.has_setup = true;
                    s.setup_time = Sampler(su, "the setup time of station '" + st.name + "'");
                    // A ZERO delay-off is not an error, it is the absence of the
                    // window: the server shuts down the instant it empties, so
                    // every busy period starts cold. Only a positive one gets a
                    // sampler; a disabled or zero entry leaves it unset.
                    if (!doff.disabled && num_traits<T>::to_double(doff.mean) > 0.0)
                        s.delayoff_time =
                            Sampler(doff, "the delay-off time of station '" + st.name + "'");
                    else
                        s.delayoff_time = Sampler();
                }
            }
        }
        {
            typename std::map<std::size_t, qn::BreakdownParam<T> >::const_iterator bp =
                sn.breakdownparam.find(i + 1);
            if (bp != sn.breakdownparam.end()) {
                s.has_breakdown = true;
                s.up = true;
                s.failure_time =
                    Sampler(bp->second.failure, "the failure time of station '" + st.name + "'");
                s.repair_time =
                    Sampler(bp->second.repair, "the repair time of station '" + st.name + "'");
                s.down_scale.assign(K, 0.0);
                s.down_rate_raw.assign(K, 0.0);
                for (std::size_t r = 0; r < K && r < bp->second.down_service_rates.size(); ++r)
                    s.down_rate_raw[r] = num_traits<T>::to_double(bp->second.down_service_rates[r]);
            }
        }
        // A BREAKDOWN IS A STATE-DEPENDENT RATE, so it reuses sd_advance /
        // sd_reschedule rather than growing a second freeze-and-retime path.
        s.state_dependent = (!s.lld.empty() || s.has_cd || s.has_breakdown);
 
        s.svc.resize(K);
        if (s.role == Role::Synchronization) continue;
        for (std::size_t r = 0; r < K; ++r) {
            s.off[r] = sn.disabled[i][r];
            if (s.off[r]) continue;
            const std::string where = "station '" + st.name + "', class '" + sn.classes[r].name + "'";
            s.svc[r] = Sampler(sn.service[i][r], where);
            s.class_mean[r] = s.svc[r].mean();
            if (!(s.class_mean[r] > 0.0))
                throw InputError("SolverLDES (native engine): " + where +
                                 " has a non-positive mean service time");
        }
 
        // ---- heterogeneous server pools ---------------------------------------
        //
        // The pools REPLACE the flat server bank rather than sitting beside it:
        // `nservers` becomes their total, so every slot loop below is already
        // over the right range and only the CHOICE of slot has to learn about
        // compatibility. A pool that declares no law of its own for a class
        // serves it at the station's class-level default, which is what the
        // reference's `generateHeteroServiceTime` falls back to.
        if (!st.server_types.empty()) {
            const std::size_t nT = st.server_types.size();
            s.has_pools = true;
            s.hetero_policy = st.hetero_policy;
            s.type_count.assign(nT, 0);
            s.type_first.assign(nT, 0);
            s.type_compat.assign(nT, std::vector<bool>(K, true));
            s.type_svc.assign(nT, std::vector<Sampler>(K));
            s.type_has_svc.assign(nT, std::vector<bool>(K, false));
            s.type_rate.assign(nT, std::vector<double>(K, 0.0));
            std::size_t total = 0;
            for (std::size_t t = 0; t < nT; ++t) {
                const typename qn::Station<T>::ServerType& pt = st.server_types[t];
                const double c = pt.count;
                if (!(c >= 1.0))
                    throw InputError("SolverLDES (native engine): station '" + st.name +
                                     "' declares server pool '" + pt.name +
                                     "' with fewer than one server");
                s.type_first[t] = total;
                s.type_count[t] = static_cast<std::size_t>(c + 0.5);
                total += s.type_count[t];
                for (std::size_t r = 0; r < K; ++r) {
                    // An empty compatibility row means every class, as the
                    // struct documents; a short one is padded the same way.
                    s.type_compat[t][r] = pt.compatible.empty() || r >= pt.compatible.size()
                                              ? true
                                              : pt.compatible[r];
                    if (s.off[r]) s.type_compat[t][r] = false;
                    if (r < pt.service.size() && !pt.service[r].disabled) {
                        const std::string where = "station '" + st.name + "', pool '" + pt.name +
                                                  "', class '" + sn.classes[r].name + "'";
                        s.type_svc[t][r] = Sampler(pt.service[r], where);
                        s.type_has_svc[t][r] = true;
                        const double mu = s.type_svc[t][r].mean();
                        if (!(mu > 0.0))
                            throw InputError("SolverLDES (native engine): " + where +
                                             " has a non-positive mean service time");
                        s.type_rate[t][r] = 1.0 / mu;
                    } else if (!s.off[r]) {
                        s.type_rate[t][r] =
                            (s.class_mean[r] > 0.0) ? 1.0 / s.class_mean[r] : 0.0;
                    }
                }
            }
            for (std::size_t r = 0; r < K; ++r) {
                if (s.off[r]) continue;
                bool served = false;
                for (std::size_t t = 0; t < nT && !served; ++t) served = s.type_compat[t][r];
                if (!served)
                    throw InputError("SolverLDES (native engine): station '" + st.name +
                                     "' declares no server pool compatible with class '" +
                                     sn.classes[r].name +
                                     "', so a job of that class would wait forever");
            }
            // The pools ARE the server bank. A declared `nservers` that
            // disagrees is the flattened view of the same pools, and taking it
            // instead would serve more or fewer jobs at once than declared.
            s.nservers = total;
            s.util_peak = static_cast<double>(total);
            if (s.sched == lang::SchedStrategy::LPS) s.lps_limit = total;
            s.server_type.assign(total, 0);
            for (std::size_t t = 0; t < nT; ++t)
                for (std::size_t j = 0; j < s.type_count[t]; ++j)
                    s.server_type[s.type_first[t] + j] = t;
            s.type_order.resize(nT);
            for (std::size_t t = 0; t < nT; ++t) s.type_order[t] = t;
            // ALFS: fewest compatible classes first, ties by declaration order,
            // so the least flexible pool is spent before the more flexible one.
            s.alfs_order = s.type_order;
            std::stable_sort(s.alfs_order.begin(), s.alfs_order.end(),
                             [&](std::size_t a, std::size_t b) {
                                 std::size_t ca = 0, cb = 0;
                                 for (std::size_t r = 0; r < K; ++r) {
                                     if (s.type_compat[a][r]) ++ca;
                                     if (s.type_compat[b][r]) ++cb;
                                 }
                                 return ca < cb;
                             });
        }
        // THE DEGRADED RATE BECOMES A MULTIPLIER, and it can only be formed
        // here, once the mean service time of each class is known: the wire
        // carries an ABSOLUTE rate while the engine's state-dependent path
        // takes a factor on the nominal one.
        if (s.has_breakdown) {
            s.down_scale.assign(K, 0.0);
            for (std::size_t r = 0; r < K; ++r) {
                if (s.off[r] || !(s.class_mean[r] > 0.0)) continue;
                if (r < s.down_rate_raw.size() && s.down_rate_raw[r] > 0.0)
                    s.down_scale[r] = s.down_rate_raw[r] * s.class_mean[r];
            }
        }
        if (!s.ps) {
            s.cmp.sched = s.sched;
            s.cmp.class_mean = &s.class_mean;
            if (s.role == Role::Queue) {
                s.server.assign(s.nservers, Job());
                s.server_busy.assign(s.nservers, false);
                s.server_start.assign(s.nservers, 0.0);
                s.server_tag.assign(s.nservers, 0);
                s.server_blocked.assign(s.nservers, false);
                s.server_held.assign(s.nservers, false);
                s.held_cls.assign(s.nservers, 0);
                s.blocked_job.assign(s.nservers, Job());
                s.blocked_dest.assign(s.nservers, 0);
                s.blocked_dest_cls.assign(s.nservers, 0);
            }
        }
    }
 
    // ---- routing, resolved once per (station, class) -----------------------
    // A Router, a ClassSwitch and a Logger are not stations. `Peff` is a
    // NODE-level matrix and keeps them, so they appear as destinations here and
    // `deliver` walks through them in one instant; a destination that is none of
    // those and none of Source/Queue/Delay/Sink/Fork/Join/Cache/Place/Transition
    // is a node this engine must not guess at.
    const std::size_t nnodes = sn.nof_nodes();
    std::vector<std::size_t> node_to_station(nnodes + 1, M);
    for (std::size_t i = 0; i < M; ++i) node_to_station[sn.station_to_node[i]] = i;
 
    /**
     * ROUTING IS RESOLVED OVER NODES, not over stations.
     *
     * A Router, a ClassSwitch and a Logger appear here as destinations: the
     * stochastic complement that removes them lives in `rt`, not in the
     * NODE-level `Peff` this table is read off. They hold nothing, so `deliver`
     * walks straight through them. A FORK and a JOIN are different again: they
     * are stateful, they change the NUMBER of jobs in flight, and no stochastic
     * complement can express that, so `deliver` replicates at a Fork and
     * synchronizes at a Join.
     */
    std::vector<std::vector<std::vector<RouteDest>>> nroute(
        nnodes + 1, std::vector<std::vector<RouteDest>>(K));
    /**
     * The dispatcher each (node, class) declares, and its round-robin pointer.
     *
     * The pointer is per (node, class) and NOT per destination: a round-robin
     * dispatcher hands out its outgoing links in turn, and one counter per link
     * would let two classes at the same node walk the cycle independently and
     * both start at the head.
     */
    std::vector<std::vector<lang::RoutingStrategy>> node_routing(
        nnodes + 1, std::vector<lang::RoutingStrategy>(K, lang::RoutingStrategy::PROB));
    std::vector<std::vector<std::size_t>> rr_counter(nnodes + 1, std::vector<std::size_t>(K, 0));
    for (std::size_t inode = 1; inode <= nnodes; ++inode)
        for (std::size_t r = 0; r < K; ++r)
            if (sn.nodes[inode - 1].routing.size() > r)
                node_routing[inode][r] = sn.nodes[inode - 1].routing[r];
    for (std::size_t inode = 1; inode <= nnodes; ++inode) {
        // THE SINK HAS NO OUTGOING PATH. The refresh writes a Sink -> Source
        // closure into the effective routing so the visit ratios of an open
        // model have a cycle to solve on; it is bookkeeping, not a route a job
        // takes. Reading it as one sends every departing job back to the
        // Source. The station-level table never saw it because the Sink is not
        // a station.
        if (sn.nodes[inode - 1].nodetype == NodeType::Sink) continue;
        // A PETRI NET DOES NOT ROUTE. Its arcs declare the structure, but the
        // tokens move by the enabling and firing matrices of the transitions,
        // so walking the routing here would move them twice.
        if (sn.nodes[inode - 1].nodetype == NodeType::Place ||
            sn.nodes[inode - 1].nodetype == NodeType::Transition)
            continue;
        for (std::size_t r = 0; r < K; ++r) {
            // ASCENDING NODE INDEX in the outer loop: this is the order the
            // round-robin pointer walks, and it is the reference's order.
            for (std::size_t j = 1; j <= nnodes; ++j) {
                RouteDest d;
                d.node = j;
                d.sink = (sn.nodes[j - 1].nodetype == NodeType::Sink);
                d.station = node_to_station[j];
                for (std::size_t s2 = 0; s2 < K; ++s2) {
                    const double p = num_traits<T>::to_double(sn.route_eff(r + 1, s2 + 1, inode, j));
                    if (!(p > 0.0)) continue;
                    d.mass += p;
                    d.cls.push_back(std::make_pair(s2, p));
                }
                if (d.cls.empty()) continue;
                // A Fork, a Join and a Cache are STATEFUL and hold jobs; a
                // ClassSwitch, a Router and a Logger do not and are walked
                // through in one step by `deliver`. Neither kind is folded away
                // by the refresh -- `Peff` is a NODE-level matrix and keeps
                // every intermediate hop -- so both must be reachable here.
                // Anything else is a node this engine would have to guess at.
                if (!d.sink && d.station >= M && sn.nodes[j - 1].nodetype != NodeType::Fork &&
                    sn.nodes[j - 1].nodetype != NodeType::Join &&
                    sn.nodes[j - 1].nodetype != NodeType::Cache &&
                    sn.nodes[j - 1].nodetype != NodeType::ClassSwitch &&
                    sn.nodes[j - 1].nodetype != NodeType::Router &&
                    sn.nodes[j - 1].nodetype != NodeType::Logger &&
                    sn.nodes[j - 1].nodetype != NodeType::Place &&
                    sn.nodes[j - 1].nodetype != NodeType::Transition)
                    throw UnsupportedError(
                        "SolverLDES (native engine): the routing crosses node '" +
                        sn.nodes[j - 1].name +
                        "', which the refresh did not fold into the station-to-station "
                        "routing");
                if (d.station < M && S[d.station].role == Role::Source)
                    throw InputError("SolverLDES (native engine): the routing sends a "
                                     "job back into the Source");
                nroute[inode][r].push_back(d);
            }
        }
    }
 
    // ---- fork and join --------------------------------------------------------
    // A Fork replicates the arriving job onto every outgoing edge; a Join holds
    // the siblings until the strategy is satisfied and releases ONE job. The
    // pair is matched by `sn.fj`, so a Fork with no Join is a model whose
    // siblings never recombine and is refused rather than run.
    std::vector<int> join_of_fork(nnodes + 1, -1);
    std::vector<bool> is_fork(nnodes + 1, false), is_join(nnodes + 1, false);
    for (const auto& fjp : sn.fj) {
        if (fjp.first > nnodes || fjp.second > nnodes) continue;
        is_fork[fjp.first] = true;
        is_join[fjp.second] = true;
        join_of_fork[fjp.first] = static_cast<int>(fjp.second);
    }
    // Whether a class routes into a Join ANYWHERE: it decides whether a phase-2
    // continuation of that class inherits the trigger's fork identity.
    std::vector<bool> spawn_joins_at(K, false);
    for (std::size_t inode = 1; inode <= nnodes; ++inode)
        for (std::size_t r = 0; r < K; ++r)
            for (const RouteDest& d : nroute[inode][r])
                if (d.node >= 1 && d.node <= nnodes &&
                    sn.nodes[d.node - 1].nodetype == NodeType::Join)
                    for (std::size_t c = 0; c < d.cls.size(); ++c)
                        spawn_joins_at[d.cls[c].first] = true;
    for (std::size_t nd = 1; nd <= nnodes; ++nd) {
        if (sn.nodes[nd - 1].nodetype == NodeType::Fork && !is_fork[nd])
            throw InputError("SolverLDES (native engine): Fork node '" + sn.nodes[nd - 1].name +
                             "' has no matching Join");
        if (sn.nodes[nd - 1].nodetype == NodeType::Join && !is_join[nd])
            throw InputError("SolverLDES (native engine): Join node '" + sn.nodes[nd - 1].name +
                             "' closes no Fork");
    }
 
    /**
     * One forked parent awaiting its siblings.
     *
     * `required` is fixed at FORK time, not at join time: the quorum is a
     * fraction of the siblings THIS fork produced, and a Join shared by two
     * forks of different fan-out would otherwise apply one fork's count to the
     * other's siblings.
     */
    struct ForkSync {
        std::size_t cls = 0;
        double t_sys = 0.0;
        int total = 0;
        int required = 0;
        /**
         * (class, instant) of every sibling this parent has parked at the Join.
         *
         * Held PER PARENT, not per Join: the Join's queue length must fall by
         * the siblings THIS synchronization releases, and a node-wide tally
         * would clear the siblings of every other parent still waiting. It also
         * dates each sibling's own wait, which is `forkedJobJoinArrivalTimes`.
         */
        std::vector<std::pair<std::size_t, double>> siblings;
    };
    std::map<std::uint64_t, ForkSync> fork_sync;
    std::uint64_t next_parent = 0;
 
 
    // ---- cache nodes ----------------------------------------------------------
    // A Cache is a ROUTING node with state: it reads an item, decides hit or
    // miss on its real list contents, and switches the job onto the hit or the
    // miss class. That makes the hit an exact sample-path event rather than a
    // rate, which is the reason to ask a simulator about a cache at all.
    std::map<std::size_t, CacheState> caches;
    for (const auto& np : sn.nodeparam) {
        const std::size_t nd = np.first;
        if (nd == 0 || nd > nnodes) continue;
        if (sn.nodes[nd - 1].nodetype != NodeType::Cache) continue;
        const auto& cp = np.second;
        CacheState cs;
        cs.node = nd;
        cs.nitems = cp.nitems;
        cs.policy = cp.replacestrat;
        for (int c : cp.itemcap)
            if (c > 0) cs.capacity.push_back(static_cast<std::size_t>(c));
        if (cs.capacity.empty()) cs.capacity.push_back(1);
        cs.lists.assign(cs.capacity.size(), std::list<std::size_t>());
        for (int v : cp.itemsize) cs.item_size.push_back(v);
        for (int v : cp.costcap) cs.cost_cap.push_back(v);
        cs.hits.assign(K, 0.0);
        cs.misses.assign(K, 0.0);
        cs.hit_class.assign(K, -1);
        cs.miss_class.assign(K, -1);
        for (std::size_t r = 0; r < K && r < cp.hitclass.size(); ++r)
            if (cp.hitclass[r] > 0) cs.hit_class[r] = static_cast<int>(cp.hitclass[r] - 1);
        for (std::size_t r = 0; r < K && r < cp.missclass.size(); ++r)
            if (cp.missclass[r] > 0) cs.miss_class[r] = static_cast<int>(cp.missclass[r] - 1);
        // The popularity is stored as a pmf per class and used as a CUMULATIVE
        // law here, so one uniform draw picks an item.
        cs.popularity.assign(K, std::vector<double>());
        for (std::size_t r = 0; r < K && r < cp.pread.size(); ++r) {
            double acc = 0.0;
            for (const T& v : cp.pread[r]) {
                acc += num_traits<T>::to_double(v);
                cs.popularity[r].push_back(acc);
            }
            if (!cs.popularity[r].empty()) cs.popularity[r].back() = 1.0;
        }
        // Delayed-hit retrieval: one class per (item, read class) carries the
        // fetch through the retrieval stations and back here.
        if (cp.retrieval_capacity > 0 && !cp.retrieval_classes.empty()) {
            cs.has_retrieval = true;
            cs.retrieval_class = cp.retrieval_classes;
            cs.retrieval_class_to_item.assign(K, -1);
            for (std::size_t it = 0; it < cs.retrieval_class.size(); ++it)
                for (std::size_t r = 0; r < cs.retrieval_class[it].size(); ++r) {
                    const std::size_t rc = cs.retrieval_class[it][r];
                    if (rc > 0 && rc <= K) cs.retrieval_class_to_item[rc - 1] = static_cast<int>(it);
                }
            cs.in_flight.assign(cs.nitems, 0);
            cs.fetch_start.assign(cs.nitems, 0.0);
            cs.held.assign(cs.nitems, std::vector<CacheState::HeldRequest>());
            cs.delayed.assign(K, 0.0);
        }
        caches[nd] = cs;
    }
 
 
    // ---- Petri net: places and transitions --------------------------------------
    // A Place is a STATION holding tokens per class; a Transition is a node that
    // moves them. The marking is the places' token counts, and it is the only
    // state the net has.
    std::vector<std::size_t> place_of(nnodes + 1, M);
    std::vector<std::size_t> place_nodes;
    for (std::size_t nd = 1; nd <= nnodes; ++nd)
        if (sn.nodes[nd - 1].nodetype == NodeType::Place) {
            place_of[nd] = place_nodes.size();
            place_nodes.push_back(nd);
        }
    std::vector<std::vector<double>> marking(place_nodes.size(), std::vector<double>(K, 0.0));
    for (std::size_t p = 0; p < place_nodes.size(); ++p) {
        auto im = sn.initmarking.find(place_nodes[p]);
        if (im != sn.initmarking.end())
            for (std::size_t r = 0; r < K && r < im->second.size(); ++r)
                marking[p][r] = num_traits<T>::to_double(im->second[r]);
    }
 
    std::vector<SpnTransition> transitions;
    for (const auto& tp : sn.transparam) {
        const std::size_t nd = tp.first;
        if (nd == 0 || nd > nnodes) continue;
        SpnTransition tr;
        tr.node = nd;
        tr.places = place_nodes;
        for (std::size_t m = 0; m < tp.second.nmodes; ++m) {
            SpnMode mode;
            // (place slot, class), flattened p*K + r, the layout `flat_marking`
            // produces and `spn_enabled`/`spn_fire` read.
            for (std::size_t p = 0; p < place_nodes.size(); ++p) {
                const std::size_t pn = place_nodes[p];
                for (std::size_t r = 0; r < K; ++r) {
                    mode.enabling.push_back(
                        (m < tp.second.enabling.size() && pn - 1 < tp.second.enabling[m].rows() &&
                         r < tp.second.enabling[m].cols())
                            ? num_traits<T>::to_double(tp.second.enabling[m](pn - 1, r))
                            : 0.0);
                    mode.inhibiting.push_back(
                        (m < tp.second.inhibiting.size() &&
                         pn - 1 < tp.second.inhibiting[m].rows() &&
                         r < tp.second.inhibiting[m].cols())
                            ? num_traits<T>::to_double(tp.second.inhibiting[m](pn - 1, r))
                            : std::numeric_limits<double>::infinity());
                    mode.firing.push_back(
                        (m < tp.second.firing.size() && pn - 1 < tp.second.firing[m].rows() &&
                         r < tp.second.firing[m].cols())
                            ? num_traits<T>::to_double(tp.second.firing[m](pn - 1, r))
                            : 0.0);
                }
            }
            if (m < tp.second.nmodeservers.size()) mode.servers = tp.second.nmodeservers[m];
            if (m < tp.second.firingprio.size()) mode.priority = tp.second.firingprio[m];
            if (m < tp.second.fireweight.size())
                mode.weight = num_traits<T>::to_double(tp.second.fireweight[m]);
            mode.immediate = (m < tp.second.timing.size() &&
                              tp.second.timing[m] == lang::TimingStrategy::IMMEDIATE);
            if (m < tp.second.firingproc.size() && !tp.second.firingproc[m].disabled) {
                const double mean = num_traits<T>::to_double(tp.second.firingproc[m].mean);
                if (mean > 0.0) mode.rate = 1.0 / mean;
            }
            if (m < tp.second.firingdep.size() && tp.second.firingdep[m]) {
                const auto f = tp.second.firingdep[m];
                const std::vector<std::size_t> pl = place_nodes;
                const std::size_t nn = nnodes, KK = K;
                // g(marking) IS NODE-INDEXED AND CLASS-SUMMED, the vector
                // `state_events.h` builds and the reader's slot lattice is keyed
                // by. Handing it the (place, class) vector instead read the
                // wrong slot whenever a place was not the node of the same
                // ordinal, and every slot once the class dimension existed.
                mode.dep = [f, pl, nn, KK](const std::vector<double>& tok) {
                    std::vector<T> mk(nn, num_traits<T>::from_int(0));
                    for (std::size_t p = 0; p < pl.size(); ++p) {
                        double s = 0.0;
                        for (std::size_t r = 0; r < KK; ++r)
                            if (p * KK + r < tok.size()) s += tok[p * KK + r];
                        if (pl[p] >= 1 && pl[p] - 1 < nn)
                            mk[pl[p] - 1] = num_traits<T>::from_double(s);
                    }
                    return num_traits<T>::to_double(f(mk));
                };
            }
            tr.modes.push_back(mode);
        }
        tr.fired.assign(tr.modes.size(), 0.0);
        transitions.push_back(tr);
    }
 
    // ---- statistics --------------------------------------------------------
    Accum acc(M, K);
    for (std::size_t i = 0; i < M; ++i) acc.util_peak[i] = S[i].util_peak;
    // Declared here and populated below: the admission lambda tracks into it,
    // and a lambda cannot capture a name introduced after it.
    BusyPeriods bp;
    std::vector<double> sys_resp_sum(K, 0.0), sys_resp_cnt(K, 0.0), sys_completed(K, 0.0);
    // Per (station, class): a model can drop at one station and balk at
    // another, and a per-class total would not say which.
    std::vector<std::vector<double>> dropped(M, std::vector<double>(K, 0.0));
    std::vector<std::vector<double>> balked(M, std::vector<double>(K, 0.0));
    std::vector<std::vector<double>> reneged(M, std::vector<double>(K, 0.0));
    std::vector<std::vector<double>> blocked_count(M, std::vector<double>(K, 0.0));
    /**
     * Every per-visit response time, when either flag asked for them.
     *
     * --trajectory counts as asking: the reference NESTS the same samples inside
     * its `transient` block, and MATLAB's `sample()` and `sampleSys()` read them
     * from there. Recording only under --respt-samples left those two returning
     * nothing on the transient path.
     */
    const bool want_respt = o.export_respt || o.export_trajectory;
    std::vector<std::vector<std::vector<double>>> resp_samples(
        want_respt ? M : 0, std::vector<std::vector<double>>(K));
    std::uint64_t job_id = 0;
    double now = 0.0;
 
    /**
     * G-NETWORK REMOVAL SIGNALS, resolved once per class.
     *
     * A removal signal is a class like any other while it is in transit and
     * stops being one the instant it reaches a station: it never joins, it
     * removes jobs already there and it is annihilated. `is_removal_signal`
     * is the whole gate, and everything the removal costs -- the registry of
     * jobs held at a Delay station, the extra draws -- is paid only where it
     * is true, so a model without signals runs exactly as before.
     */
    std::vector<bool> is_removal_signal(K, false);
    bool has_removal_signal = false;
    for (std::size_t r = 0; r < K; ++r) {
        if (r >= sn.issignal.size() || !sn.issignal[r]) continue;
        const lang::SignalType kind =
            (r < sn.signaltype.size()) ? sn.signaltype[r] : lang::SignalType::NEGATIVE;
        if (kind == lang::SignalType::NEGATIVE || kind == lang::SignalType::CATASTROPHE) {
            is_removal_signal[r] = true;
            has_removal_signal = true;
        }
    }
    /**
     * The jobs a Delay station currently holds, in arrival order.
     *
     * A Delay keeps no buffer and no server: its occupants live only in the
     * departure events already scheduled for them, and the event list has no
     * handle to cancel one. The registry is that handle -- a job removed by a
     * signal is erased from it, and the departure that arrives later finds no
     * entry and is discarded, which is the same tag trick the buffered
     * stations use against a preempted departure.
     */
    std::vector<std::vector<Job>> delay_live(has_removal_signal ? M : 0);
 
    /**
     * SYNCHRONOUS CALLS, `sn.syncreply`: the class a caller expects its reply in.
     *
     * A caller that expects one KEEPS ITS SERVER when it leaves for the callee
     * and gets it back only when the matching REPLY class returns, which is the
     * whole difference between a synchronous call and an ordinary hop. The
     * pending calls are keyed by an identity that rides on the job, because the
     * reply comes back from a station the caller never names.
     */
    std::vector<std::size_t> sync_reply(K, K);
    std::vector<bool> is_reply_signal(K, false);
    bool has_sync_call = false;
    for (std::size_t r = 0; r < K; ++r) {
        if (r < sn.syncreply.size() && sn.syncreply[r] >= 1 && sn.syncreply[r] <= K) {
            sync_reply[r] = sn.syncreply[r] - 1;
            has_sync_call = true;
        }
        if (r < sn.issignal.size() && sn.issignal[r] && r < sn.signaltype.size() &&
            sn.signaltype[r] == lang::SignalType::REPLY)
            is_reply_signal[r] = true;
    }
    /** One outstanding call: the slot it holds and the class that holds it. */
    struct PendingCall {
        std::size_t station = 0;
        std::size_t slot = 0;
        std::size_t cls = 0;
        double since = 0.0;
    };
    std::map<std::uint64_t, PendingCall> pending_reply;
    std::uint64_t next_call = 0;
 
    /** Whether any (station, class) declares immediate feedback at all. */
    bool has_immfeed = false;
    for (std::size_t i = 0; i < M && i < sn.immfeed.size(); ++i)
        for (std::size_t r = 0; r < sn.immfeed[i].size(); ++r)
            if (sn.immfeed[i][r]) has_immfeed = true;
 
    /**
     * SPAWN ON COMPLETION, `sn.classspawn`: completing a job of class r injects
     * a FRESH job of the mapped class at the same station.
     *
     * It is an LQN phase-2 continuation: the work the server owes after the
     * reply has gone back. The continuation is a new passage, so it starts its
     * own response-time clock and is not the completing job under another name.
     */
    std::vector<std::size_t> spawn_of(K, K);
    bool has_spawn = false;
    for (std::size_t r = 0; r < K; ++r)
        if (sn.classes[r].spawn >= 1 && sn.classes[r].spawn <= K) {
            spawn_of[r] = sn.classes[r].spawn - 1;
            has_spawn = true;
        }
 
    // ---- event list --------------------------------------------------------
    std::priority_queue<Event, std::vector<Event>, EventLater> evq;
    std::uint64_t seq = 0, ps_tag = 0;
 
    auto push = [&](Event e) {
        e.seq = seq++;
        evq.push(e);
    };
 
 
    /**
     * Push onto, and pop from, a station's waiting room.
     *
     * FSP takes the linear branch: its key is recomputed against the residual
     * work present at the station RIGHT NOW, over the waiting jobs and the
     * ones in service alike, so no stamped ordering can stand in for it. Every
     * other discipline has a per-job constant key and rides the heap.
     */
    auto buffer_push = [&](std::size_t i, const Job& job) {
        StationState& s = S[i];
        s.buffer.push_back(job);
        if (s.sched != SchedStrategy::FSP)
            std::push_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
    };
    auto buffer_pop = [&](std::size_t i) -> Job {
        StationState& s = S[i];
        if (s.sched != SchedStrategy::FSP) {
            std::pop_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
            Job job = s.buffer.back();
            s.buffer.pop_back();
            return job;
        }
        std::vector<double> works;
        for (const Job& j : s.buffer) works.push_back(j.remaining);
        for (std::size_t sl = 0; sl < s.nservers; ++sl)
            if (s.server_busy[sl])
                works.push_back(std::max(0.0, s.server[sl].remaining -
                                                  (now - s.server_start[sl])));
        std::size_t best = 0;
        double best_vft = std::numeric_limits<double>::infinity();
        for (std::size_t j = 0; j < s.buffer.size(); ++j) {
            const double v = fsp_virtual_finish(works, s.buffer[j].remaining,
                                                static_cast<double>(s.nservers), now);
            if (v < best_vft) {
                best_vft = v;
                best = j;
            }
        }
        Job job = s.buffer[best];
        s.buffer.erase(s.buffer.begin() + static_cast<std::ptrdiff_t>(best));
        return job;
    };
 
    // ---- heterogeneous server pools: slot choice and buffer matching ---------
    //
    // Two questions the flat server bank never had to ask. WHICH SLOT can take
    // this job, and WHICH WAITING JOB can this slot take. Everything else in the
    // engine addresses a server by its slot and is untouched.
    //
    // On a station with no pools both reduce to what was there before -- the
    // first free slot, and the head of the waiting room -- so the pool-free path
    // keeps its exact draw sequence and a seeded homogeneous run is unchanged.
 
    /** Free, unblocked, unheld, and its pool accepts `cls`. */
    auto slot_ok_for = [&](const StationState& s, std::size_t sl, std::size_t cls) -> bool {
        if (s.server_busy[sl] || s.server_blocked[sl] || s.server_held[sl]) return false;
        if (!s.has_pools) return true;
        return s.type_compat[s.server_type[sl]][cls];
    };
 
    /**
     * The free slot a job of class `cls` may enter, or `s.nservers` for none.
     *
     * The POOL is chosen first, by the station's heterogeneous scheduling
     * policy, and the slot is then the first free one inside it. That is the
     * reference's order and it matters: picking the globally-first free slot
     * would make every policy behave as ORDER.
     */
    auto free_slot_for = [&](std::size_t i, std::size_t cls) -> std::size_t {
        StationState& s = S[i];
        if (!s.has_pools) {
            for (std::size_t sl = 0; sl < s.nservers; ++sl)
                if (slot_ok_for(s, sl, cls)) return sl;
            return s.nservers;
        }
        // The pools that accept the class AND still hold a free slot.
        std::vector<std::size_t> cand;
        for (std::size_t t = 0; t < s.type_count.size(); ++t) {
            if (!s.type_compat[t][cls]) continue;
            for (std::size_t j = 0; j < s.type_count[t]; ++j)
                if (slot_ok_for(s, s.type_first[t] + j, cls)) {
                    cand.push_back(t);
                    break;
                }
        }
        if (cand.empty()) return s.nservers;
        std::size_t chosen = cand[0];
        if (cand.size() > 1) {
            switch (s.hetero_policy) {
                case lang::HeteroSchedPolicy::ALIS:
                case lang::HeteroSchedPolicy::FAIRNESS: {
                    // Round robin over the pools: the one just used goes to the
                    // back, so a steady stream spreads over the compatible
                    // pools instead of filling the first one.
                    for (std::size_t k = 0; k < s.type_order.size(); ++k) {
                        const std::size_t t = s.type_order[k];
                        if (std::find(cand.begin(), cand.end(), t) == cand.end()) continue;
                        chosen = t;
                        s.type_order.erase(s.type_order.begin() +
                                           static_cast<std::ptrdiff_t>(k));
                        s.type_order.push_back(t);
                        break;
                    }
                    break;
                }
                case lang::HeteroSchedPolicy::ALFS: {
                    // Least flexible first: spend the pool that fewest classes
                    // can use before the one that many can.
                    for (std::size_t k = 0; k < s.alfs_order.size(); ++k)
                        if (std::find(cand.begin(), cand.end(), s.alfs_order[k]) != cand.end()) {
                            chosen = s.alfs_order[k];
                            break;
                        }
                    break;
                }
                case lang::HeteroSchedPolicy::FSF: {
                    // Fastest server first, by this class's rate at the pool.
                    double best = -1.0;
                    for (std::size_t k = 0; k < cand.size(); ++k) {
                        const double rt = s.type_rate[cand[k]][cls];
                        if (rt > best) {
                            best = rt;
                            chosen = cand[k];
                        }
                    }
                    break;
                }
                case lang::HeteroSchedPolicy::RAIS: {
                    // A uniform draw over the candidate pools, off the ROUTING
                    // stream as the reference draws it.
                    std::size_t at = static_cast<std::size_t>(
                        uniform01(g_routing) * static_cast<double>(cand.size()));
                    if (at >= cand.size()) at = cand.size() - 1;
                    chosen = cand[at];
                    break;
                }
                case lang::HeteroSchedPolicy::ORDER:
                default: {
                    // Declaration order, which `cand` is already built in.
                    chosen = cand[0];
                    break;
                }
            }
        }
        for (std::size_t j = 0; j < s.type_count[chosen]; ++j) {
            const std::size_t sl = s.type_first[chosen] + j;
            if (slot_ok_for(s, sl, cls)) return sl;
        }
        return s.nservers;
    };
 
    /** Is any waiting job servable by slot `sl`? */
    auto buffer_has_for_slot = [&](std::size_t i, std::size_t sl) -> bool {
        StationState& s = S[i];
        if (s.buffer.empty()) return false;
        if (!s.has_pools) return true;
        const std::vector<bool>& ok = s.type_compat[s.server_type[sl]];
        for (std::size_t j = 0; j < s.buffer.size(); ++j)
            if (ok[s.buffer[j].cls]) return true;
        return false;
    };
 
    /**
     * Pop the job slot `sl` should serve next.
     *
     * Off a pooled station this IS `buffer_pop`. On one it is the job the
     * station's own order would have taken FIRST AMONG THOSE THE SLOT CAN
     * SERVE, which is what a pool means: a server does not skip the queue, it
     * skips the jobs it is not compatible with. The heap is drained and rebuilt
     * because the discipline's order lives in the comparator, not in the
     * sequence, so scanning the vector would answer in heap order rather than
     * in service order.
     */
    auto buffer_pop_for_slot = [&](std::size_t i, std::size_t sl) -> Job {
        StationState& s = S[i];
        if (!s.has_pools) return buffer_pop(i);
        const std::vector<bool>& ok = s.type_compat[s.server_type[sl]];
        std::vector<Job> skipped;
        Job job;
        bool got = false;
        while (!s.buffer.empty()) {
            Job cand = buffer_pop(i);
            if (ok[cand.cls]) {
                job = cand;
                got = true;
                break;
            }
            skipped.push_back(cand);
        }
        for (std::size_t j = 0; j < skipped.size(); ++j) buffer_push(i, skipped[j]);
        if (!got)
            throw InputError("SolverLDES (native engine): buffer_pop_for_slot was asked for a "
                             "job no pool of this slot can serve; guard with buffer_has_for_slot");
        return job;
    };
 
 
    /**
     * The RATE multiplier a station applies to a class at its current
     * population: load dependence times class dependence.
     *
     * A sampled requirement is DIVIDED by this, so a multiplier above one is a
     * faster station. `lldscaling` is indexed by the total population minus
     * one, saturating at the table's end, which is the reference's convention
     * and not a guard: the table is declared up to the station's capacity and a
     * closed model can hold exactly that many.
     *
     * A SHARING station does NOT take the load-dependent factor here. Its
     * `lldscaling` is read instead as the effective SERVER COUNT in
     * `ps_shares`, because c servers and a c-fold rate are the same thing only
     * when the capacity is shared -- applying both would scale it twice.
     */
    // The load-dependent speed alone, which is both a factor of the service
    // rate below and the scale the busy integral is measured in. Sharing one
    // definition is what keeps `U = T*E[S]/peak` exact on the sample path.
    auto lld_factor = [&](std::size_t i) -> double {
        StationState& s = S[i];
        if (s.lld.empty() || s.ps) return 1.0;
        double total = 0.0;
        for (std::size_t k = 0; k < K; ++k) total += acc.qlen[i][k];
        if (!(total > 0.0)) return 1.0;
        const std::size_t idx = std::min(static_cast<std::size_t>(total) - 1, s.lld.size() - 1);
        return s.lld[idx] > 0.0 ? s.lld[idx] : 1.0;
    };
    // The class-dependent speed alone. It is a factor of the RATE, so every
    // path that turns a population into a departure time has to apply it -- the
    // sharing disciplines included, which is what `ps_advance`/`ps_reschedule`
    // below missed: they built their shares from `ps_shares` and the breakdown
    // clock and never read the handle, so a PS station carrying a `cdscaling`
    // ran its whole sample path at the UNSCALED rates. On a closed
    // Delay+Queue(PS), N=3, beta(n)=min(n,2), that was Tput 0.935 against the
    // exact 1.412, and utilization inherited the error through T*E[S]/peak.
    auto cd_factor = [&](std::size_t i, std::size_t r) -> double {
        StationState& s = S[i];
        if (!s.has_cd) return 1.0;
        std::vector<double> nvec(K, 0.0);
        for (std::size_t k = 0; k < K; ++k) nvec[k] = acc.qlen[i][k];
        const std::vector<double> beta = s.cd(nvec);
        const double b = (beta.size() > 1) ? beta[r] : (beta.empty() ? 1.0 : beta[0]);
        return (b > 1e-10) ? b : 1.0;
    };
    auto rate_scaling = [&](std::size_t i, std::size_t r) -> double {
        StationState& s = S[i];
        double factor = lld_factor(i) * cd_factor(i, r);
        // A DOWN SERVER IS A ZERO FACTOR, not a skipped station: the job in
        // service keeps its residual and simply stops accruing, which is what
        // makes the outage an interruption rather than an eviction. A positive
        // `down_scale` is the degraded server of `downServiceRates`.
        if (s.has_breakdown && !s.up)
            factor *= (r < s.down_scale.size() ? s.down_scale[r] : 0.0);
        return factor;
    };
 
    /**
     * The effective server count of a SHARING station, which is where its
     * load-dependence table is read. Transcribes `getEffectivePSServerCount`:
     * with no jobs present the count is 1, not the table's first entry, so the
     * very first arrival gets a whole server.
     */
    auto ps_servers = [&](std::size_t i) -> double {
        StationState& s = S[i];
        if (s.lld.empty()) return static_cast<double>(s.nservers);
        double total = 0.0;
        for (std::size_t k = 0; k < K; ++k) total += acc.qlen[i][k];
        if (!(total > 0.0)) return 1.0;
        const std::size_t idx = std::min(static_cast<std::size_t>(total) - 1, s.lld.size() - 1);
        return s.lld[idx];
    };
 
    // ---- the sharing disciplines -------------------------------------------
    // Every population change at a PS station invalidates every departure it
    // has scheduled, so the pair below is always called together: integrate the
    // work delivered since the last change, then reschedule everything. A
    // stale event is not cancelled (the event list has no handle); it is
    // ignored on arrival by comparing its tag against the job's current one.
    auto ps_advance = [&](std::size_t i) {
        StationState& s = S[i];
        const double dt = now - s.ps_last_update;
        s.ps_last_update = now;
        // THE BUSY TIME OF A SHARING STATION IS THE INTEGRAL OF THE SHARES, not
        // of the head count. Every job at a PS station is "in service", so
        // counting heads makes the busy time n*dt and reports a utilization of
        // n/c -- above one as soon as the station queues at all. The reference
        // accumulates `rate * elapsed` per job for exactly this reason, which
        // also makes U = T*E[S] hold automatically, since the shares sum to the
        // work actually delivered.
        for (std::size_t r = 0; r < K; ++r) acc.last_busy[i][r] = now;
        if (!(dt > 0.0) || s.ps_jobs.empty()) return;
        std::vector<double> rates = ps_shares(s.sched, s.ps_jobs, ps_servers(i), s.weight, K);
        if (s.has_cd)
            for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
                rates[j] *= s.ps_cd[s.ps_jobs[j].cls];
        if (s.has_breakdown && !s.up)
            for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
                rates[j] *= (s.ps_jobs[j].cls < s.down_scale.size() ? s.down_scale[s.ps_jobs[j].cls]
                                                                   : 0.0);
        for (std::size_t j = 0; j < s.ps_jobs.size(); ++j)
            if (rates[j] > 0.0) {
                acc.tot_busy[i][s.ps_jobs[j].cls] += rates[j] * dt;
                s.ps_jobs[j].remaining = std::max(0.0, s.ps_jobs[j].remaining - rates[j] * dt);
            }
    };
    auto ps_reschedule = [&](std::size_t i) {
        StationState& s = S[i];
        // The new speed takes effect HERE, on the population the caller has just
        // written, and the departures re-timed below are what it applies to.
        if (s.has_cd)
            for (std::size_t r = 0; r < K; ++r) s.ps_cd[r] = cd_factor(i, r);
        if (s.ps_jobs.empty()) return;
        std::vector<double> rates = ps_shares(s.sched, s.ps_jobs, ps_servers(i), s.weight, K);
        if (s.has_cd)
            for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
                rates[j] *= s.ps_cd[s.ps_jobs[j].cls];
        if (s.has_breakdown && !s.up)
            for (std::size_t j = 0; j < rates.size() && j < s.ps_jobs.size(); ++j)
                rates[j] *= (s.ps_jobs[j].cls < s.down_scale.size() ? s.down_scale[s.ps_jobs[j].cls]
                                                                   : 0.0);
        for (std::size_t j = 0; j < s.ps_jobs.size(); ++j) {
            PsJob& pj = s.ps_jobs[j];
            // Retag the STARVED jobs too: the tag is what invalidates a departure already in the list.
            pj.tag = ++ps_tag;
            const double rate = rates[j];
            if (!(rate > 0.0)) continue;
            Event e;
            e.t = now + ((pj.remaining <= 1e-12) ? 1e-12 : pj.remaining / rate);
            e.kind = EV_DEPARTURE;
            e.station = i;
            e.cls = pj.cls;
            e.tag = pj.tag;
            push(e);
        }
    };
 
 
    /**
     * Integrate the work delivered at a state-dependent station since its last
     * change, then re-time every completion it holds.
     *
     * `sd_advance` must run BEFORE the population is written and `sd_reschedule`
     * after, because the rate in force over the elapsed interval is the one the
     * OLD population implied. The stale departures are neutralised by the tag,
     * exactly as a preemption's are.
     */
    auto sd_advance = [&](std::size_t i) {
        StationState& s = S[i];
        if (!s.state_dependent || s.ps) return;
        const double dt = now - s.sd_last_update;
        s.sd_last_update = now;
        if (!(dt > 0.0)) return;
        for (std::size_t sl = 0; sl < s.nservers; ++sl)
            if (s.server_busy[sl]) {
                const double scale = rate_scaling(i, s.server[sl].cls);
                s.server[sl].remaining = std::max(0.0, s.server[sl].remaining - dt * scale);
                s.server[sl].elapsed += dt * scale;
                s.server_start[sl] = now;
            }
    };
    auto sd_reschedule = [&](std::size_t i) {
        StationState& s = S[i];
        if (!s.state_dependent || s.ps) return;
        // The population has just been written, so this is where the new speed
        // takes effect -- for the busy integral exactly as for the departures
        // re-timed below. `set_busy_scale` credits the interval that just ended
        // at the OLD speed before installing the new one.
        acc.set_busy_scale(i, lld_factor(i), now);
        for (std::size_t sl = 0; sl < s.nservers; ++sl)
            if (s.server_busy[sl]) {
                const double scale = rate_scaling(i, s.server[sl].cls);
                // RETAG BEFORE THE ZERO TEST, as ps_reschedule does. The tag is
                // the only thing that invalidates a departure already in the
                // list, so a slot whose rate has fallen to zero must be retagged
                // too -- otherwise the departure scheduled at the OLD rate stays
                // valid and the job leaves in the middle of the outage that was
                // supposed to stop it. Measured on M/M/1/8 with a failing
                // server: mean queue 1.486 against the exact chain's 1.922, with
                // throughput 1.6% HIGH, because the stalled work was delivered
                // anyway.
                s.server_tag[sl] = ++ps_tag;
                if (!(scale > 0.0)) continue;
                Event e;
                e.t = now + s.server[sl].remaining / scale;
                e.kind = EV_DEPARTURE;
                e.station = i;
                e.cls = s.server[sl].cls;
                e.slot = sl;
                e.tag = s.server_tag[sl];
                push(e);
            }
    };
 
 
    /**
     * Whether station `j` can admit one more job of class `r` RIGHT NOW,
     * counting the jobs already blocked in front of it.
     *
     * THE CAPACITY TEST EXCLUDES THE BLOCKED JOBS, and reporting includes
     * them. `acc.qlen` carries both, because the reference's
     * effectiveQueueLength adds basBlockedAtDest/bbsBlockedAtDest for the
     * statistics, while `destinationHasCapacity` reads
     * getTotalCustomersAtStation -- the jobs actually AT the station. Testing
     * against the inclusive count DEADLOCKS a station of capacity one: the job
     * blocked in front of it occupies the very slot it is waiting for, so the
     * station can never make room and the upstream server is never released.
     */
    auto dest_has_room = [&](std::size_t j, std::size_t r) -> bool {
        StationState& d = S[j];
        if (d.role != Role::Queue) return true;
        double total = 0.0;
        for (std::size_t k = 0; k < K; ++k) total += acc.qlen[j][k] - d.blocked_at[k];
        return (total + 1.0 <= d.cap) &&
               (acc.qlen[j][r] - d.blocked_at[r] + 1.0 <= d.classcap[r]);
    };
 
    /** The blocking policy the DESTINATION declares for an arriving class. */
    auto dest_policy = [&](std::size_t j, std::size_t r) -> lang::DropStrategy {
        return (S[j].role == Role::Queue && r < S[j].droprule.size()) ? S[j].droprule[r]
                                                                     : lang::DropStrategy::DROP;
    };
 
    // ---- finite capacity regions --------------------------------------------
    // At most ONE region per station, which is the reference's own restriction:
    // `fcRegionIndices` is a single index per station, and overlapping regions
    // would need a join over their constraints that no codebase implements.
    std::vector<Region> regions;
    std::vector<int> region_of(M, -1);
    for (const auto& rg : sn.regions) {
        Region R;
        R.name = rg.name;
        R.members.assign(M, false);
        R.class_cap.assign(K, -1.0);
        R.class_size.assign(K, 1.0);
        R.class_weight.assign(K, 1.0);
        R.rule.assign(K, lang::DropStrategy::DROP);
        R.jobs.assign(K, 0.0);
        R.blocked.assign(K, 0.0);
        R.dropped.assign(K, 0.0);
        R.tot_jobs.assign(K, 0.0);
        R.tot_weight.assign(K, 0.0);
        R.tot_mem.assign(K, 0.0);
        R.completed.assign(K, 0.0);
        R.resp_sum.assign(K, 0.0);
        R.resp_cnt.assign(K, 0.0);
        for (std::size_t i = 0; i < M && i < rg.members.size(); ++i)
            if (rg.members[i]) {
                R.members[i] = true;
                if (region_of[i] >= 0)
                    throw UnsupportedError(
                        "SolverLDES (native engine): station '" + sn.stations[i].name +
                        "' belongs to more than one finite capacity region");
                region_of[i] = static_cast<int>(regions.size());
            }
        // `cap(i,r)` is the class-r bound at member i and `cap(i,K)` the global
        // one; -1 is the reference's unbounded sentinel and is kept, because it
        // is compared with != -1 to test MEMBERSHIP as well as boundedness.
        for (std::size_t i = 0; i < rg.cap.size(); ++i) {
            if (i >= M || !R.members[i]) continue;
            for (std::size_t r = 0; r < K && r < rg.cap[i].size(); ++r)
                if (rg.cap[i][r] >= 0.0)
                    R.class_cap[r] = (R.class_cap[r] < 0.0) ? rg.cap[i][r]
                                                            : std::min(R.class_cap[r], rg.cap[i][r]);
            if (rg.cap[i].size() > K && rg.cap[i][K] >= 0.0)
                R.global_cap = (R.global_cap < 0.0) ? rg.cap[i][K]
                                                    : std::min(R.global_cap, rg.cap[i][K]);
        }
        for (double mm : rg.maxmem)
            if (mm >= 0.0) R.max_mem = (R.max_mem < 0.0) ? mm : std::min(R.max_mem, mm);
        for (std::size_t r = 0; r < K && r < rg.rule.size(); ++r) R.rule[r] = rg.rule[r];
        for (std::size_t r = 0; r < K && r < rg.size.size(); ++r)
            R.class_size[r] = num_traits<T>::to_double(rg.size[r]);
        for (std::size_t r = 0; r < K && r < rg.weight.size(); ++r)
            R.class_weight[r] = num_traits<T>::to_double(rg.weight[r]);
        for (std::size_t c = 0; c < rg.lincon_A.rows(); ++c) {
            std::vector<double> row;
            for (std::size_t r = 0; r < rg.lincon_A.cols(); ++r)
                row.push_back(num_traits<T>::to_double(rg.lincon_A(c, r)));
            R.lincon_A.push_back(row);
            R.lincon_b.push_back(c < rg.lincon_b.size() ? num_traits<T>::to_double(rg.lincon_b[c])
                                                        : 0.0);
        }
        regions.push_back(R);
    }
 
    /**
     * Jobs parked by a WAITQ region: they have NOT entered the region, so they
     * are outside every station's buffer, but they are still in the system and
     * are charged to the station they were trying to enter.
     */
    struct RegionWaiter {
        std::size_t station = 0;
        Job job;
    };
    std::vector<std::vector<RegionWaiter>> region_wait(regions.size());
 
    auto start_service = [&](std::size_t i, std::size_t slot, const Job& job_in) {
        StationState& s = S[i];
        Job job = job_in;
        // A RESUMED job owes its residual; a RESTARTED one draws afresh and
        // discards what it had done, which is the whole of the PR/PI split.
        if (s.preemptive && !s.resume && job.elapsed > 0.0) {
            job.service = slot_snap(s.svc[job.cls].next_at(g_svc[i][job.cls], now), "service time");
            job.remaining = job.service;
            job.elapsed = 0.0;
        }
        // A POOLED STATION DRAWS THE SERVICE FROM THE POOL, and it can only be
        // drawn here: the law depends on which server takes the job, and that
        // is not known until the slot is assigned. The draw made on arrival is
        // the station's class-level default and stands as the NOMINAL size the
        // size-based disciplines order the waiting room by; the duration
        // actually served is this one. A pool that declares no law of its own
        // for the class keeps that nominal draw, so it costs no extra variate.
        if (s.has_pools && !(s.preemptive && s.resume && job.elapsed > 0.0)) {
            const std::size_t ty = s.server_type[slot];
            if (s.type_has_svc[ty][job.cls]) {
                job.service = slot_snap(
                    s.type_svc[ty][job.cls].next_at(g_hsvc[i][ty][job.cls], now),
                    "service time");
                job.remaining = job.service;
                job.elapsed = 0.0;
            }
        }
        s.server[slot] = job;
        s.server_busy[slot] = true;
        s.server_start[slot] = now;
        s.server_tag[slot] = ++ps_tag;
        acc.update_busy(i, job.cls, now);
        acc.busy[i][job.cls] += 1.0;
        // A ZERO RATE SCHEDULES NOTHING, it does not schedule a departure at
        // t = infinity. A job that enters service while the server is down
        // holds the slot with its residual intact and gets its departure back
        // from the repair's reschedule; pushing `now + remaining/0` instead
        // puts an infinite timestamp in the event list, which the run then
        // pops and processes with now = inf.
        const double scale0 = rate_scaling(i, job.cls);
        if (!(scale0 > 0.0)) return;
        Event e;
        e.t = now + job.remaining / scale0;
        e.kind = EV_DEPARTURE;
        e.station = i;
        e.cls = job.cls;
        e.slot = slot;
        e.tag = s.server_tag[slot];
        push(e);
    };
 
    /**
     * Interrupt the job in `slot` and return it to the waiting room.
     *
     * Its residual work and attained service are carried on the job, so a
     * resume needs no side table; the tag is left stale, which is how the
     * departure event already in the list is neutralised.
     */
    auto preempt = [&](std::size_t i, std::size_t slot) {
        StationState& s = S[i];
        Job job = s.server[slot];
        const double served = now - s.server_start[slot];
        job.remaining = std::max(0.0, job.remaining - served);
        job.elapsed += served;
        s.server_busy[slot] = false;
        s.server_tag[slot] = 0;
        acc.update_busy(i, job.cls, now);
        acc.busy[i][job.cls] -= 1.0;
        // The job stays in the system: its queue length is unchanged, and it
        // re-enters the waiting room under the station's own order.
        buffer_push(i, job);
    };
 
    /**
     * Advance a polling server one leg, or park it.
     *
     * `poll_serve` starts the next job at the buffer the server stands at;
     * `poll_advance` walks it to the next buffer, paying that leg's switchover.
     * The walk stops after a full lap of zero-time legs that found no work, or
     * the server would spin forever at the same simulated instant.
     */
    std::function<void(std::size_t)> poll_serve;
    std::function<void(std::size_t)> poll_advance;
 
    poll_serve = [&](std::size_t i) {
        StationState& s = S[i];
        if (s.server_busy[0] || s.server_blocked[0] || s.server_held[0]) return;
        // Only the class the server currently stands at may be served.
        std::size_t at = s.buffer.size();
        for (std::size_t j = 0; j < s.buffer.size(); ++j)
            if (s.buffer[j].cls == s.poll_at &&
                (at == s.buffer.size() || s.buffer[j].t_arr < s.buffer[at].t_arr))
                at = j;
        if (at == s.buffer.size() || s.poll_budget == 0) {
            poll_advance(i);
            return;
        }
        Job job = s.buffer[at];
        s.buffer.erase(s.buffer.begin() + static_cast<std::ptrdiff_t>(at));
        std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
        if (s.poll_type != lang::PollingType::EXHAUSTIVE) --s.poll_budget;
        start_service(i, 0, job);
    };
 
    poll_advance = [&](std::size_t i) {
        StationState& s = S[i];
        for (std::size_t step = 0; step < K; ++step) {
            const std::size_t from = s.poll_at;
            const std::size_t next = (from + 1) % K;
            // Takagi's r_i: the leg OUT of `from` is charged the switchover
            // stored against `from`, not against the buffer being entered.
            if (s.has_switchover[from]) {
                s.poll_switching = true;
                Event e;
                e.t = now + s.switchover[from].next(g_aux[i]);
                e.kind = EV_SWITCHOVER;
                e.station = i;
                e.cls = next;
                push(e);
                return;
            }
            s.poll_at = next;
            s.poll_budget = (s.poll_type == lang::PollingType::EXHAUSTIVE)
                                ? std::numeric_limits<std::size_t>::max()
                                : ((s.poll_type == lang::PollingType::KLIMITED) ? s.poll_k : 1);
            bool has_work = false;
            for (const Job& j : s.buffer)
                if (j.cls == next) {
                    has_work = true;
                    break;
                }
            if (has_work) {
                poll_serve(i);
                return;
            }
        }
        // A whole lap of zero-time legs found nothing. Parking is the only
        // alternative to spinning at this instant, and it is unobservable: the
        // next arrival restarts the server.
        s.poll_parked = true;
    };
 
    // ---- admission ---------------------------------------------------------
 
    /**
     * Reschedule the aggregate completion of a pass-and-swap station.
     *
     * ONE CLOCK FOR THE WHOLE STATION, at the total rate mu(c) of its ordered
     * list, and the position that completes is drawn from the per-position
     * INCREMENTS. Scheduling one clock per job instead would need each job's
     * marginal rate to be constant, which is exactly what an order-independent
     * station does not have: every arrival and departure changes every other
     * job's rate.
     */
    std::function<void(std::size_t)> pas_reschedule = [&](std::size_t i) {
        StationState& s = S[i];
        if (s.pas_list.empty()) return;
        std::vector<std::size_t> seq;
        for (const Job& j : s.pas_list) seq.push_back(j.cls);
        const double total = s.pas_rate(seq);
        if (!(total > 0.0)) return;
        s.pas_tag = ++ps_tag;
        Event e;
        e.t = now + (-std::log(uniform01(g_aux[i])) / total);
        e.kind = EV_DEPARTURE;
        e.station = i;
        e.cls = 0;
        e.tag = s.pas_tag;
        e.slot = std::numeric_limits<std::size_t>::max();
        push(e);
    };
 
    /**
     * The position that completes, drawn from the per-position rate increments,
     * and the job that actually departs after the swap chain.
     */
    auto pas_complete = [&](std::size_t i) -> Job {
        StationState& s = S[i];
        std::vector<std::size_t> seq;
        for (const Job& j : s.pas_list) seq.push_back(j.cls);
        std::vector<double> inc(seq.size(), 0.0);
        double prev = 0.0, total = 0.0;
        for (std::size_t p = 0; p < seq.size(); ++p) {
            std::vector<std::size_t> pre(seq.begin(),
                                         seq.begin() + static_cast<std::ptrdiff_t>(p + 1));
            const double cur = s.pas_rate(pre);
            inc[p] = std::max(0.0, cur - prev);
            prev = cur;
            total += inc[p];
        }
        std::size_t pos = 0;
        if (total > 0.0) {
            const double u = uniform01(g_routing) * total;
            double acc2 = 0.0;
            for (std::size_t p = 0; p < inc.size(); ++p) {
                acc2 += inc[p];
                if (u < acc2) {
                    pos = p;
                    break;
                }
            }
        }
        // The swap chain: from `pos`, walk forward to the last position the swap
        // graph can reach. That job departs and the chain shifts back into the
        // hole at `pos`. An empty graph makes the chain trivial, which is the
        // plain order-independent station.
        std::vector<std::size_t> chain(1, pos);
        std::size_t moving = seq[pos], cur = pos;
        while (true) {
            std::size_t nxt = seq.size();
            for (std::size_t j = cur + 1; j < seq.size(); ++j) {
                const bool swappable =
                    s.pas_swap.empty() ||
                    (moving < s.pas_swap.size() && seq[j] < s.pas_swap[moving].size() &&
                     s.pas_swap[moving][seq[j]]);
                if (swappable) {
                    nxt = j;
                    break;
                }
            }
            if (nxt >= seq.size()) break;
            chain.push_back(nxt);
            moving = seq[nxt];
            cur = nxt;
        }
        Job departing = s.pas_list[chain.back()];
        std::vector<Job> arr = s.pas_list;
        for (std::size_t k = 0; k + 1 < chain.size(); ++k)
            arr[chain[k + 1]] = s.pas_list[chain[k]];
        std::vector<Job> survivors;
        for (std::size_t k = 0; k < arr.size(); ++k)
            if (k != chain.front()) survivors.push_back(arr[k]);
        s.pas_list = survivors;
        return departing;
    };
    /**
     * Admit a job at station `i`, or drop it when the station is full.
     *
     * `from` is the station the job comes from, or M for an arrival that enters
     * from outside the network. IT IS NOT BOOKKEEPING: a hop BETWEEN TWO
     * MEMBERS of the same finite capacity region does not cross the region
     * boundary, so the job must neither be re-tested against the region's caps
     * nor counted into it again. Without that test an intra-region hop adds one
     * to the region on arrival and removes nothing on departure, so occupancy
     * grows by one per hop until every arrival is refused and the run never
     * completes another job.
     */
    std::function<bool(std::size_t, Job, std::size_t)> admit =
        [&](std::size_t i, Job job, std::size_t from) -> bool {
        StationState& s = S[i];
        if (s.off[job.cls])
            throw InputError("SolverLDES (native engine): a job of class '" +
                             sn.classes[job.cls].name + "' reached station '" +
                             sn.stations[i].name + "', which does not serve it");
        // THE REGION IS CONSULTED FIRST. Its cap spans stations, so a job may
        // be refused entry to a station that has room of its own; testing the
        // station first and the region second would admit into the buffer and
        // then have to unwind it. The test is on MEMBERSHIP ALONE, exactly as
        // on the departure side: a Delay is capped by the region like any other
        // member, and gating entry on Role::Queue while the exit is ungated
        // drains the occupancy below zero, after which nothing is ever refused.
        int rgi = region_of[i];
        if (rgi >= 0 && from < M && region_of[from] == rgi) rgi = -1;  // stays inside
        if (rgi >= 0) {
            Region& R = regions[static_cast<std::size_t>(rgi)];
            if (R.would_exceed(job.cls)) {
                R.update(now);
                if (R.drops(job.cls)) {
                    R.dropped[job.cls] += 1.0;
                    dropped[i][job.cls] += 1.0;
                    return false;
                }
                // WAITQ: parked outside the region, still in the system.
                RegionWaiter w;
                w.station = i;
                w.job = job;
                w.job.t_arr = now;
                w.job.id = ++job_id;
                region_wait[static_cast<std::size_t>(rgi)].push_back(w);
                R.blocked[job.cls] += 1.0;
                return true;
            }
        }
 
        // Kendall's K counts the jobs HELD by the station, waiting and in
        // service alike; a per-class buffer caps its own class on top of that.
        if (s.role == Role::Queue) {
            double total = 0.0;
            for (std::size_t r = 0; r < K; ++r) total += acc.qlen[i][r] - s.blocked_at[r];
            if (total + 1.0 > s.cap ||
                acc.qlen[i][job.cls] - s.blocked_at[job.cls] + 1.0 > s.classcap[job.cls]) {
                // A RETRIAL class does not drop: it joins the orbit and
                // re-offers itself later. The orbit is a separate population --
                // not at the station, not in its capacity, not in its queue
                // length -- so nothing here is charged to the station.
                if (s.has_retrial[job.cls] &&
                    (s.max_attempts[job.cls] <= 0 ||
                     job.attempts < s.max_attempts[job.cls])) {
                    const double dt = now - s.orbit_last;
                    if (dt > 0.0)
                        for (std::size_t r = 0; r < K; ++r) s.tot_orbit[r] += s.orbit_size[r] * dt;
                    s.orbit_last = now;
                    s.orbit_size[job.cls] += 1.0;
                    Event e;
                    e.t = now + s.retrial[job.cls].next(g_svc[i][job.cls]);
                    e.kind = EV_RETRIAL;
                    e.station = i;
                    e.cls = job.cls;
                    e.job = job;
                    e.job.attempts = job.attempts + 1;
                    push(e);
                    return true;
                }
                if (s.has_retrial[job.cls]) s.retrial_lost[job.cls] += 1.0;
                dropped[i][job.cls] += 1.0;
                return false;
            }
            // BALKING is a refusal to JOIN, decided on what the arrival sees.
            // It is not a drop for want of room: the station may have plenty,
            // and the job declines anyway. Counted separately for that reason.
            if (!s.balk[job.cls].empty()) {
                double p = 0.0;
                for (const StationState::BalkRule& br : s.balk[job.cls])
                    if (total >= br.min_jobs && (br.max_jobs < 0.0 || total <= br.max_jobs))
                        p = std::max(p, br.probability);
                if (p > 0.0 && uniform01(g_routing) < p) {
                    balked[i][job.cls] += 1.0;
                    return false;
                }
            }
        }
        job.t_arr = now;
        job.priority = classprio[job.cls];
        // A TIME-INHOMOGENEOUS service is drawn FROM THE CURRENT INSTANT: its
        // rate depends on absolute time, so a duration alone does not
        // determine it. Every other family ignores the argument.
        job.service = slot_snap(S[i].svc[job.cls].next_at(g_svc[i][job.cls], now), "service time");
        job.remaining = job.service;
        job.elapsed = 0.0;
        job.rank = uniform01(g_routing);
        job.deadline = now + classdeadline[job.cls];
        job.id = ++job_id;
        if (rgi >= 0) {
            Region& R = regions[static_cast<std::size_t>(rgi)];
            R.update(now);
            R.enter(job.cls);
        }
        sd_advance(i);
        acc.update_qlen(i, job.cls, now);
        acc.qlen[i][job.cls] += 1.0;
        // Counted where the job JOINS the station, so a drop, a balk or a
        // reneging leaves AN above TN, which is the difference AN exists to show.
        acc.arrived[i][job.cls] += 1.0;
        bp.track(i, job.cls, +1, now);
 
        if (s.role == Role::Delay) {
            if (has_removal_signal) delay_live[i].push_back(job);
            Event e;
            e.t = now + job.service;
            e.kind = EV_DEPARTURE;
            e.station = i;
            e.cls = job.cls;
            e.slot = std::numeric_limits<std::size_t>::max();
            e.job = job;
            push(e);
            return true;
        }
        if (s.ps) {
            ps_advance(i);
            // LPS admits at most `limit` jobs into the shared server; the rest
            // wait in the buffer under the same order a PS station would use.
            if (s.lps_limit > 0 && s.ps_jobs.size() >= s.lps_limit) {
                buffer_push(i, job);
                return true;
            }
            PsJob pj;
            pj.cls = job.cls;
            pj.priority = job.priority;
            pj.t_arr = job.t_arr;
            pj.t_sys = job.t_sys;
            pj.total = job.service;
            pj.remaining = job.service;
            pj.parent = job.parent;
            s.ps_jobs.push_back(pj);
            ps_reschedule(i);
            return true;
        }
        // A station that powers down: the first arrival to a cold server pays
        // the setup before any service starts, and one arriving inside the
        // delay-off window pays nothing because the server never shut down.
        if (s.has_setup && !s.setup_on) {
            buffer_push(i, job);
            if (!s.setup_running) {
                s.setup_running = true;
                Event e;
                e.t = now + s.setup_time.next(g_aux[i]);
                e.kind = EV_SETUP;
                e.station = i;
                e.cls = job.cls;
                push(e);
            }
            return true;
        }
        if (s.pas) {
            s.pas_list.push_back(job);
            acc.update_busy(i, job.cls, now);
            acc.busy[i][job.cls] += 1.0;
            pas_reschedule(i);
            return true;
        }
        if (s.polling) {
            buffer_push(i, job);
            if (s.poll_parked && !s.poll_switching) {
                s.poll_parked = false;
                if (job.cls == s.poll_at)
                    poll_serve(i);
                else
                    poll_advance(i);
            } else if (!s.server_busy[0] && !s.poll_switching) {
                poll_serve(i);
            }
            return true;
        }
        {
            const std::size_t sl = free_slot_for(i, job.cls);
            if (sl < s.nservers) {
                start_service(i, sl, job);
                sd_reschedule(i);
                return true;
            }
        }
        // Every server is busy, or on a pooled station every server THIS CLASS
        // MAY USE is. A preemptive station now asks whether this arrival
        // displaces one of them; the freed slot is taken by the arrival itself,
        // not by the head of the queue, because the victim has just been pushed
        // there and could otherwise be picked straight back.
        if (s.preemptive) {
            std::vector<Job> held;
            std::vector<std::size_t> slot_of;
            for (std::size_t sl = 0; sl < s.nservers; ++sl)
                if (s.server_busy[sl] && !s.server_blocked[sl]) {
                    // ONLY A SLOT THE ARRIVAL COULD THEN OCCUPY is a candidate
                    // victim: preempting a job on a pool that does not accept
                    // this class frees a server the arrival cannot use, so the
                    // victim would be displaced for nothing.
                    if (s.has_pools && !s.type_compat[s.server_type[sl]][job.cls]) continue;
                    Job h = s.server[sl];
                    // The residual as of NOW, which is what every size-based
                    // rule compares against; the stored value is as of the
                    // instant service started.
                    h.remaining = std::max(0.0, h.remaining - (now - s.server_start[sl]));
                    h.elapsed += now - s.server_start[sl];
                    held.push_back(h);
                    slot_of.push_back(sl);
                }
            const std::size_t v = preemption_victim(s.sched, held, job,
                                                    static_cast<double>(s.nservers), now);
            if (v < held.size()) {
                const std::size_t sl = slot_of[v];
                preempt(i, sl);
                start_service(i, sl, job);
                sd_reschedule(i);
                return true;
            }
        }
        buffer_push(i, job);
        // RENEGING times out a job that is WAITING, not one in service: a job
        // that reached a server has been served and no longer abandons. The
        // timer is therefore armed here and nowhere else, and it is re-armed
        // if a preemption sends the job back to the waiting room.
        if (s.has_patience[job.cls]) {
            Event e;
            e.t = now + s.patience[job.cls].next(g_svc[i][job.cls]);
            e.kind = EV_RENEGE;
            e.station = i;
            e.cls = job.cls;
            e.tag = job.id;
            push(e);
        }
        sd_reschedule(i);
        return true;
    };
 
    const std::map<std::size_t, double> kNoWeights;
    /** Every index of `tab`, so JSQ can be told to consider all of them. */
    std::vector<std::size_t> tab_all;
 
    /**
     * PROB: one draw against the static row.
     *
     * The mass is NOT renormalized before the draw. The reference scales the
     * uniform by the row total instead, so a row summing to one only up to
     * rounding still lands, and the two engines consume the routing stream
     * identically.
     */
    auto draw_prob = [&](const std::vector<RouteDest>& tab) -> std::size_t {
        double total = 0.0;
        for (std::size_t k = 0; k < tab.size(); ++k) total += tab[k].mass;
        const double u = uniform01(g_routing) * total;
        double cum = 0.0;
        std::size_t pick = 0;
        for (pick = 0; pick + 1 < tab.size(); ++pick) {
            cum += tab[pick].mass;
            if (u <= cum) break;
        }
        return pick;
    };
 
    // ---- state-dependent routing (Krzesinski 1987) ----------------------------
    // The coefficients are derived ONCE: `pfqn_sdrcoeff` validates the
    // declaration and its population bounds, and repeating that at every
    // departure would make the validation the cost of the run.
    //
    // INDEX SPACE. The probabilities are read off the STATION-indexed structure
    // against station populations, and the destination split is mapped back
    // through the NODE-indexed twin, exactly as `state.h`'s `rt_state` does. The
    // two spaces must not be mixed: `pfqn_sdrprob` uses `sdr.branch[b][k]` as
    // indices into the population vector, so the structure and the vector have
    // to agree, while `entryOf`/`departure` are compared against node indices.
    const bool has_sdr = !sn.sdr.empty() && !sn.sdr_nodes.empty();
    pfqn::SdrCoeff sdr_coeff;
    if (has_sdr) sdr_coeff = pfqn::pfqn_sdrcoeff(sn.sdr);
    std::vector<double> sdr_pop, sdr_mass;
 
    /**
     * Eq. (10) of Krzesinski (1987) at the CURRENT marking: the split out of the
     * entry centre into the branch entries, plus the residual that returns the
     * customer to the departure centre.
     *
     * A Source is reported as zero rather than as the infinite reservoir its
     * encoding carries, which is the `station_populations` rule and the one the
     * enumerating solvers use, so the sample path here has an exact counterpart
     * in SolverNC rather than merely a similar one.
     */
    auto draw_sdr = [&](const std::vector<RouteDest>& tab) -> std::size_t {
        sdr_pop.assign(M, 0.0);
        for (std::size_t i = 0; i < M; ++i) {
            if (S[i].role == Role::Source) continue;
            double held = 0.0;
            for (std::size_t q = 0; q < K; ++q) held += acc.qlen[i][q];
            sdr_pop[i] = held;
        }
        const std::vector<double> Pb = pfqn::pfqn_sdrprob(sdr_coeff, sdr_pop);
        const double Ped = pfqn::pfqn_sdrped(Pb);
 
        sdr_mass.assign(tab.size(), 0.0);
        double total = 0.0;
        for (std::size_t i = 0; i < tab.size(); ++i) {
            const std::size_t dnode = tab[i].node - 1;  // 0-based, as sdr_nodes is
            double p = 0.0;
            for (std::size_t b = 1; b < sn.sdr_nodes.branch.size(); ++b)
                if (sn.sdr_nodes.entryOf[b] == dnode) p += Pb[b];
            if (sn.sdr_nodes.departure == dnode) p += Ped;
            // Eq. (10) returns zero at a branch closed to new arrivals; the
            // clamp only guards the residual against round-off.
            if (p < 0.0) p = 0.0;
            sdr_mass[i] = p;
            total += p;
        }
        // The branch probabilities and the residual sum to one over the branch
        // entries and the departure centre, so this normalization is an identity
        // whenever the declaration and the links agree. It is kept so the draw
        // stays a distribution if they do not, and the fallback below keeps a
        // model whose links reach neither from dividing by zero.
        if (!(total > 0.0)) return draw_prob(tab);
        const double u = uniform01(g_routing) * total;
        double cum = 0.0;
        std::size_t pick = 0;
        for (pick = 0; pick + 1 < tab.size(); ++pick) {
            cum += sdr_mass[pick];
            if (u <= cum) break;
        }
        return pick;
    };
 
    /**
     * The candidate holding the FEWEST jobs, ties broken uniformly.
     *
     * A destination that is not a service station holds nothing to compare, so
     * it is skipped rather than counted as empty -- a Sink would otherwise win
     * every comparison and absorb the whole stream. When no candidate is a
     * station the first one is taken, which is the reference's fallback.
     */
    auto shortest_queue = [&](const std::vector<RouteDest>& tab,
                              const std::vector<std::size_t>& cand) -> std::size_t {
        double best = std::numeric_limits<double>::infinity();
        std::vector<std::size_t> tied;
        for (std::size_t c = 0; c < cand.size(); ++c) {
            const std::size_t k = cand[c];
            if (tab[k].station >= M || S[tab[k].station].role == Role::Source) continue;
            double held = 0.0;
            for (std::size_t q = 0; q < K; ++q) held += acc.qlen[tab[k].station][q];
            if (held < best) {
                best = held;
                tied.clear();
                tied.push_back(k);
            } else if (held == best) {
                tied.push_back(k);
            }
        }
        if (tied.empty()) return cand.empty() ? 0 : cand[0];
        if (tied.size() == 1) return tied[0];
        std::size_t at = static_cast<std::size_t>(uniform01(g_routing) *
                                                  static_cast<double>(tied.size()));
        if (at >= tied.size()) at = tied.size() - 1;
        return tied[at];
    };
 
    /**
     * One hop out of node `inode` in class `r`: the DISPATCHER picks the node,
     * the class-switch row picks the class.
     *
     * NEITHER STAGE DRAWS WHEN ITS OUTCOME IS DETERMINED. A single destination
     * and a single arriving class each cost nothing, which is the reference's
     * rule and not an optimization: a routing stream that advanced on a forced
     * hop would put every subsequent draw out of step with the Java engine and
     * make a seeded run of a tandem irreproducible across the two.
     */
    auto draw_node_route = [&](std::size_t inode, std::size_t r) -> RouteEntry {
        const std::vector<RouteDest>& tab = nroute[inode][r];
        if (tab.empty())
            throw InputError("SolverLDES (native engine): node '" + sn.nodes[inode - 1].name +
                             "' has no routing for class '" + sn.classes[r].name + "'");
        std::size_t pick = 0;
        if (tab.size() > 1) {
            if (tab_all.size() != tab.size()) {
                tab_all.resize(tab.size());
                for (std::size_t k = 0; k < tab.size(); ++k) tab_all[k] = k;
            }
            switch (node_routing[inode][r]) {
                case lang::RoutingStrategy::RAND: {
                    pick = static_cast<std::size_t>(uniform01(g_routing) *
                                                    static_cast<double>(tab.size()));
                    if (pick >= tab.size()) pick = tab.size() - 1;
                    break;
                }
                case lang::RoutingStrategy::RROBIN: {
                    pick = rr_counter[inode][r] % tab.size();
                    ++rr_counter[inode][r];
                    break;
                }
                case lang::RoutingStrategy::WRROBIN: {
                    // The weights become an INTEGER schedule and the counter
                    // walks it, so a 3:1 split really is three then one and not
                    // a 0.75 coin. A fractional weight is scaled by 1000 and a
                    // positive weight never rounds to zero, or a link declared
                    // open would never be offered a job.
                    const std::map<std::size_t, double>& w =
                        (sn.nodes[inode - 1].routing_weights.size() > r)
                            ? sn.nodes[inode - 1].routing_weights[r]
                            : kNoWeights;
                    std::vector<double> raw(tab.size(), 0.0);
                    bool integral = true;
                    for (std::size_t k = 0; k < tab.size(); ++k) {
                        const std::map<std::size_t, double>::const_iterator it = w.find(tab[k].node);
                        raw[k] = (it == w.end()) ? 0.0 : it->second;
                        if (raw[k] < 0.0 || raw[k] != std::floor(raw[k])) integral = false;
                    }
                    std::vector<std::size_t> quota(tab.size(), 0);
                    std::size_t total = 0;
                    for (std::size_t k = 0; k < tab.size(); ++k) {
                        double v = integral ? raw[k] : raw[k] * 1000.0;
                        std::size_t q = (v > 0.0) ? static_cast<std::size_t>(v) : 0;
                        if (raw[k] > 0.0 && q < 1) q = 1;
                        quota[k] = q;
                        total += q;
                    }
                    if (total == 0) {
                        // No usable weights: the dispatcher degenerates to plain
                        // round robin, which is the reference's fallback.
                        pick = rr_counter[inode][r] % tab.size();
                        ++rr_counter[inode][r];
                        break;
                    }
                    std::size_t counter = rr_counter[inode][r] % total;
                    ++rr_counter[inode][r];
                    std::size_t cum = 0;
                    for (pick = 0; pick + 1 < tab.size(); ++pick) {
                        cum += quota[pick];
                        if (counter < cum) break;
                    }
                    break;
                }
                case lang::RoutingStrategy::SDR: {
                    // Unlike JSQ and SQ(d) this routing keeps a product form of
                    // its own, so this sample path has an exact counterpart in
                    // solver_nc_sdr. The two agree only if both read the same
                    // populations; see draw_sdr. A routing named SDR with no
                    // declaration behind it is a malformed model every reader
                    // refuses, and draw_prob is what it falls back to rather
                    // than dividing by an absent structure.
                    pick = has_sdr ? draw_sdr(tab) : draw_prob(tab);
                    break;
                }
                case lang::RoutingStrategy::JSQ: {
                    pick = shortest_queue(tab, tab_all);
                    break;
                }
                case lang::RoutingStrategy::SQ: {
                    // POWER OF d: sample d candidates WITHOUT replacement, then
                    // take the shortest among them. With d at least the number
                    // of destinations the sampling is vacuous and it IS JSQ.
                    std::size_t d = (sn.nodes[inode - 1].routing_param.size() > r &&
                                     sn.nodes[inode - 1].routing_param[r] > 0)
                                        ? static_cast<std::size_t>(
                                              sn.nodes[inode - 1].routing_param[r])
                                        : 2;
                    if (tab.size() <= d) {
                        pick = shortest_queue(tab, tab_all);
                        break;
                    }
                    std::vector<std::size_t> avail(tab.size());
                    for (std::size_t k = 0; k < tab.size(); ++k) avail[k] = k;
                    std::vector<std::size_t> cand;
                    for (std::size_t k = 0; k < d && !avail.empty(); ++k) {
                        std::size_t at = static_cast<std::size_t>(
                            uniform01(g_routing) * static_cast<double>(avail.size()));
                        if (at >= avail.size()) at = avail.size() - 1;
                        cand.push_back(avail[at]);
                        avail.erase(avail.begin() + static_cast<std::ptrdiff_t>(at));
                    }
                    pick = shortest_queue(tab, cand);
                    break;
                }
                default: {
                    pick = draw_prob(tab);
                    break;
                }
            }
        }
        const RouteDest& d = tab[pick];
        RouteEntry e;
        e.node = d.node;
        e.station = d.station;
        e.sink = d.sink;
        e.cls = d.cls[0].first;
        if (d.cls.size() > 1) {
            double total = 0.0;
            for (std::size_t k = 0; k < d.cls.size(); ++k) total += d.cls[k].second;
            const double u = uniform01(g_routing) * total;
            double cum = 0.0;
            for (std::size_t k = 0; k < d.cls.size(); ++k) {
                cum += d.cls[k].second;
                if (u <= cum || k + 1 == d.cls.size()) {
                    e.cls = d.cls[k].first;
                    break;
                }
            }
        }
        return e;
    };
 
    /**
     * The class a job routed out of `node` under `cls` will ARRIVE the true
     * destination station under, resolving through any ClassSwitch/Router/
     * Logger chain `draw_node_route` itself stops one hop short of.
     *
     * `draw_node_route(i, r).cls` names the class on the FIRST edge out of a
     * station, and when `link()` synthesized a ClassSwitch node for that edge
     * (any route whose class differs from station to station), that first
     * edge's class is the PRE-switch one -- the switch happens INSIDE the
     * synthesized node, one hop later. The synchronous-call park check below
     * needs the class the job actually lands under, or a callee's own reply
     * (its edge into the synthesized node still reads as the CALL class) is
     * misread as a brand-new call: the reply's `job.call` is then overwritten
     * with a fresh id, the caller's real call is never answered, and its
     * server stays held while the buffer behind it grows without bound.
     *
     * Consumes no extra draw when every hop is forced, which a link()-
     * synthesized pass-through node always is (`draw_node_route` never spends
     * a draw on a single-destination table), so the routing stream stays in
     * step with a caller that already resolved through this same node.
     */
    auto resolve_final_cls = [&](std::size_t node, std::size_t cls) -> std::size_t {
        std::size_t at_node = node, at_cls = cls;
        for (std::size_t hop = 0; hop <= nnodes && at_node != 0; ++hop) {
            const NodeType dt = sn.nodes[at_node - 1].nodetype;
            if (dt != NodeType::ClassSwitch && dt != NodeType::Router && dt != NodeType::Logger)
                break;
            const RouteEntry e = draw_node_route(at_node, at_cls);
            at_cls = e.cls;
            at_node = e.node;
        }
        return at_cls;
    };
 
    // ---- busy periods -------------------------------------------------------
    // Every station is a target set of its own, every declared subnetwork is a
    // set, and each is measured aggregated and then per class, which is the
    // order `initializeBusyPeriods` builds them in and the order the result
    // rows carry.
    if (o.busy_period_orders > 0) {
        std::vector<std::vector<std::size_t>> sets;
        std::vector<std::string> set_names;
        for (std::size_t i = 0; i < M; ++i) {
            if (S[i].role == Role::Source || S[i].role == Role::Synchronization) continue;
            sets.push_back(std::vector<std::size_t>(1, i));
            set_names.push_back(sn.stations[i].name);
        }
        for (const std::vector<std::size_t>& sub : o.busy_period_subnets) {
            std::string nm;
            for (std::size_t s2 : sub) {
                if (s2 >= M || S[s2].role == Role::Source)
                    throw InputError("SolverLDES (native engine): station " +
                                     std::to_string(s2) +
                                     " of a busy period subnetwork is not a service station");
                nm += (nm.empty() ? "" : "+") + sn.stations[s2].name;
            }
            sets.push_back(sub);
            set_names.push_back(nm);
        }
        std::vector<BpTarget> targets;
        for (std::size_t si = 0; si < sets.size(); ++si) {
            BpTarget t;
            t.stations = sets[si];
            t.job_class = -1;
            t.name = set_names[si];
            targets.push_back(t);
            for (std::size_t r = 0; r < K; ++r) {
                BpTarget tc;
                tc.stations = sets[si];
                tc.job_class = static_cast<int>(r);
                tc.name = set_names[si] + ":" + sn.classes[r].name;
                targets.push_back(tc);
            }
        }
        bp.init(targets, o.busy_period_orders, M);
    }
 
 
    /**
     * Hand on the jobs blocked in front of `j`, oldest first, now that it has
     * room.
     *
     * OLDEST FIRST is the reference's rule and it matters: releasing the newest
     * would starve a server that has been blocked longest, which no blocking
     * discipline does and which shows up as an unbounded blocked time at one
     * station while another cycles freely.
     */
    std::function<void(std::size_t)> release_blocked = [&](std::size_t j) {
        bool progress = true;
        while (progress) {
            progress = false;
            std::size_t best_i = M, best_sl = 0;
            double oldest = std::numeric_limits<double>::infinity();
            for (std::size_t a = 0; a < M; ++a) {
                // A SHARING station has no per-server slots at all, so it has
                // no slot to block and none of these vectors is allocated for
                // it. Indexing them by `nservers` reads off the end.
                if (S[a].role != Role::Queue || S[a].ps || S[a].server_blocked.empty()) continue;
                for (std::size_t sl = 0; sl < S[a].nservers; ++sl)
                    if (S[a].server_blocked[sl] && S[a].blocked_dest[sl] == j &&
                        S[a].blocked_job[sl].t_arr < oldest) {
                        oldest = S[a].blocked_job[sl].t_arr;
                        best_i = a;
                        best_sl = sl;
                    }
            }
            if (best_i >= M) break;
            const std::size_t dcls = S[best_i].blocked_dest_cls[best_sl];
            if (!dest_has_room(j, dcls)) break;
 
            StationState& src = S[best_i];
            Job moved = src.blocked_job[best_sl];
            src.server_blocked[best_sl] = false;
            src.server_busy[best_sl] = false;
            src.server_tag[best_sl] = 0;
            acc.update_qlen(j, dcls, now);
            acc.qlen[j][dcls] -= 1.0;
            S[j].blocked_at[dcls] -= 1.0;
 
            // The job now leaves the upstream station for good.
            acc.update_busy(best_i, moved.cls, now);
            if (acc.busy[best_i][moved.cls] > 0.0) acc.busy[best_i][moved.cls] -= 1.0;
 
            Job fresh;
            fresh.cls = dcls;
            fresh.t_sys = moved.t_sys;
            fresh.parent = moved.parent;
            admit(j, fresh, best_i);
 
            // The freed server takes its own next waiting job -- on a pooled
            // station, the first one in service order that its pool accepts.
            if (buffer_has_for_slot(best_i, best_sl)) {
                Job nextjob = buffer_pop_for_slot(best_i, best_sl);
                start_service(best_i, best_sl, nextjob);
            }
            progress = true;
        }
    };
 
 
    /**
     * Let the jobs parked by region `rg` in, oldest first, now that it has
     * room.
     *
     * A waiter re-enters through `admit`, which re-tests every constraint: the
     * region may have freed one slot while the waiter's own class cap or the
     * linear constraints still bind, and admitting on the global count alone
     * would breach them.
     */
    std::function<void(std::size_t)> release_region = [&](std::size_t rg) {
        std::vector<RegionWaiter>& q = region_wait[rg];
        bool progress = true;
        while (progress && !q.empty()) {
            progress = false;
            for (std::size_t w = 0; w < q.size(); ++w) {
                Region& R = regions[rg];
                if (R.would_exceed(q[w].job.cls)) continue;
                RegionWaiter taken = q[w];
                R.update(now);
                R.blocked[taken.job.cls] -= 1.0;
                q.erase(q.begin() + static_cast<std::ptrdiff_t>(w));
                Job fresh;
                fresh.cls = taken.job.cls;
                fresh.t_sys = taken.job.t_sys;
                admit(taken.station, fresh, M);
                progress = true;
                break;
            }
        }
    };
 
    // ---- G-network signals --------------------------------------------------
 
    /** The jobs station `i` actually HOLDS, waiting and in service alike. */
    auto signal_held = [&](std::size_t i) -> std::size_t {
        const StationState& s = S[i];
        if (s.role == Role::Delay) return delay_live[i].size();
        if (s.pas) return s.pas_list.size();
        std::size_t n = s.buffer.size();
        if (s.ps) return n + s.ps_jobs.size();
        for (std::size_t sl = 0; sl < s.nservers; ++sl)
            if (s.server_busy[sl]) ++n;
        return n;
    };
 
    /**
     * Remove ONE job from station `i` under `policy`, returning its class or
     * `K` when the station holds nothing.
     *
     * THE POLICY IS TWO-TIER, as in the reference: FCFS and LCFS rank the
     * WAITING line by age and only reach into the servers once nobody waits,
     * while RANDOM draws uniformly over waiting and in-service jobs together.
     * A station whose occupants are all in service -- Delay, PS, PAS -- has no
     * waiting line for the age tiers to rank, so there the three policies
     * differ only in which of the equally-served jobs is taken.
     *
     * The victim is drawn WITHOUT REGARD TO CLASS. That is Gelenbe's negative
     * customer and what the Java engine and SolverMAM both do; `sn.signaltarget`
     * narrows the eligible class only in the state-space solvers, whose state
     * carries the per-class composition this one draws over directly.
     */
    auto signal_remove_one = [&](std::size_t i, lang::RemovalPolicy policy) -> std::size_t {
        StationState& s = S[i];
 
        // A DELAY holds every occupant in service, in the registry.
        if (s.role == Role::Delay) {
            std::vector<Job>& live = delay_live[i];
            if (live.empty()) return K;
            std::size_t at = 0;
            if (policy == lang::RemovalPolicy::FCFS) {
                for (std::size_t j = 1; j < live.size(); ++j)
                    if (live[j].t_arr < live[at].t_arr) at = j;
            } else if (policy == lang::RemovalPolicy::LCFS) {
                for (std::size_t j = 1; j < live.size(); ++j)
                    if (live[j].t_arr > live[at].t_arr) at = j;
            } else {
                at = static_cast<std::size_t>(uniform01(g_routing) *
                                              static_cast<double>(live.size()));
                if (at >= live.size()) at = live.size() - 1;
            }
            const std::size_t rc = live[at].cls;
            live.erase(live.begin() + static_cast<std::ptrdiff_t>(at));
            return rc;
        }
 
        // A PASS-AND-SWAP list: one aggregate clock covers the whole station,
        // so the survivors must be re-timed rather than left on a rate the
        // list no longer has.
        if (s.pas) {
            if (s.pas_list.empty()) return K;
            std::size_t at = 0;
            if (policy == lang::RemovalPolicy::FCFS) {
                for (std::size_t j = 1; j < s.pas_list.size(); ++j)
                    if (s.pas_list[j].t_arr < s.pas_list[at].t_arr) at = j;
            } else if (policy == lang::RemovalPolicy::LCFS) {
                for (std::size_t j = 1; j < s.pas_list.size(); ++j)
                    if (s.pas_list[j].t_arr > s.pas_list[at].t_arr) at = j;
            } else {
                at = static_cast<std::size_t>(uniform01(g_routing) *
                                              static_cast<double>(s.pas_list.size()));
                if (at >= s.pas_list.size()) at = s.pas_list.size() - 1;
            }
            const std::size_t rc = s.pas_list[at].cls;
            s.pas_list.erase(s.pas_list.begin() + static_cast<std::ptrdiff_t>(at));
            acc.update_busy(i, rc, now);
            if (acc.busy[i][rc] > 0.0) acc.busy[i][rc] -= 1.0;
            pas_reschedule(i);
            return rc;
        }
 
        const std::size_t waiting = s.buffer.size();
 
        // A SHARING station: `ps_jobs` are in service, `buffer` holds only the
        // jobs an LPS admission cap has left outside.
        if (s.ps) {
            const std::size_t inservice = s.ps_jobs.size();
            if (waiting + inservice == 0) return K;
            bool from_wait;
            if (policy == lang::RemovalPolicy::RANDOM) {
                const std::size_t pick = static_cast<std::size_t>(
                    uniform01(g_routing) * static_cast<double>(waiting + inservice));
                from_wait = pick < waiting;
            } else {
                from_wait = waiting > 0;
            }
            if (from_wait) {
                std::size_t at = 0;
                if (policy == lang::RemovalPolicy::FCFS) {
                    for (std::size_t j = 1; j < s.buffer.size(); ++j)
                        if (s.buffer[j].t_arr < s.buffer[at].t_arr) at = j;
                } else if (policy == lang::RemovalPolicy::LCFS) {
                    for (std::size_t j = 1; j < s.buffer.size(); ++j)
                        if (s.buffer[j].t_arr > s.buffer[at].t_arr) at = j;
                } else {
                    at = std::min(waiting - 1,
                                  static_cast<std::size_t>(uniform01(g_routing) *
                                                           static_cast<double>(waiting)));
                }
                const std::size_t rc = s.buffer[at].cls;
                s.buffer.erase(s.buffer.begin() + static_cast<std::ptrdiff_t>(at));
                if (s.sched != SchedStrategy::FSP)
                    std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
                return rc;
            }
            ps_advance(i);
            std::size_t at = 0;
            if (policy == lang::RemovalPolicy::FCFS) {
                for (std::size_t j = 1; j < s.ps_jobs.size(); ++j)
                    if (s.ps_jobs[j].t_arr < s.ps_jobs[at].t_arr) at = j;
            } else if (policy == lang::RemovalPolicy::LCFS) {
                for (std::size_t j = 1; j < s.ps_jobs.size(); ++j)
                    if (s.ps_jobs[j].t_arr > s.ps_jobs[at].t_arr) at = j;
            } else {
                at = std::min(inservice - 1,
                              static_cast<std::size_t>(uniform01(g_routing) *
                                                       static_cast<double>(inservice)));
            }
            const std::size_t rc = s.ps_jobs[at].cls;
            s.ps_jobs.erase(s.ps_jobs.begin() + static_cast<std::ptrdiff_t>(at));
            // The freed share is taken by whatever an LPS cap was holding out.
            if (s.lps_limit > 0 && !s.buffer.empty() && s.ps_jobs.size() < s.lps_limit) {
                Job nextjob = buffer_pop(i);
                PsJob pj;
                pj.cls = nextjob.cls;
                pj.priority = nextjob.priority;
                pj.t_arr = nextjob.t_arr;
                pj.t_sys = nextjob.t_sys;
                pj.total = nextjob.service;
                pj.remaining = nextjob.service;
                pj.parent = nextjob.parent;
                s.ps_jobs.push_back(pj);
            }
            ps_reschedule(i);
            return rc;
        }
 
        // A BUFFERED station. A blocked slot is NOT a candidate: the job it
        // holds has already completed here and is charged to the destination
        // it is queueing for, so removing it would decrement the wrong station.
        std::size_t inservice = 0;
        for (std::size_t sl = 0; sl < s.nservers; ++sl)
            if (s.server_busy[sl]) ++inservice;
        if (waiting + inservice == 0) return K;
        bool from_wait;
        if (policy == lang::RemovalPolicy::RANDOM) {
            const std::size_t pick = static_cast<std::size_t>(
                uniform01(g_routing) * static_cast<double>(waiting + inservice));
            from_wait = pick < waiting;
        } else {
            from_wait = waiting > 0;
        }
 
        if (from_wait) {
            std::size_t at = 0;
            if (policy == lang::RemovalPolicy::FCFS) {
                for (std::size_t j = 1; j < s.buffer.size(); ++j)
                    if (s.buffer[j].t_arr < s.buffer[at].t_arr) at = j;
            } else if (policy == lang::RemovalPolicy::LCFS) {
                for (std::size_t j = 1; j < s.buffer.size(); ++j)
                    if (s.buffer[j].t_arr > s.buffer[at].t_arr) at = j;
            } else {
                at = std::min(waiting - 1,
                              static_cast<std::size_t>(uniform01(g_routing) *
                                                       static_cast<double>(waiting)));
            }
            const std::size_t rc = s.buffer[at].cls;
            s.buffer.erase(s.buffer.begin() + static_cast<std::ptrdiff_t>(at));
            if (s.sched != SchedStrategy::FSP)
                std::make_heap(s.buffer.begin(), s.buffer.end(), s.cmp);
            return rc;
        }
 
        std::size_t slot = s.nservers;
        for (std::size_t sl = 0; sl < s.nservers; ++sl) {
            if (!s.server_busy[sl]) continue;
            if (slot == s.nservers) {
                slot = sl;
                continue;
            }
            if (policy == lang::RemovalPolicy::FCFS) {
                if (s.server[sl].t_arr < s.server[slot].t_arr) slot = sl;
            } else if (policy == lang::RemovalPolicy::LCFS) {
                if (s.server[sl].t_arr > s.server[slot].t_arr) slot = sl;
            }
        }
        if (policy == lang::RemovalPolicy::RANDOM) {
            std::size_t pick = std::min(inservice - 1,
                                        static_cast<std::size_t>(uniform01(g_routing) *
                                                                 static_cast<double>(inservice)));
            for (std::size_t sl = 0; sl < s.nservers; ++sl)
                if (s.server_busy[sl]) {
                    if (pick == 0) {
                        slot = sl;
                        break;
                    }
                    --pick;
                }
        }
        if (slot == s.nservers) return K;
        sd_advance(i);
        const std::size_t rc = s.server[slot].cls;
        // Freeing the slot and clearing its tag is what neutralises the
        // departure already in the list; there is no handle to cancel it with.
        s.server_busy[slot] = false;
        s.server_tag[slot] = 0;
        acc.update_busy(i, rc, now);
        acc.busy[i][rc] -= 1.0;
        // The vacated server takes the next job, under the station's own order:
        // the polling controller decides for itself, everyone else pops.
        if (s.polling) {
            poll_serve(i);
        } else if (!s.server_blocked[slot] && !s.server_held[slot] &&
                   buffer_has_for_slot(i, slot)) {
            Job nextjob = buffer_pop_for_slot(i, slot);
            start_service(i, slot, nextjob);
        }
        sd_reschedule(i);
        return rc;
    };
 
    /**
     * A G-NETWORK REMOVAL SIGNAL arriving at station `i`, and annihilated there.
     *
     * The signal NEVER JOINS the station and never continues along its routing
     * chain: it removes a batch of the jobs it finds, and one that finds the
     * station empty is simply lost. That is Gelenbe's semantics and the one
     * `State.afterEventStationSignal` and the Java engine both implement;
     * routing the signal onward instead would make a single signal fire once
     * per downstream station, which is what the tandem regression measures.
     *
     * A CATASTROPHE removes EVERY job held, in-service ones included, and
     * ignores the batch-size law entirely -- removing all of them is what the
     * word means. A NEGATIVE signal draws its batch size from `signalremdist`
     * (absent: exactly one) and CLIPS it at what the station holds, so an
     * oversized batch drains the station instead of driving it negative.
     */
    auto apply_signal = [&](std::size_t i, const Job& sig) {
        StationState& s = S[i];
        const std::size_t held = signal_held(i);
        if (held > 0) {
            std::size_t to_remove = 1;
            if (sig.cls < sn.signaltype.size() &&
                sn.signaltype[sig.cls] == lang::SignalType::CATASTROPHE) {
                to_remove = held;
            } else if (sig.cls < sn.signalremdist.size() && !sn.signalremdist[sig.cls].empty()) {
                const std::vector<T>& pmf = sn.signalremdist[sig.cls];
                const double u = uniform01(g_routing);
                double cum = 0.0;
                std::size_t b = 0;
                for (; b + 1 < pmf.size(); ++b) {
                    cum += static_cast<double>(pmf[b]);
                    if (u < cum) break;
                }
                to_remove = std::min(b, held);
            } else {
                to_remove = std::min<std::size_t>(1, held);
            }
            lang::RemovalPolicy policy = lang::RemovalPolicy::RANDOM;
            if (sig.cls < sn.signalrempolicy.size()) policy = sn.signalrempolicy[sig.cls];
            for (std::size_t rep = 0; rep < to_remove; ++rep) {
                if (signal_held(i) == 0) break;
                const std::size_t rc = signal_remove_one(i, policy);
                if (rc >= K) break;
                sd_advance(i);
                acc.update_qlen(i, rc, now);
                acc.qlen[i][rc] -= 1.0;
                bp.track(i, rc, -1, now);
                // A job destroyed inside a finite capacity region frees the
                // slot it occupied; leaving it charged would shrink the region
                // by one for the rest of the run.
                if (region_of[i] >= 0) {
                    const std::size_t rg = static_cast<std::size_t>(region_of[i]);
                    regions[rg].update(now);
                    regions[rg].leave(rc);
                    release_region(rg);
                }
            }
            release_blocked(i);
        }
        // The signal is annihilated here: it leaves the system without ever
        // having been at a station, so it is a system completion and nothing
        // else -- no arrival, no queue length and no station throughput.
        sys_resp_sum[sig.cls] += now - sig.t_sys;
        sys_resp_cnt[sig.cls] += 1.0;
        sys_completed[sig.cls] += 1.0;
    };
 
 
    /**
     * The completion budget, spent by anything the model actually DOES.
     *
     * Declared ahead of `deliver` because a service completion is not the only
     * thing that spends it: a Petri net has no services of its own and counts a
     * FIRING, and a cache-only model (Source -> Cache -> Sink) has none either
     * and counts the cache's delivery to the Sink. Without the latter the loop
     * at the bottom of this function never terminates on such a model, since
     * the Source keeps the event queue non-empty forever. `Solver_ssj` spends
     * the budget at the same point (`deliverCacheCompletion` calls
     * `checkEventCountStop` on the sink branch), so this is the reference's
     * rule and not a local convention.
     */
    std::uint64_t total_completions = 0;
 
    /**
     * Free every request parked while `item` was being fetched.
     *
     * Declared ahead of `deliver` and assigned after it because the two are
     * mutually recursive: a release routes each freed request onward, and a
     * delivery to a cache is what triggers a release. Each freed request reads
     * the now-cached item and completes at once, counted as a DELAYED HIT and
     * leaving under its own hit class; its original `t_sys` is preserved so the
     * response time still spans the wait it actually served.
     */
    std::function<void(CacheState&, std::size_t, std::size_t)> release_delayed_hits;
 
    /**
     * Hand a job to a NODE, replicating at a Fork and synchronizing at a Join.
     *
     * A FORK sends a sibling down EVERY outgoing edge, not one drawn from the
     * routing probabilities: the edges of a fork are branches taken together,
     * and drawing among them would turn a fork into a probabilistic split with
     * the same picture and a completely different response time. `tasksPerLink`
     * multiplies each branch.
     *
     * A JOIN releases ONE job once its strategy is satisfied and DISCARDS the
     * siblings that arrive afterwards -- under a quorum they are the losers of
     * the race, and counting them as completions would inflate the throughput
     * of everything downstream.
     */
    std::function<void(std::size_t, Job, std::size_t)> deliver =
        [&](std::size_t node, Job job, std::size_t from) {
        const NodeType nt = sn.nodes[node - 1].nodetype;
 
        if (nt == NodeType::Sink) {
            sys_resp_sum[job.cls] += now - job.t_sys;
            sys_resp_cnt[job.cls] += 1.0;
            sys_completed[job.cls] += 1.0;
            return;
        }
 
        if (nt == NodeType::Cache) {
            auto ci = caches.find(node);
            if (ci == caches.end())
                throw InputError("SolverLDES (native engine): Cache node '" +
                                 sn.nodes[node - 1].name + "' carries no parameters");
            CacheState& cs = ci->second;
            Job out = job;
 
            // A RETURNING FETCH, not a read. The job wearing a retrieval class
            // is the fetch this cache sent out on an earlier miss: it inserts
            // the item, clears the in-flight flag and frees everyone parked
            // behind it. The miss was counted when the fetch was TRIGGERED, so
            // counting it again here would double every miss.
            if (cs.has_retrieval) {
                const int fetched = cs.fetches_item(job.cls);
                if (fetched >= 0) {
                    const std::size_t fi = static_cast<std::size_t>(fetched);
                    cache_miss(cs, fi, uniform01(g_routing), uniform01(g_routing));
                    cs.in_flight[fi] = 0;
                    cs.total_fetch_time += now - cs.fetch_start[fi];
                    cs.completed_fetches += 1.0;
                    release_delayed_hits(cs, fi, node);
                    if (cs.miss_class[job.cls] >= 0)
                        out.cls = static_cast<std::size_t>(cs.miss_class[job.cls]);
                    const RouteEntry re = draw_node_route(node, out.cls);
                    out.cls = re.cls;
                    if (re.node != 0 && sn.nodes[re.node - 1].nodetype == NodeType::Sink)
                        ++total_completions;
                    deliver(re.node, out, M);
                    return;
                }
            }
 
            const std::vector<double>& pop = cs.popularity[job.cls];
            if (pop.empty())
                throw InputError("SolverLDES (native engine): class '" +
                                 sn.classes[job.cls].name + "' reads cache '" +
                                 sn.nodes[node - 1].name + "' with no item popularity");
            const double u = uniform01(g_routing);
            std::size_t item = 0;
            while (item + 1 < pop.size() && u >= pop[item]) ++item;
 
            const int at = cs.find(item);
            if (at >= 0) {
                cs.hits[job.cls] += 1.0;
                cache_hit(cs, item, static_cast<std::size_t>(at), uniform01(g_routing));
                if (cs.hit_class[job.cls] >= 0)
                    out.cls = static_cast<std::size_t>(cs.hit_class[job.cls]);
            } else if (cs.has_retrieval && cs.retrieval_of(item, job.cls) > 0) {
                if (cs.in_flight[item]) {
                    // A fetch for this item is already out: PARK the request. It
                    // leaves the event stream entirely and is released when the
                    // fetch returns, so it is neither a hit nor a miss yet.
                    CacheState::HeldRequest hr;
                    hr.cls = job.cls;
                    hr.hold_time = now;
                    hr.t_sys = job.t_sys;
                    cs.held[item].push_back(hr);
                    return;
                }
                // Trigger the fetch: count the miss now and leave wearing the
                // per-item retrieval class, which the routing carries through
                // the retrieval stations and back to this cache.
                cs.misses[job.cls] += 1.0;
                cs.in_flight[item] = 1;
                cs.fetch_start[item] = now;
                out.cls = cs.retrieval_of(item, job.cls) - 1;
            } else {
                cs.misses[job.cls] += 1.0;
                cache_miss(cs, item, uniform01(g_routing), uniform01(g_routing));
                if (cs.miss_class[job.cls] >= 0)
                    out.cls = static_cast<std::size_t>(cs.miss_class[job.cls]);
            }
            // The hit and miss classes are what the POST_CACHE branch routes on,
            // so the job leaves the cache under its new class.
            const RouteEntry e = draw_node_route(node, out.cls);
            out.cls = e.cls;
            // A cache delivering to a Sink is this model's only completion when
            // it has no service station, so it spends the budget (see above).
            if (e.node != 0 && sn.nodes[e.node - 1].nodetype == NodeType::Sink)
                ++total_completions;
            deliver(e.node, out, M);
            return;
        }
 
        if (nt == NodeType::ClassSwitch || nt == NodeType::Router || nt == NodeType::Logger) {
            // A PASS-THROUGH NODE HOLDS NOTHING. It changes the class (a
            // ClassSwitch), picks a destination (a Router) or records the
            // crossing (a Logger) and the job continues in the SAME instant, so
            // there is no queue, no service and no row to fill.
            //
            // `link()` synthesizes `CS_<i>_to_<j>` for every link that switches
            // class, so this branch carries every class switch the routing
            // matrix used to hold as an edge attribute. Such a node has ONE
            // outgoing destination, so `draw_node_route` spends no draw on it
            // and the sample path stays in step with the Java engine.
            //
            // The walk is iterative, not recursive, and bounded: a routing that
            // sends a job round a cycle of pass-through nodes never reaches a
            // station and would otherwise hang the run.
            Job out = job;
            std::size_t at = node;
            for (std::size_t hop = 0; hop <= nnodes; ++hop) {
                const RouteEntry e = draw_node_route(at, out.cls);
                out.cls = e.cls;
                if (e.node == 0) return;
                const NodeType dt = sn.nodes[e.node - 1].nodetype;
                if (dt == NodeType::ClassSwitch || dt == NodeType::Router ||
                    dt == NodeType::Logger) {
                    at = e.node;
                    continue;
                }
                deliver(e.node, out, from);
                return;
            }
            throw InputError("SolverLDES (native engine): the routing out of node '" +
                             sn.nodes[node - 1].name +
                             "' cycles through pass-through nodes without reaching a station");
        }
 
        if (nt == NodeType::Fork) {
            // A FORK TAKES EVERY EDGE, so it flattens the two-stage table back
            // into one sibling per (destination, arriving class): there is no
            // dispatcher decision to make and no draw to spend.
            const std::vector<RouteDest>& tab = nroute[node][job.cls];
            if (tab.empty()) return;
            std::vector<RouteEntry> links;
            for (std::size_t di = 0; di < tab.size(); ++di)
                for (std::size_t ci = 0; ci < tab[di].cls.size(); ++ci) {
                    RouteEntry e;
                    e.node = tab[di].node;
                    e.station = tab[di].station;
                    e.sink = tab[di].sink;
                    e.cls = tab[di].cls[ci].first;
                    links.push_back(e);
                }
            const qn::NodeDef& fk = sn.nodes[node - 1];
            const qn::ForkParam<double>* fp = sn.fork_param_of(node);
            const int per_link = std::max(1, static_cast<int>(fk.tasks_per_link + 0.5));
 
            // How many tasks each branch emits on THIS firing. A plain fork
            // takes `per_link` on every link and draws NOTHING -- which is what
            // keeps every seeded golden and the bit-exact agreement with the
            // Java engine intact. A variable one draws the branch activation
            // and the degree per link, so the count the Join waits for is known
            // only after the draws.
            std::vector<int> per_branch(links.size(), per_link);
            int total = 0;
            const bool variable = fork_is_variable(node);
            for (std::size_t li = 0; li < links.size(); ++li) {
                if (variable) {
                    const std::size_t k0 = links[li].node - 1, r0 = job.cls;
                    const double p = fp->fan_out_prob(k0, r0);
                    if (p < 1.0 && g_fork.aux.next_double() >= p) {
                        per_branch[li] = 0;
                    } else if (!fp->fan_out_dist[k0][r0].disabled) {
                        per_branch[li] = sample_fork_degree(fp->fan_out_dist[k0][r0]);
                    } else {
                        per_branch[li] = static_cast<int>(fp->fan_out_link(k0, r0) + 0.5);
                    }
                }
                total += per_branch[li];
            }
            if (total <= 0)
                throw InputError("SolverLDES (native engine): Fork '" + fk.name +
                                 "' emitted no sibling; at least one branch must be certain "
                                 "to emit at least one task");
 
            ForkSync fs;
            fs.cls = job.cls;
            fs.t_sys = job.t_sys;
            fs.total = total;
            fs.required = fs.total;
            const int jn = join_of_fork[node];
            if (jn > 0) {
                auto jd = sn.joindecl.find(static_cast<std::size_t>(jn));
                if (jd != sn.joindecl.end() &&
                    jd->second.strategy != lang::JoinStrategy::STD && jd->second.quorum > 0.0) {
                    // The quorum is a count of THIS fork's siblings, fixed here
                    // rather than at the Join: a Join shared by two forks of
                    // different fan-out would otherwise apply one fork's count
                    // to the other's siblings.
                    fs.required = std::min(fs.total,
                                           std::max(1, static_cast<int>(jd->second.quorum + 0.5)));
                }
            }
            const std::uint64_t pid = ++next_parent;
            fork_sync[pid] = fs;
            for (std::size_t li = 0; li < links.size(); ++li)
                for (int t = 0; t < per_branch[li]; ++t) {
                    Job sib;
                    sib.cls = links[li].cls;
                    sib.t_sys = job.t_sys;
                    sib.parent = pid;
                    deliver(links[li].node, sib, M);
                }
            return;
        }
 
        if (nt == NodeType::Join) {
            auto it = fork_sync.find(job.parent);
            if (it == fork_sync.end()) {
                // Its parent already synchronized: this sibling lost the quorum
                // race and is discarded rather than released a second time. It
                // goes into the SAME accumulators every other station uses, so
                // it is warmup-truncated with them; see Accum::join_dropped.
                //
                // IT STILL ARRIVED. `arrived` is the OFFERED rate, and the whole
                // point of a quorum is that N siblings are offered and only K
                // consumed, so leaving the straggler out reports AN = K*TN and
                // hides the very loss the measurement exists for. The reference
                // counts it the same way: its AN at the join is exactly N times
                // the fork rate under a quorum, not K times.
                const std::size_t jd = node_to_station[node];
                if (jd < M) {
                    acc.arrived[jd][job.cls] += 1.0;
                    acc.join_dropped[jd][job.cls] += 1.0;
                }
                return;
            }
            ForkSync& fs = it->second;
            // A JOIN IS A STATION LIKE ANY OTHER. Its siblings are its queue
            // length and its synchronizations are its completions, so they go
            // into the same accumulators every other station uses and are
            // warmup-truncated with them rather than measured on a private
            // tally the result assembly never reads.
            // `js == M` is a Join the refresh registered as a node and not as a
            // station, which has no row to fill; `handleJoinArrival` gates on
            // the same bound. The synchronization itself runs either way.
            const std::size_t js = node_to_station[node];
            if (js < M) {
                acc.update_qlen(js, job.cls, now);
                acc.qlen[js][job.cls] += 1.0;
                acc.arrived[js][job.cls] += 1.0;
            }
            fs.siblings.push_back(std::make_pair(job.cls, now));
            if (static_cast<int>(fs.siblings.size()) < fs.required) return;
 
            // Satisfied: the siblings leave the Join together and one job goes on.
            if (js < M) {
                for (std::size_t si = 0; si < fs.siblings.size(); ++si) {
                    const std::size_t sc = fs.siblings[si].first;
                    acc.update_qlen(js, sc, now);
                    acc.qlen[js][sc] -= 1.0;
                    // Each sibling waited from ITS OWN arrival, not the first:
                    // the Join charges the delay it imposed on each.
                    acc.resp_sum[js][sc] += now - fs.siblings[si].second;
                    acc.resp_cnt[js][sc] += 1.0;
                }
                acc.completed[js][fs.cls] += 1.0;
            }
            Job merged;
            merged.cls = fs.cls;
            merged.t_sys = fs.t_sys;
            fork_sync.erase(it);
            const RouteEntry e = draw_node_route(node, merged.cls);
            merged.cls = e.cls;
            deliver(e.node, merged, M);
            return;
        }
 
        // An ordinary station.
        const std::size_t j = node_to_station[node];
        // A REMOVAL SIGNAL acts on the station and is annihilated: it is
        // intercepted here, before any of the accounting an arrival gets,
        // because it never becomes one.
        if (has_removal_signal && j < M && is_removal_signal[job.cls]) {
            apply_signal(j, job);
            return;
        }
        /**
         * A REPLY ANSWERS ITS CALL AND TRAVELS ON. It releases the server the
         * caller has been holding, lets that station take its next waiting job,
         * and then CONTINUES along its own routing -- which is what separates it
         * from a removal signal, annihilated where it lands.
         *
         * A reply with no matching call is not an error: the call may have been
         * answered already, or the class may be routed as a reply without any
         * caller having parked. It routes on, having released nothing.
         */
        if (has_sync_call && j < M && is_reply_signal[job.cls]) {
            acc.completed[j][job.cls] += 1.0;
            typename std::map<std::uint64_t, PendingCall>::iterator pc =
                pending_reply.find(job.call);
            if (pc != pending_reply.end()) {
                const std::size_t bi = pc->second.station, bs = pc->second.slot;
                const std::size_t bc = pc->second.cls;
                pending_reply.erase(pc);
                StationState& bst = S[bi];
                bst.server_held[bs] = false;
                acc.update_busy(bi, bc, now);
                if (acc.busy[bi][bc] > 0.0) acc.busy[bi][bc] -= 1.0;
                if (!bst.server_blocked[bs] && buffer_has_for_slot(bi, bs)) {
                    Job nextjob = buffer_pop_for_slot(bi, bs);
                    start_service(bi, bs, nextjob);
                    sd_reschedule(bi);
                }
                release_blocked(bi);
            }
            const RouteEntry re2 = draw_node_route(node, job.cls);
            Job onward;
            onward.cls = re2.cls;
            onward.t_sys = job.t_sys;
            onward.parent = job.parent;
            deliver(re2.node, onward, M);
            return;
        }
        Job moved = job;
        // A closed class completes its passage on reaching its reference
        // station: that crossing IS the system completion, and the elapsed time
        // is the cycle time.
        if (j < M && sn.classes[job.cls].type != lang::JobClassType::OPEN &&
            j + 1 == sn.classes[job.cls].refstat) {
            sys_resp_sum[job.cls] += now - job.t_sys;
            sys_resp_cnt[job.cls] += 1.0;
            sys_completed[job.cls] += 1.0;
            moved.t_sys = now;
        }
        if (!admit(j, moved, from) && sn.classes[job.cls].type != lang::JobClassType::OPEN)
            throw InputError("SolverLDES (native engine): a closed job was dropped at "
                             "station '" + sn.stations[j].name +
                             "'; a closed class cannot lose population");
    };
 
    release_delayed_hits = [&](CacheState& cs, std::size_t item, std::size_t node) {
        if (item >= cs.held.size() || cs.held[item].empty()) return;
        std::vector<CacheState::HeldRequest> freed;
        freed.swap(cs.held[item]);
        for (std::size_t i = 0; i < freed.size(); ++i) {
            const CacheState::HeldRequest& hr = freed[i];
            cs.delayed[hr.cls] += 1.0;
            cs.delayed_wait += now - hr.hold_time;
            Job out;
            out.cls = (cs.hit_class[hr.cls] >= 0) ? static_cast<std::size_t>(cs.hit_class[hr.cls])
                                                  : hr.cls;
            out.t_sys = hr.t_sys;
            const RouteEntry e = draw_node_route(node, out.cls);
            out.cls = e.cls;
            if (e.node != 0 && sn.nodes[e.node - 1].nodetype == NodeType::Sink)
                ++total_completions;
            deliver(e.node, out, M);
        }
    };
 
 
    /**
     * The marking as one flat vector over (place slot, class), `p * K + r`,
     * which is what a mode's arcs are indexed by.
     *
     * SUMMED OVER CLASSES UNTIL 2026-08-12, which made every arc colourless and
     * let a token of one class satisfy another's pre-arc.
     */
    auto flat_marking = [&]() {
        std::vector<double> tok(place_nodes.size() * K, 0.0);
        for (std::size_t p = 0; p < place_nodes.size(); ++p)
            for (std::size_t r = 0; r < K; ++r) tok[p * K + r] = marking[p][r];
        return tok;
    };
 
    /**
     * Fire every enabled IMMEDIATE mode, then arm the timed ones.
     *
     * The immediate modes are exhausted FIRST and to a fixed point: a marking
     * that enables one is not a marking the net rests in, so letting a timed
     * mode fire from it would visit a state the model does not have. The guard
     * bounds a net whose immediate modes form a cycle, which is a modelling
     * error rather than something to simulate forever.
     */
    std::uint64_t fire_tag = 0;
    std::vector<std::uint64_t> live_fire(transitions.size(), 0);
    std::function<void()> spn_settle = [&]() {
        if (transitions.empty()) return;
        for (int guard = 0; guard < 1000000; ++guard) {
            bool fired_any = false;
            std::vector<double> tok = flat_marking();
            for (std::size_t t = 0; t < transitions.size(); ++t) {
                const int k = spn_pick_immediate(transitions[t].modes, tok, uniform01(g_routing));
                if (k < 0) continue;
                const SpnMode& m = transitions[t].modes[static_cast<std::size_t>(k)];
                // The tokens move in THEIR OWN CLASS: an arc names a (place,
                // class) pair, so a Class2 pre-arc consumes Class2 tokens.
                for (std::size_t p = 0; p < place_nodes.size(); ++p)
                    for (std::size_t r = 0; r < K; ++r) {
                        marking[p][r] -= m.enabling[p * K + r];
                        marking[p][r] += m.firing[p * K + r];
                    }
                transitions[t].fired[static_cast<std::size_t>(k)] += 1.0;
                ++total_completions;
                fired_any = true;
                break;
            }
            if (!fired_any) break;
        }
        // Arm one clock per transition on its fastest enabled timed mode. The
        // tag invalidates it as soon as the marking moves, because a
        // marking-dependent rate must be RESAMPLED and not merely rescaled.
        const std::vector<double> tok = flat_marking();
        for (std::size_t t = 0; t < transitions.size(); ++t) {
            double best = std::numeric_limits<double>::infinity();
            int best_mode = -1;
            for (std::size_t k = 0; k < transitions[t].modes.size(); ++k) {
                const SpnMode& m = transitions[t].modes[k];
                if (m.immediate || !spn_enabled(m, tok)) continue;
                const double rate = spn_rate(m, tok);
                if (!(rate > 0.0)) continue;
                const double d = -std::log(uniform01(g_spn)) / rate;
                if (d < best) {
                    best = d;
                    best_mode = static_cast<int>(k);
                }
            }
            live_fire[t] = ++fire_tag;
            if (best_mode < 0) continue;
            Event e;
            e.t = now + best;
            e.kind = EV_FIRING;
            e.station = t;
            e.cls = static_cast<std::size_t>(best_mode);
            e.tag = live_fire[t];
            push(e);
        }
    };
 
    // ---- initial state -----------------------------------------------------
    /**
     * WARM START, `--initsol`: the placement a companion solver's steady state
     * implies, as a STATION-MAJOR vector [st0_cl0, ..., stM-1_clK-1].
     *
     * It REPLACES the reference-station placement of the closed classes; the
     * open arrival streams are armed either way. A Petri net is exempt: its
     * tokens live in Places and not in the service queues this vector
     * describes, so a station-based placement would put nothing anywhere and
     * then fail the conservation check.
     *
     * CLASSES DESCEND WITHIN A STATION. Every initial job arrives at time 0, so
     * under an order-preserving discipline the LAST one injected holds the
     * server; ascending order puts the highest-index class in service, which is
     * the mirror of what `State.initDefault` encodes and selects the other
     * closed communicating class on a chain made reducible by non-overtaking
     * routing -- a different stationary distribution, not a different tie-break.
     */
    const bool warm_start = !o.init_sol.empty() && place_nodes.empty();
    if (warm_start) {
        for (std::size_t i = 0; i < M; ++i) {
            if (S[i].role != Role::Queue && S[i].role != Role::Delay) continue;
            for (std::size_t k = K; k-- > 0;) {
                const std::size_t idx = i * K + k;
                if (idx >= o.init_sol.size()) continue;
                const double v = o.init_sol[idx];
                if (!(v > 0.0)) continue;
                const std::size_t count = static_cast<std::size_t>(v);
                for (std::size_t j = 0; j < count; ++j) {
                    Job job;
                    job.cls = k;
                    job.t_sys = 0.0;
                    if (!admit(i, job, M))
                        throw InputError("SolverLDES (native engine): the warm-start placement "
                                         "of class '" + sn.classes[k].name + "' does not fit "
                                         "station '" + sn.stations[i].name + "'");
                }
            }
        }
        // A closed class must be placed in full: a warm start that loses jobs is
        // a different model, and the loss is invisible in every mean it reports.
        for (std::size_t r = 0; r < K; ++r) {
            if (sn.classes[r].type == lang::JobClassType::OPEN) continue;
            double held = 0.0;
            for (std::size_t i = 0; i < M; ++i) held += acc.qlen[i][r];
            const double want = sn.classes[r].population;
            if (std::fabs(held - want) > 0.5)
                throw InputError("SolverLDES (native engine): the warm-start placement holds " +
                                 std::to_string(static_cast<long>(held + 0.5)) + " jobs of class '" +
                                 sn.classes[r].name + "' against a population of " +
                                 std::to_string(static_cast<long>(want + 0.5)));
        }
    }
    // Every breakdown-carrying station starts UP with its first failure armed,
    // which is the reference's initial condition: a model that started down
    // would report a transient nobody asked for.
    for (std::size_t i = 0; i < M; ++i) {
        if (!S[i].has_breakdown) continue;
        Event e;
        e.t = S[i].failure_time.next_at(g_aux[i], 0.0);
        e.kind = EV_BREAKDOWN;
        e.station = i;
        push(e);
    }
 
    std::vector<Sampler> arrival(K);
    std::vector<double> lambda(K, 0.0);
    for (std::size_t r = 0; r < K; ++r) {
        if (sn.classes[r].type == lang::JobClassType::OPEN) {
            if (source_st >= M)
                throw InputError("SolverLDES (native engine): an open class with no Source");
            if (S[source_st].off[r]) continue;
            arrival[r] = S[source_st].svc[r];
            lambda[r] = 1.0 / arrival[r].mean();
            Event e;
            e.t = slot_snap(arrival[r].next_at(g_arr[r], 0.0), "interarrival time");
            e.kind = EV_ARRIVAL;
            e.station = source_st;
            e.cls = r;
            push(e);
        } else {
            if (warm_start) continue;  // already placed from the warm start
            const double n = sn.classes[r].population;
            if (!(n >= 0.0) || !std::isfinite(n))
                throw InputError("SolverLDES (native engine): class '" + sn.classes[r].name +
                                 "' has a population that is not a finite count");
            const std::size_t refst = sn.classes[r].refstat;
            if (refst == 0 || refst > M)
                throw InputError("SolverLDES (native engine): class '" + sn.classes[r].name +
                                 "' has no reference station");
            const std::size_t count = static_cast<std::size_t>(n + 0.5);
            for (std::size_t j = 0; j < count; ++j) {
                Job job;
                job.cls = r;
                job.t_sys = 0.0;
                if (!admit(refst - 1, job, M))
                    throw InputError("SolverLDES (native engine): the initial population of "
                                     "class '" + sn.classes[r].name +
                                     "' does not fit its reference station's capacity");
            }
        }
    }
 
    // ---- convergence-based stopping -----------------------------------------
    Convergence cnvg;
    cnvg.init(M, K, o, max_events);
    std::vector<std::size_t> servers_of(M, 1);
    std::vector<std::vector<bool>> off_of(M, std::vector<bool>(K, true));
    for (std::size_t i = 0; i < M; ++i) {
        servers_of[i] = S[i].nservers;
        for (std::size_t r = 0; r < K; ++r)
            off_of[i][r] = S[i].off[r] || S[i].role == Role::Source ||
                           S[i].role == Role::Synchronization;
    }
    std::uint64_t last_cnvg_events = 0;
    bool converged = false;
 
    // The immediate modes settle the initial marking before any clock runs.
    spn_settle();
 
    // ---- MSER-5 and batch-means bookkeeping --------------------------------
    Observations obs(M, K, o.tranfilter == "mser5", o.mserbatch > 0 ? o.mserbatch : 5);
    for (std::size_t i = 0; i < M; ++i)
        if (S[i].role == Role::Source || S[i].role == Role::Synchronization) obs.in_mser[i] = 0;
    std::uint64_t mser_interval = max_events / 1000ULL;
    if (mser_interval < 1) mser_interval = 1;
    std::uint64_t last_mser_events = 0;
 
 
    // ---- transient sampling and the joint-state histogram ---------------------
    /**
     * A TRANSIENT run is bounded by SIMULATED TIME, not by completions.
     *
     * `options.timespan = [t0, t1]` sets the horizon and the completion budget
     * is then ignored, exactly as the reference does: a transient trajectory is
     * a function of time, so stopping on an event count would end it at a
     * different instant on every path and make the series incomparable across
     * seeds.
     */
    const bool transient_run = o.has_timespan && std::isfinite(o.t1) && o.t1 > o.t0;
    const double horizon = transient_run ? o.t1 : std::numeric_limits<double>::infinity();
    // ~1000 points, the reference's target: enough to draw and few enough that
    // the series is not itself a memory problem on a long horizon.
    const double tran_interval = transient_run ? (o.t1 - o.t0) / 1000.0 : 0.0;
    double next_tran_sample = transient_run ? o.t0 + tran_interval : 0.0;
    double last_tran_time = transient_run ? o.t0 : 0.0;
    std::vector<double> tran_times;
    std::vector<std::vector<std::vector<double>>> tran_q(M, std::vector<std::vector<double>>(K));
    std::vector<std::vector<std::vector<double>>> tran_u(M, std::vector<std::vector<double>>(K));
    std::vector<std::vector<std::vector<double>>> tran_t(M, std::vector<std::vector<double>>(K));
    std::vector<std::vector<double>> last_tran_q(M, std::vector<double>(K, 0.0));
    std::vector<std::vector<double>> last_tran_b(M, std::vector<double>(K, 0.0));
    std::vector<std::vector<double>> last_tran_c(M, std::vector<double>(K, 0.0));
 
    /**
     * The EXACT joint-state residence time, keyed by the aggregate state row.
     *
     * `histogram_space` plus `histogram_time` let a caller evaluate its own
     * reward on the empirical distribution, INCLUDING a nonlinear one:
     * E[r] = sum_s (t_s / sum t) r(state_s) is exact only because the residence
     * time of each distinct state is kept, which a trajectory of means cannot
     * reconstruct.
     */
    std::map<std::vector<int>, double> histogram;
    std::vector<std::pair<double, std::vector<int>>> trajectory;
    double hist_last = 0.0;
    // A transient run needs the state walk anyway; a steady-state one pays for
    // it only when asked.
    const bool want_traj = transient_run || o.export_trajectory;
    const bool want_hist = transient_run || o.export_histogram || want_traj;
 
    auto joint_state = [&]() {
        std::vector<int> row(M * K, 0);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r)
                row[i * K + r] = static_cast<int>(acc.qlen[i][r] + 0.5);
        return row;
    };
    auto hist_accumulate = [&]() {
        if (!want_hist) return;
        const double dt = now - hist_last;
        if (dt > 0.0) {
            const std::vector<int> row = joint_state();
            histogram[row] += dt;
            if (want_traj) trajectory.push_back(std::make_pair(hist_last, row));
        }
        hist_last = now;
    };
 
    auto tran_sample = [&]() {
        const double dt = now - last_tran_time;
        tran_times.push_back(now);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                if (dt > 0.0) {
                    tran_q[i][r].push_back((acc.tot_qlen[i][r] - last_tran_q[i][r]) / dt);
                    const double c = S[i].util_peak;
                    tran_u[i][r].push_back(
                        (S[i].role == Role::Delay)
                            ? (acc.tot_qlen[i][r] - last_tran_q[i][r]) / dt
                            : (acc.tot_busy[i][r] - last_tran_b[i][r]) / (dt * c));
                    tran_t[i][r].push_back((acc.completed[i][r] - last_tran_c[i][r]) / dt);
                } else {
                    tran_q[i][r].push_back(acc.qlen[i][r]);
                    tran_u[i][r].push_back(0.0);
                    tran_t[i][r].push_back(0.0);
                }
                last_tran_q[i][r] = acc.tot_qlen[i][r];
                last_tran_b[i][r] = acc.tot_busy[i][r];
                last_tran_c[i][r] = acc.completed[i][r];
            }
        last_tran_time = now;
    };
 
    // ---- the event loop ----------------------------------------------------
    while (!evq.empty() && (transient_run || total_completions < max_events)) {
        const Event ev = evq.top();
        if (transient_run && ev.t > horizon) break;
        evq.pop();
        now = ev.t;
 
        // THE STATE IS WALKED BEFORE THE EVENT IS APPLIED, so the interval just
        // ended is charged to the state that held during it. A steady-state run
        // asked for the histogram takes the same walk and nothing else.
        if (want_hist) hist_accumulate();
        if (transient_run) {
            while (now >= next_tran_sample && next_tran_sample <= horizon) {
                const double save = now;
                now = next_tran_sample;
                for (std::size_t a = 0; a < M; ++a) {
                    if (S[a].ps) ps_advance(a);
                    for (std::size_t r = 0; r < K; ++r) {
                        acc.update_qlen(a, r, now);
                        acc.update_busy(a, r, now);
                    }
                }
                tran_sample();
                next_tran_sample += tran_interval;
                now = save;
            }
        }
 
        if (ev.kind == EV_BREAKDOWN) {
            // FREEZE AT THE OLD RATE, FLIP, RE-TIME AT THE NEW, in that order,
            // for the same reason sd_advance runs before a population change:
            // the rate in force over the elapsed interval is the one the OLD
            // state implied. The stale departures are neutralised by the tag.
            StationState& s = S[ev.station];
            if (s.ps)
                ps_advance(ev.station);
            else
                sd_advance(ev.station);
            s.up = !s.up;
            if (s.ps)
                ps_reschedule(ev.station);
            else
                sd_reschedule(ev.station);
            // The next transition is drawn from the clock now in force. A job
            // stalled by the outage has no departure event at all; the repair's
            // reschedule is what gives it one back, with its residual intact.
            // `next_at` returns a DURATION, not an instant: it takes `from`
            // only so a time-inhomogeneous family can read its own schedule.
            // Assigning it straight to `nxt.t` puts the next flip in the PAST,
            // which the loop pops immediately and re-arms -- an infinite
            // cascade at one simulated instant.
            Event nxt;
            nxt.t = now + (s.up ? s.failure_time.next_at(g_aux[ev.station], now)
                                : s.repair_time.next_at(g_aux[ev.station], now));
            nxt.kind = EV_BREAKDOWN;
            nxt.station = ev.station;
            push(nxt);
            continue;
        }
 
        if (ev.kind == EV_FIRING) {
            // A stale clock: the marking moved after this one was armed, so its
            // rate no longer describes the net.
            if (ev.station >= transitions.size() || live_fire[ev.station] != ev.tag) continue;
            SpnTransition& tr = transitions[ev.station];
            const std::vector<double> tok = flat_marking();
            const SpnMode& m = tr.modes[ev.cls];
            if (!spn_enabled(m, tok)) continue;
            for (std::size_t p = 0; p < place_nodes.size(); ++p)
                for (std::size_t r = 0; r < K; ++r) {
                    marking[p][r] -= m.enabling[p * K + r];
                    marking[p][r] += m.firing[p * K + r];
                }
            tr.fired[ev.cls] += 1.0;
            ++total_completions;
            spn_settle();
            continue;
        }
 
        if (ev.kind == EV_SWITCHOVER) {
            StationState& sp = S[ev.station];
            sp.poll_switching = false;
            sp.poll_at = ev.cls;
            sp.poll_budget = (sp.poll_type == lang::PollingType::EXHAUSTIVE)
                                 ? std::numeric_limits<std::size_t>::max()
                                 : ((sp.poll_type == lang::PollingType::KLIMITED) ? sp.poll_k : 1);
            bool work = false;
            for (const Job& j : sp.buffer)
                if (j.cls == sp.poll_at) {
                    work = true;
                    break;
                }
            if (work)
                poll_serve(ev.station);
            else
                poll_advance(ev.station);
            continue;
        }
 
        if (ev.kind == EV_RETRIAL) {
            StationState& sr = S[ev.station];
            const double dt = now - sr.orbit_last;
            if (dt > 0.0)
                for (std::size_t r = 0; r < K; ++r) sr.tot_orbit[r] += sr.orbit_size[r] * dt;
            sr.orbit_last = now;
            sr.orbit_size[ev.cls] -= 1.0;
            sr.retried[ev.cls] += 1.0;
            // The retry is an ordinary admission: it succeeds if the station has
            // room by then and rejoins the orbit if it does not.
            admit(ev.station, ev.job, M);
            continue;
        }
 
        if (ev.kind == EV_SETUP) {
            StationState& ss = S[ev.station];
            if (ev.slot == 1) {
                // A power-down. It is CANCELLED if a job arrived meanwhile: the
                // server never went cold, so that job must not pay a setup.
                double held = 0.0;
                for (std::size_t r = 0; r < K; ++r) held += acc.qlen[ev.station][r];
                if (held > 0.0 || now + 1e-12 < ss.delayoff_at) continue;
                ss.setup_on = false;
                ss.delayoff_at = std::numeric_limits<double>::infinity();
                continue;
            }
            if (ss.setup_on) continue;  // superseded
            ss.setup_running = false;
            ss.setup_on = true;
            ss.delayoff_at = std::numeric_limits<double>::infinity();
            for (std::size_t sl = 0; sl < ss.nservers && !ss.buffer.empty(); ++sl)
                if (!ss.server_busy[sl] && !ss.server_blocked[sl] && !ss.server_held[sl] &&
                    buffer_has_for_slot(ev.station, sl)) {
                    Job nextjob = buffer_pop_for_slot(ev.station, sl);
                    start_service(ev.station, sl, nextjob);
                }
            sd_reschedule(ev.station);
            continue;
        }
 
        if (ev.kind == EV_RENEGE) {
            // The job abandons only if it is STILL WAITING. Having entered
            // service, departed, or been dropped all leave the timer with
            // nothing to find, which is the ordinary case and not an error.
            StationState& sr = S[ev.station];
            std::size_t at = sr.buffer.size();
            for (std::size_t j = 0; j < sr.buffer.size(); ++j)
                if (sr.buffer[j].id == ev.tag) {
                    at = j;
                    break;
                }
            if (at == sr.buffer.size()) continue;
            const std::size_t rc = sr.buffer[at].cls;
            sd_advance(ev.station);
            acc.update_qlen(ev.station, rc, now);
            acc.qlen[ev.station][rc] -= 1.0;
            bp.track(ev.station, rc, -1, now);
            sr.buffer.erase(sr.buffer.begin() + static_cast<std::ptrdiff_t>(at));
            if (sr.sched != SchedStrategy::FSP)
                std::make_heap(sr.buffer.begin(), sr.buffer.end(), sr.cmp);
            reneged[ev.station][rc] += 1.0;
            sd_reschedule(ev.station);
            continue;
        }
 
        if (ev.kind == EV_ARRIVAL) {
            Event nxt;
            const double gap = slot_snap(arrival[ev.cls].next_at(g_arr[ev.cls], now),
                                        "interarrival time");
            // A non-cyclic schedule past its horizon yields NO further arrival:
            // the stream has ended, and re-arming it at rate zero would loop.
            if (!(gap > 0.0) && arrival[ev.cls].time_varying()) continue;
            nxt.t = now + gap;
            nxt.kind = EV_ARRIVAL;
            nxt.station = ev.station;
            nxt.cls = ev.cls;
            push(nxt);
 
            const RouteEntry e = draw_node_route(sn.station_to_node[ev.station], ev.cls);
            Job job;
            job.cls = e.cls;
            job.t_sys = now;
            deliver(e.node, job, M);
            continue;
        }
 
        // ---- a service completion ------------------------------------------
        const std::size_t i = ev.station;
        StationState& s = S[i];
        Job job;
 
        if (s.role == Role::Delay) {
            job = ev.job;
            if (has_removal_signal) {
                // The registry is the Delay's cancellation handle: no entry
                // means a signal already destroyed this job, and serving the
                // departure anyway would put it back into circulation.
                std::size_t at = delay_live[i].size();
                for (std::size_t j = 0; j < delay_live[i].size(); ++j)
                    if (delay_live[i][j].id == job.id) {
                        at = j;
                        break;
                    }
                if (at == delay_live[i].size()) continue;
                delay_live[i].erase(delay_live[i].begin() + static_cast<std::ptrdiff_t>(at));
            }
        } else if (s.ps) {
            ps_advance(i);
            std::size_t idx = s.ps_jobs.size();
            for (std::size_t j = 0; j < s.ps_jobs.size(); ++j)
                if (s.ps_jobs[j].tag == ev.tag) {
                    idx = j;
                    break;
                }
            // A superseded departure: the job it named has since been
            // rescheduled by a population change. Dropping it here is the
            // engine's substitute for cancelling the event.
            if (idx == s.ps_jobs.size()) continue;
            const PsJob& pj = s.ps_jobs[idx];
            job.cls = pj.cls;
            job.t_arr = pj.t_arr;
            job.t_sys = pj.t_sys;
            job.parent = pj.parent;
            s.ps_jobs.erase(s.ps_jobs.begin() + static_cast<std::ptrdiff_t>(idx));
        } else if (s.pas) {
            // A stale aggregate clock: the list changed after it was armed, so
            // its total rate no longer describes the station.
            if (s.pas_tag != ev.tag || s.pas_list.empty()) continue;
            job = pas_complete(i);
            acc.update_busy(i, job.cls, now);
            if (acc.busy[i][job.cls] > 0.0) acc.busy[i][job.cls] -= 1.0;
            pas_reschedule(i);
        } else {
            // A preemption leaves its victim's departure in the list with no
            // way to cancel it; the tag mismatch is what neutralises it.
            if (!s.server_busy[ev.slot] || s.server_tag[ev.slot] != ev.tag) continue;
            sd_advance(i);
            job = s.server[ev.slot];
            s.server_busy[ev.slot] = false;
            s.server_tag[ev.slot] = 0;
            acc.update_busy(i, job.cls, now);
            acc.busy[i][job.cls] -= 1.0;
        }
 
        sd_advance(i);
        acc.update_qlen(i, job.cls, now);
        acc.qlen[i][job.cls] -= 1.0;
        bp.track(i, job.cls, -1, now);
        acc.completed[i][job.cls] += 1.0;
        acc.resp_sum[i][job.cls] += now - job.t_arr;
        if (want_respt) resp_samples[i][job.cls].push_back(now - job.t_arr);
        acc.resp_cnt[i][job.cls] += 1.0;
        ++total_completions;
 
        if (s.ps) {
            if (s.lps_limit > 0 && !s.buffer.empty() && s.ps_jobs.size() < s.lps_limit) {
                Job nextjob = buffer_pop(i);
                PsJob pj;
                pj.cls = nextjob.cls;
                pj.priority = nextjob.priority;
                pj.t_arr = nextjob.t_arr;
                pj.t_sys = nextjob.t_sys;
                pj.total = nextjob.service;
                pj.remaining = nextjob.service;
                pj.parent = nextjob.parent;
                s.ps_jobs.push_back(pj);
            }
            ps_reschedule(i);
        }
 
        const RouteEntry re = draw_node_route(sn.station_to_node[i], job.cls);
        const bool to_sink = re.sink;
        const std::size_t dst = re.station, dcls = re.cls;
 
        /**
         * IMMEDIATE FEEDBACK: the completing job goes straight back into THE
         * SAME SERVER under the destination class, holding the slot instead of
         * re-queueing behind whoever is waiting.
         *
         * That is the whole content of the feature and it is not the same model
         * as a self-loop: a self-loop puts the job at the tail and lets the head
         * of the queue in, while feedback lets one job hold the server for as
         * many services as its routing keeps returning it. It also never leaves
         * the station, so it exits no region, releases no blocked upstream
         * server and is not a system completion.
         */
        if (has_immfeed && !to_sink && dst == i && s.role == Role::Queue && !s.ps && !s.pas &&
            ev.slot < s.nservers && dcls < sn.immfeed[i].size() && sn.immfeed[i][dcls]) {
            Job fed;
            fed.cls = dcls;
            fed.t_sys = job.t_sys;
            fed.parent = job.parent;
            fed.t_arr = now;
            fed.priority = classprio[dcls];
            fed.service = slot_snap(s.svc[dcls].next_at(g_svc[i][dcls], now), "service time");
            fed.remaining = fed.service;
            fed.elapsed = 0.0;
            fed.rank = uniform01(g_routing);
            fed.deadline = now + classdeadline[dcls];
            fed.id = ++job_id;
            sd_advance(i);
            acc.update_qlen(i, dcls, now);
            acc.qlen[i][dcls] += 1.0;
            bp.track(i, dcls, +1, now);
            start_service(i, ev.slot, fed);
            sd_reschedule(i);
            if (total_completions >= max_events) break;
            continue;
        }
 
        /**
         * A SYNCHRONOUS CALL PARKS THE SERVER. The completing job leaves for the
         * callee and its slot stays held until the matching reply comes back:
         * no waiter is promoted into it, and the time it spends held counts as
         * BUSY, because the server is unavailable and reporting it idle would
         * credit the station with capacity it does not have. That is the same
         * rule the BAS branch below applies for the same reason.
         *
         * A hop that ALREADY switched into the reply class is not a call: it is
         * the answer, and parking on it would wait for a reply to a reply.
         * That switch may not show on `dcls` yet -- see `resolve_final_cls` --
         * when the routing matrix moves this edge through a synthesized
         * ClassSwitch node, so the check resolves through it rather than
         * reading `dcls` directly.
         */
        if (has_sync_call && sync_reply[job.cls] < K &&
            !is_reply_signal[resolve_final_cls(re.node, dcls)] &&
            s.role == Role::Queue && !s.ps && !s.pas && ev.slot < s.nservers) {
            const std::uint64_t call = ++next_call;
            PendingCall pc;
            pc.station = i;
            pc.slot = ev.slot;
            pc.cls = job.cls;
            pc.since = now;
            pending_reply[call] = pc;
            s.server_held[ev.slot] = true;
            s.held_cls[ev.slot] = job.cls;
            acc.update_busy(i, job.cls, now);
            acc.busy[i][job.cls] += 1.0;
            Job moved;
            moved.cls = re.cls;
            moved.t_sys = job.t_sys;
            moved.parent = job.parent;
            moved.call = call;
            deliver(re.node, moved, i);
            if (total_completions >= max_events) break;
            continue;
        }
 
        // BLOCKING IS DECIDED BEFORE THE SLOT IS REUSED. If the destination is
        // full and declares BAS or BBS, this server keeps the completed job and
        // stops serving; taking the next waiting job first would let the station
        // run at full rate while its output is stopped, which is precisely the
        // behaviour blocking removes.
        if (!to_sink && dst < M && s.role == Role::Queue && !s.ps && ev.slot < s.nservers &&
            !dest_has_room(dst, dcls)) {
            const lang::DropStrategy rule = dest_policy(dst, dcls);
            if (rule == lang::DropStrategy::BAS || rule == lang::DropStrategy::BBS) {
                Job held = job;
                held.cls = dcls;
                s.server_blocked[ev.slot] = true;
                s.blocked_job[ev.slot] = held;
                s.blocked_dest[ev.slot] = dst;
                s.blocked_dest_cls[ev.slot] = dcls;
                // Charged to the DESTINATION it is queueing for, not to the
                // station whose server it occupies.
                // The blocked job IS part of the destination's queue length --
                // the reference's effectiveQueueLength adds basBlockedAtDest and
                // bbsBlockedAtDest to it -- so it is carried in `qlen` and the
                // counter beside it is bookkeeping, not a second population.
                acc.update_qlen(dst, dcls, now);
                acc.qlen[dst][dcls] += 1.0;
                S[dst].blocked_at[dcls] += 1.0;
                // The upstream server is BLOCKED, not idle. Utilization counts
                // blocking time as busy time (the reference divides
                // busy + blocking by the horizon), because the server is
                // unavailable either way and reporting it idle would credit the
                // station with capacity it does not have.
                acc.update_busy(i, job.cls, now);
                acc.busy[i][job.cls] += 1.0;
                blocked_count[i][job.cls] += 1.0;
                continue;
            }
        }
 
        // ONLY NOW may the freed slot take the next waiting job: a slot that
        // blocked above must not also start serving, or the station holds two
        // jobs in one server and the blocked one is duplicated on release.
        if (s.polling) {
            // The controller decides what happens next, not the buffer order:
            // it may keep serving this buffer, or walk on under its budget.
            poll_serve(i);
        } else if (s.role == Role::Queue && !s.ps && !s.pas && !s.server_blocked.empty() &&
                   ev.slot < s.nservers && !s.server_blocked[ev.slot] &&
                   !s.server_held[ev.slot] && buffer_has_for_slot(i, ev.slot)) {
            Job nextjob = buffer_pop_for_slot(i, ev.slot);
            start_service(i, ev.slot, nextjob);
            sd_reschedule(i);
        }
        // THE FORK IDENTITY IS HANDED OVER, NOT COPIED. A continuation that
        // stands in for the trigger at a Join must take its sibling slot; if
        // the trigger kept it too, one fork would present two siblings and the
        // Join would release early and discard the second.
        std::uint64_t spawn_parent = 0;
        if (has_spawn && spawn_of[job.cls] < K && spawn_joins_at[spawn_of[job.cls]]) {
            spawn_parent = job.parent;
            job.parent = 0;
        }
        {
            Job moved;
            moved.cls = re.cls;
            moved.t_sys = job.t_sys;
            moved.parent = job.parent;
            // The identity travels WITH the job, through the callee and back on
            // the reply: a call answered two stations downstream is still this
            // caller's call, and dropping the tag here would strand its server.
            moved.call = job.call;
            deliver(re.node, moved, i);
        }
 
        // A job that left the region frees a slot for whatever is parked on it.
        bool region_to_release = false;
        std::size_t released_region = 0;
        if (region_of[i] >= 0) {
            const std::size_t rg = static_cast<std::size_t>(region_of[i]);
            // The hop is INSIDE the region when the destination is a member:
            // the job never crosses the boundary and the region's occupancy is
            // unchanged, so releasing a waiter there would over-admit.
            // A Fork or Join destination is not a station and has no region,
            // so it is a crossing; `dst` is M there and must not be indexed.
            const bool stays_inside = !to_sink && dst < M && region_of[dst] == region_of[i];
            if (!stays_inside) {
                regions[rg].update(now);
                regions[rg].leave(job.cls);
                regions[rg].completed[job.cls] += 1.0;
                region_to_release = true;
                released_region = rg;
            }
        }
        // THE PHASE-2 CONTINUATION IS INJECTED HERE, between the completing
        // job's release of its region slot and the release of whatever is
        // parked on that region: the continuation takes over the slot the
        // trigger just freed, ahead of the blocked queue, which is the order
        // `maybeSpawnOnCompletion` is called in.
        if (has_spawn && spawn_of[job.cls] < K) {
            const std::size_t scls = spawn_of[job.cls];
            Job spawned;
            spawned.cls = scls;
            spawned.t_sys = now;
            // A continuation aimed at a Join stands in for the trigger at the
            // fork it belongs to (phase 2 at an AND-join branch tail), so it
            // inherits the fork identity rather than arriving as an orphan the
            // Join would discard.
            spawned.parent = spawn_parent;
            admit(i, spawned, M);
        }
        if (region_to_release) release_region(released_region);
        // The station has emptied: arm the delay-off timer, after which the
        // server shuts down and the next arrival pays a setup.
        if (s.has_setup && s.setup_on) {
            double held = 0.0;
            for (std::size_t r = 0; r < K; ++r) held += acc.qlen[i][r];
            if (!(held > 0.0)) {
                const double idle =
                    s.delayoff_time.disabled() ? 0.0 : s.delayoff_time.next(g_aux[i]);
                s.delayoff_at = now + idle;
                Event e;
                e.t = s.delayoff_at;
                e.kind = EV_SETUP;  // reuse the tag; the handler below powers down
                e.station = i;
                e.cls = job.cls;
                e.slot = 1;  // marks a power-down rather than a setup completion
                push(e);
            }
        }
        release_blocked(i);
 
        if (total_completions >= max_events) break;
        if (cnvg.enabled() && (total_completions - last_cnvg_events) >= cnvg.interval()) {
            for (std::size_t a = 0; a < M; ++a) {
                if (S[a].ps) ps_advance(a);
                for (std::size_t r = 0; r < K; ++r) {
                    acc.update_qlen(a, r, now);
                    acc.update_busy(a, r, now);
                }
            }
            cnvg.finalize_batch(acc, servers_of, now);
            last_cnvg_events = total_completions;
            if (cnvg.converged(off_of)) {
                converged = true;
                break;
            }
        }
        if ((total_completions - last_mser_events) >= mser_interval) {
            for (std::size_t a = 0; a < M; ++a) {
                if (S[a].ps) ps_advance(a);
                for (std::size_t r = 0; r < K; ++r) {
                    acc.update_qlen(a, r, now);
                    acc.update_busy(a, r, now);
                }
            }
            obs.collect(acc, now);
            last_mser_events = total_completions;
        }
    }
 
    // Close every open integral at the final instant.
    for (std::size_t i = 0; i < M; ++i) {
        if (S[i].ps) ps_advance(i);
        for (std::size_t r = 0; r < K; ++r) {
            acc.update_qlen(i, r, now);
            acc.update_busy(i, r, now);
        }
    }
 
    bp.commit();
    const Truncation tr = obs.truncate();
    const double sim_time = now - tr.warmup_end;
    const double elapsed = tr.applied ? (now - obs.time[tr.index]) : sim_time;
 
    // ---- result -------------------------------------------------------------
    LdesResult res;
    res.nstations = M;
    res.nclasses = K;
    res.nchains = sn.nchains;
    for (std::size_t i = 0; i < M; ++i) res.station_names.push_back(sn.stations[i].name);
    for (std::size_t r = 0; r < K; ++r) res.class_names.push_back(sn.classes[r].name);
    res.QN = Matrix<double>(M, K, 0.0);
    res.UN = Matrix<double>(M, K, 0.0);
    res.RN = Matrix<double>(M, K, 0.0);
    res.TN = Matrix<double>(M, K, 0.0);
    res.CN = Matrix<double>(1, K, 0.0);
    res.XN = Matrix<double>(1, K, 0.0);
    res.AN = Matrix<double>(M, K, 0.0);
    res.WN = Matrix<double>(M, K, 0.0);
 
    for (std::size_t i = 0; i < M; ++i) {
        for (std::size_t r = 0; r < K; ++r) {
            if (S[i].role == Role::Source) {
                res.TN(i, r) = lambda[r];
                continue;
            }
            if (S[i].off[r]) continue;
            if (elapsed > 0.0) {
                const double q0 = tr.applied ? obs.qt[i][r][tr.index] : 0.0;
                const double b0 = tr.applied ? obs.bt[i][r][tr.index] : 0.0;
                const double c0 = tr.applied ? obs.cmp[i][r][tr.index] : 0.0;
                res.QN(i, r) = (acc.tot_qlen[i][r] - q0) / elapsed;
                res.TN(i, r) = (acc.completed[i][r] - c0) / elapsed;
                if (S[i].role == Role::Delay) {
                    // An infinite server has no capacity to be busy against: the
                    // reference reports the traffic intensity T/mu, which is what
                    // SolverNC and SolverMVA return for a Delay.
                    res.UN(i, r) = res.TN(i, r) * S[i].class_mean[r];
                } else if (S[i].has_cd) {
                    // A CLASS- OR JOINT-DEPENDENT STATION IS READ THROUGH THE
                    // LAW, not the busy integral, which is what the Java engine
                    // does (`getUtilization` returns T/(mu*peakScaling) here).
                    // beta_r(n) scales the RATE, so a busy server delivers
                    // beta times the work; the head-count integral would report
                    // the fraction of TIME occupied and divide it by a peak that
                    // counts work, mixing the two conventions in one number. The
                    // load-dependent case keeps the integral because
                    // `Accum::busy_scale` already weights it by alpha(n).
                    res.UN(i, r) = (S[i].util_peak > 0.0)
                                       ? res.TN(i, r) * S[i].class_mean[r] / S[i].util_peak
                                       : 0.0;
                } else {
                    res.UN(i, r) = (acc.tot_busy[i][r] - b0) / (elapsed * S[i].util_peak);
                }
            }
            if (acc.resp_cnt[i][r] > 0.0) res.RN(i, r) = acc.resp_sum[i][r] / acc.resp_cnt[i][r];
            // AN is the OFFERED rate, arrivedCustomers / simTime, as
            // `getArrivalRate` measures it; WN is the residence time, which for
            // this engine's single-visit accounting is the response time. Both
            // are part of the result contract every client reads, and leaving
            // them at zero reads as "no arrivals" rather than "not measured".
            if (elapsed > 0.0) res.AN(i, r) = acc.arrived[i][r] / elapsed;
            res.WN(i, r) = res.RN(i, r);
        }
    }
    // Sibling drops at a Join, on the Join rows only. Emitted whenever the model
    // HAS a fork-join, drops or none: a MEASURED zero is what a standard join
    // loses, and it is a better answer than the client's derived AN - K*TN,
    // which carries the finite-sample gap between the two rates. An absent
    // matrix means "this engine does not count", not "nothing was lost".
    if (!sn.fj.empty()) {
        res.DropRateJoin = Matrix<double>(M, K, 0.0);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                const double d0 = tr.applied ? obs.drp[i][r][tr.index] : 0.0;
                if (elapsed > 0.0)
                    res.DropRateJoin(i, r) = (acc.join_dropped[i][r] - d0) / elapsed;
            }
    }
    for (std::size_t r = 0; r < K; ++r) {
        if (sim_time > 0.0) res.XN(0, r) = sys_completed[r] / sim_time;
        if (sys_resp_cnt[r] > 0.0) res.CN(0, r) = sys_resp_sum[r] / sys_resp_cnt[r];
    }
 
    // ---- busy periods --------------------------------------------------------
    for (std::size_t ti = 0; ti < bp.targets().size(); ++ti) {
        LdesResult::BusyPeriodTarget out;
        out.name = bp.targets()[ti].name;
        out.stations = bp.targets()[ti].stations;
        out.job_class = bp.targets()[ti].job_class;
        for (std::size_t oi = 0; oi < static_cast<std::size_t>(bp.orders()); ++oi) {
            out.mean.push_back(bp.mean(ti, oi));
            out.count.push_back(bp.count(ti, oi));
        }
        res.busy_periods.push_back(out);
    }
 
    // ---- transient trajectory and histogram -----------------------------------
    if (transient_run) {
        now = std::min(now, horizon);
        hist_accumulate();
        res.t = tran_times;
        res.QNt.assign(M, std::vector<Matrix<double>>(K));
        res.UNt.assign(M, std::vector<Matrix<double>>(K));
        res.TNt.assign(M, std::vector<Matrix<double>>(K));
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                const std::size_t n = tran_times.size();
                // Columns are [value, time], which is the layout the reference
                // emits and the clients read.
                Matrix<double> q(n, 2, 0.0), u(n, 2, 0.0), t(n, 2, 0.0);
                for (std::size_t k = 0; k < n; ++k) {
                    q(k, 0) = tran_q[i][r][k];
                    q(k, 1) = tran_times[k];
                    u(k, 0) = tran_u[i][r][k];
                    u(k, 1) = tran_times[k];
                    t(k, 0) = tran_t[i][r][k];
                    t(k, 1) = tran_times[k];
                }
                res.QNt[i][r] = q;
                res.UNt[i][r] = u;
                res.TNt[i][r] = t;
            }
        res.stopping_reason = "max_time";
    }
 
    if (want_respt) {
        res.respTimeSamples.assign(M, std::vector<std::vector<double>>(K));
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) res.respTimeSamples[i][r] = resp_samples[i][r];
    }
 
    // ---- the joint-state histogram and the state path -------------------------
    // OUTSIDE the transient branch: a steady-state run that asked for either
    // fills it too, which is what `--export-histogram` and `--trajectory` mean.
    // The last interval is closed here, or the state the run ended in would
    // carry no residence time at all.
    if (want_hist && !transient_run) hist_accumulate();
    if (!histogram.empty()) {
        res.histogram_space = Matrix<double>(histogram.size(), M * K, 0.0);
        res.histogram_time = Matrix<double>(histogram.size(), 1, 0.0);
        std::size_t row = 0;
        for (const auto& kv : histogram) {
            for (std::size_t c = 0; c < kv.first.size(); ++c)
                res.histogram_space(row, c) = kv.first[c];
            res.histogram_time(row, 0) = kv.second;
            ++row;
        }
    }
    if (!trajectory.empty()) {
        res.traj_space = Matrix<double>(trajectory.size(), M * K, 0.0);
        res.traj_time = Matrix<double>(trajectory.size(), 1, 0.0);
        for (std::size_t k = 0; k < trajectory.size(); ++k) {
            for (std::size_t c = 0; c < trajectory[k].second.size(); ++c)
                res.traj_space(k, c) = trajectory[k].second[c];
            res.traj_time(k, 0) = trajectory[k].first;
        }
    }
 
    // ---- cache metrics --------------------------------------------------------
    for (const auto& kv : caches) {
        const CacheState& cs = kv.second;
        LdesCacheMetrics cm;
        cm.hit = Matrix<double>(1, K, 0.0);
        cm.miss = Matrix<double>(1, K, 0.0);
        // A DELAYED HIT IS ITS OWN OUTCOME, so it enters the denominator beside
        // the other two: with a retrieval system the three fractions partition
        // the reads, and dividing by hits + misses alone would report fractions
        // that sum above one.
        if (cs.has_retrieval) cm.delayed = Matrix<double>(1, K, 0.0);
        for (std::size_t r = 0; r < K; ++r) {
            const double dl = cs.has_retrieval ? cs.delayed[r] : 0.0;
            const double tot = cs.hits[r] + cs.misses[r] + dl;
            if (tot > 0.0) {
                cm.hit(0, r) = cs.hits[r] / tot;
                cm.miss(0, r) = cs.misses[r] / tot;
                if (cs.has_retrieval) cm.delayed(0, r) = dl / tot;
            }
        }
        // The measured retrieval latency: the mean time a request WAITED on the
        // retrieval system, averaged over BOTH populations that wait -- the one
        // that triggered the fetch and every request parked behind it. A parked
        // request waits the RESIDUAL of the fetch period, not the whole of it,
        // so averaging the fetch sojourns alone answers a different question:
        // it agrees only where the fetch period is memoryless, and reads 6.5%
        // high on `retrieval_chain`, whose fetch is a two-stage hypoexponential.
        //
        // BOTH POPULATIONS IS WHAT MAKES THE PAIR (ArvR, ResidT) MEAN ANYTHING.
        // `getAvgCacheTable` reports ArvR = lambda*(miss + delayed), the rate
        // INTO the retrieval system, so the residence time beside it must be
        // Little-consistent with that rate: E[N] = miss*W_fetch + delayed*W_parked
        // over the same arrivals. That is also what the analytic reference
        // computes -- `retrieval_fpi_latency`'s Z = sum(phi + d) / sum(lambda *
        // (phi + pi0)) is Little's law over the same two populations -- and what
        // `Solver_ssj.exportCacheResults` measures. Reported only where measured.
        const double released = std::accumulate(cs.delayed.begin(), cs.delayed.end(), 0.0);
        if (cs.has_retrieval && cs.completed_fetches + released > 0.0) {
            cm.latency = Matrix<double>(1, K, 0.0);
            const double mean_wait = (cs.total_fetch_time + cs.delayed_wait)
                                     / (cs.completed_fetches + released);
            for (std::size_t r = 0; r < K; ++r)
                cm.latency(0, r) = (cs.hits[r] + cs.misses[r] + cs.delayed[r] > 0.0)
                                       ? mean_wait
                                       : std::numeric_limits<double>::quiet_NaN();
        }
        res.cache_metrics[sn.nodes[cs.node - 1].name] = cm;
    }
 
    // ---- finite capacity regions ---------------------------------------------
    if (!regions.empty()) {
        const std::size_t NR = regions.size();
        res.nregions = NR;
        res.QNfcr = Matrix<double>(NR, K, 0.0);
        res.TNfcr = Matrix<double>(NR, K, 0.0);
        res.WeightNfcr = Matrix<double>(NR, K, 0.0);
        res.MemOccNfcr = Matrix<double>(NR, K, 0.0);
        res.DropRateNfcr = Matrix<double>(NR, K, 0.0);
        for (std::size_t g = 0; g < NR; ++g) {
            regions[g].update(now);
            for (std::size_t r = 0; r < K; ++r) {
                if (sim_time > 0.0) {
                    // JMT's convention: the region's queue length counts only
                    // the jobs INSIDE it. A WAITQ job parked outside has not
                    // entered and is not part of the region's occupancy.
                    res.QNfcr(g, r) = regions[g].tot_jobs[r] / sim_time;
                    res.WeightNfcr(g, r) = regions[g].tot_weight[r] / sim_time;
                    res.MemOccNfcr(g, r) = regions[g].tot_mem[r] / sim_time;
                    res.TNfcr(g, r) = regions[g].completed[r] / sim_time;
                    res.DropRateNfcr(g, r) = regions[g].dropped[r] / sim_time;
                }
            }
        }
    }
 
    // ---- confidence intervals ------------------------------------------------
    // The event count is set BEFORE the intervals, not after: `batch_means_ci`
    // plans a run length from it when the caller asked for one, and reading a
    // zero there would silently fall back to the BUDGET, which an early
    // convergence stop never spent.
    res.total_simulated_events = static_cast<long long>(total_completions);
    batch_means_ci(obs, tr, o, res);
 
    // The impatience block is reported only when the model declares one: an
    // all-zero matrix would read as "measured, and none happened".
    double any_balk = 0.0, any_renege = 0.0;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) {
            any_balk += balked[i][r];
            any_renege += reneged[i][r];
        }
    double any_orbit = 0.0;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) any_orbit += S[i].retried[r];
    if (any_orbit > 0.0) {
        res.retriedCustomers = Matrix<double>(M, K, 0.0);
        res.retrialDropped = Matrix<double>(M, K, 0.0);
        res.avgOrbitSize = Matrix<double>(M, K, 0.0);
        for (std::size_t i = 0; i < M; ++i) {
            const double dt = now - S[i].orbit_last;
            for (std::size_t r = 0; r < K; ++r) {
                double tot = S[i].tot_orbit[r];
                if (dt > 0.0) tot += S[i].orbit_size[r] * dt;
                res.retriedCustomers(i, r) = S[i].retried[r];
                res.retrialDropped(i, r) = S[i].retrial_lost[r];
                if (sim_time > 0.0) res.avgOrbitSize(i, r) = tot / sim_time;
            }
        }
    }
 
    if (any_balk > 0.0 || any_renege > 0.0) {
        res.balkedCustomers = Matrix<double>(M, K, 0.0);
        res.renegedCustomers = Matrix<double>(M, K, 0.0);
        res.renegingRate = Matrix<double>(M, K, 0.0);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                res.balkedCustomers(i, r) = balked[i][r];
                res.renegedCustomers(i, r) = reneged[i][r];
                if (sim_time > 0.0) res.renegingRate(i, r) = reneged[i][r] / sim_time;
            }
    }
 
    res.converged = converged;
    res.convergence_batches = cnvg.batches();
    if (!transient_run) res.stopping_reason = converged ? "convergence" : "max_events";
    res.engine = "native";
    res.method = o.method;
    return res;
}
 
/**
 * Simulate `sn`, over independent REPLICATIONS when `options.replications > 1`.
 *
 * Each replication is one path with its OWN seed, `seed + r`, exactly as
 * `Solver_ldes_analyzer_parallel` derives them. That is what makes the paths
 * independent while keeping the whole run reproducible; sharing a seed would
 * make every replication identical and collapse the interval to zero, and
 * drawing them randomly would make the run irreproducible.
 *
 * THE INTERVAL COMES FROM THE SPREAD BETWEEN REPLICATIONS, not from within a
 * path, and it is the honest one: independent paths need no assumption about
 * the correlation structure that batch means has to model. The cost is that
 * `replications - 1` degrees of freedom is few, so the t critical value is
 * large and the interval wide -- which is the estimate, not a defect of it.
 *
 * A single replication yields NO interval rather than a zero-width one: with
 * one path there is no spread to measure, and reporting zero would read as
 * perfect precision.
 */
template <class T>
LdesResult ldes_engine_solve(const qn::NetworkStruct<T>& sn, const LdesOptions& o) {
    using namespace engine;
    const int R = o.replications;
    if (R <= 1) return ldes_engine_solve_one(sn, o);
 
    std::vector<LdesResult> reps;
    reps.reserve(static_cast<std::size_t>(R));
    for (int r = 0; r < R; ++r) {
        LdesOptions ro = o;
        ro.replications = 0;
        // A replication measures its OWN interval nowhere: the spread across
        // paths is the estimate, so the within-path machinery is switched off.
        ro.cimethod = "none";
        ro.seed = (o.seed >= 0) ? (o.seed + r) : o.seed;
        reps.push_back(ldes_engine_solve_one(sn, ro));
    }
 
    LdesResult res = reps[0];
    res.QNCI = Matrix<double>();
    res.UNCI = Matrix<double>();
    res.RNCI = Matrix<double>();
    res.TNCI = Matrix<double>();
 
    const std::size_t M = res.nstations, K = res.nclasses;
    Matrix<double>* means[4] = {&res.QN, &res.UN, &res.RN, &res.TN};
    Matrix<double>* cis[4] = {&res.QNCI, &res.UNCI, &res.RNCI, &res.TNCI};
    const Matrix<double>* per_rep[4];
 
    for (int q = 0; q < 4; ++q) {
        *means[q] = Matrix<double>(M, K, 0.0);
        *cis[q] = Matrix<double>(M, K, 0.0);
    }
    const double tcrit = t_critical(o.confint > 0.0 ? o.confint : 0.95, R - 1);
 
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t k = 0; k < K; ++k)
            for (int q = 0; q < 4; ++q) {
                double sum = 0.0;
                for (int r = 0; r < R; ++r) {
                    per_rep[q] = (q == 0)   ? &reps[static_cast<std::size_t>(r)].QN
                                 : (q == 1) ? &reps[static_cast<std::size_t>(r)].UN
                                 : (q == 2) ? &reps[static_cast<std::size_t>(r)].RN
                                            : &reps[static_cast<std::size_t>(r)].TN;
                    sum += (*per_rep[q])(i, k);
                }
                const double mean = sum / R;
                (*means[q])(i, k) = mean;
                double ss = 0.0;
                for (int r = 0; r < R; ++r) {
                    per_rep[q] = (q == 0)   ? &reps[static_cast<std::size_t>(r)].QN
                                 : (q == 1) ? &reps[static_cast<std::size_t>(r)].UN
                                 : (q == 2) ? &reps[static_cast<std::size_t>(r)].RN
                                            : &reps[static_cast<std::size_t>(r)].TN;
                    const double d = (*per_rep[q])(i, k) - mean;
                    ss += d * d;
                }
                (*cis[q])(i, k) = tcrit * std::sqrt((ss / (R - 1)) / R);
            }
 
    // The system-level rows average the same way; they carry no interval,
    // matching what the single-path engine reports for them.
    res.XN = Matrix<double>(1, K, 0.0);
    res.CN = Matrix<double>(1, K, 0.0);
    for (std::size_t k = 0; k < K; ++k) {
        double sx = 0.0, sc = 0.0;
        for (int r = 0; r < R; ++r) {
            sx += reps[static_cast<std::size_t>(r)].XN(0, k);
            sc += reps[static_cast<std::size_t>(r)].CN(0, k);
        }
        res.XN(0, k) = sx / R;
        res.CN(0, k) = sc / R;
    }
 
    long long total = 0;
    for (int r = 0; r < R; ++r) total += reps[static_cast<std::size_t>(r)].total_simulated_events;
    res.total_simulated_events = total;
    return res;
}
 
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_ENGINE_H