ldes_ln_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_LN_ENGINE_H
#define LINE_SOLVERS_LDES_LDES_LN_ENGINE_H
 
/**
 * @file
 * @ingroup line_solvers
 * The NATIVE LDES engine for LAYERED (LQN) models, the C++ twin of
 * `jline/solvers/ldes/handlers/Solver_ssj_ln.java`.
 *
 * IT IS A DIFFERENT SIMULATOR FROM THE FLAT ONE, not a wrapper around it. A
 * layered model has no jobs circulating over a routing matrix: it has
 * REFERENCE TASKS that think and then invoke an entry, ACTIVITIES that consume
 * host demand, and CALLS that suspend the caller until the callee replies. The
 * decomposition into layers that `SolverLN` performs is an APPROXIMATION; this
 * engine simulates the layered semantics directly and is therefore the
 * reference the layered solvers are checked against.
 *
 * THE TWO RESOURCES ARE HELD AT ONCE, and that is the whole content of a
 * layered model. An activity needs a THREAD of its task and a SERVER of its
 * host, and it holds the thread across a synchronous call while the callee runs
 * on a different host entirely. Releasing the thread during the call would turn
 * every synchronous call into an asynchronous one and remove the layered
 * contention the model exists to represent -- the answer stays plausible and
 * every utilization falls.
 *
 * WHAT IT COVERS: reference tasks with think time and multiplicity; task
 * multiplicity as a thread semaphore, with setup/delay-off power cycling;
 * host demand on c-server processors under FCFS, LCFS, SIRO, HOL, FCFSPRIO,
 * LCFSPRIO, FCFSPRPRIO, FCFSPIPRIO, LCFSPRPRIO, LCFSPIPRIO, PS, PSPRIO and INF,
 * and task queueing under FCFS, LCFS, SIRO, HOL, FCFSPRIO, LCFSPRIO and INF,
 * each served as Solver_ssj_ln serves it. The priority disciplines order by
 * `lsn.prio` on the lqns scale, a LARGER number first and the default 0 lowest;
 * HOL and FCFSPRIO are FIFO within a level, the LCFS variants LIFO; the PR/PI
 * variants preempt, resuming the residual work (PR) or drawing a fresh demand
 * (PI); PS is an egalitarian share min(1, c/n) of the c servers and PSPRIO
 * gives the servers to the levels from the top down; INF never queues, whatever
 * the multiplicity. Also synchronous calls with a mean call multiplicity,
 * asynchronous calls and open arrivals (each a request nobody waits on),
 * activity think times (after the demand, before the calls, holding the thread
 * but not the processor), REPLICATION of a processor or a task, one queue per
 * copy, and, as `Solver_ssj_ln` runs them:
 *
 *  - ACTIVITY GRAPHS walked as graphs, not as a list. An entry starts at its
 *    bound activity and follows `graph`: one successor continues, an OR fork
 *    or a loop exit draws one branch, an AND fork splits the invocation into
 *    concurrent LINES sharing its thread, and an AND join fires on the root
 *    line once every PRE_AND input has arrived. A phase-1 activity followed by
 *    phase 2 replies early, and the thread is held through phase 2.
 *  - CACHE TASKS: one content per task replica, an ItemEntry read drawing an
 *    item from its pmf and continuing on the hit or the miss branch, the
 *    replacement rules RR/FIFO/SFIFO/LRU/HLRU/CLIMB/QLRU over the level lists,
 *    and delayed-hit retrieval (a read of an in-flight item is parked until
 *    the fetcher replies, then continues on the hit branch).
 *  - ADMISSION CONSTRAINTS A n <= b: a host row blocks a host request outside
 *    the processor (thread kept, residence running), a task row blocks a
 *    request before the task (counted nowhere), FIFO release on a departure,
 *    and a run that stops with a blocked request is reported as a deadlock.
 *  - HETEROGENEOUS SERVER POOLS on a processor, served as OI rate sharing: each
 *    pool puts its whole capacity count*rate on the first held job it serves.
 *
 * An activity runs its host demand FIRST, then its think time, then its calls,
 * as LQN2QN builds it and Solver_ssj_ln runs it. An entry's response time is
 * its service time, from taking a thread of its task to its reply: the wait for
 * the thread belongs to the caller's call. Its occupancy runs from the same
 * instant to the end of its last phase.
 *
 * WHAT IT REFUSES by name: every discipline the Java engine refuses (PS or
 * preemption on a task, which holds threads and does not divide or preempt
 * them, and any weighted discipline anywhere), a setup on an infinite-thread
 * task, server pools on a task, a partial AND-join quorum (the Java engine waits
 * for every input), and a cache read whose ItemEntry has no popularity or no
 * hit/miss pair.
 *
 * WHAT IT STILL DOES NOT SIMULATE, silently, as before: forwarding.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <functional>
#include <limits>
#include <list>
#include <map>
#include <queue>
#include <set>
#include <string>
#include <vector>
 
#include "line/lang/lqn/lqn_struct.h"
#include "line/solvers/ldes/ldes_sampler.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 {
 
/** One scheduled event of the layered engine. */
struct LnEvent {
    double t = 0.0;
    /**
     * 0 = host completion (`who` = the LINE), 1 = think completion (`who` = the
     * reference customer), 2 = a task thread finished its cold start, 3 = an
     * idle task thread's delay-off countdown expired, 4 = the next completion
     * on a pooled host replica (`who` = the host slot), 5 = the next completion
     * on a PS host replica (`who` = the host slot), 6 = line `who` finished the
     * think time of its activity, 7 = an open arrival at entry `who` on replica
     * `aux` of its task. A kind-0 event carries the line's host generation, so a
     * preemption makes it stale.
     *
     * For kinds 2 and 3 `who` is the task SLOT and `aux` the thread within it,
     * and `gen` is the generation the timer was armed with: a countdown caught
     * by an arriving request is cancelled by bumping the generation, so the
     * event still fires but finds itself stale and does nothing. Kinds 4 and 5 use the
     * same device per host slot, since every change of the job set moves the
     * next completion. The event queue has no removal.
     */
    int kind = 0;
    std::size_t who = 0;
    std::size_t aux = 0;
    std::uint64_t gen = 0;
    std::uint64_t seq = 0;
};
 
struct LnEventLater {
    bool operator()(const LnEvent& a, const LnEvent& 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;
    }
};
 
/** Index meaning "none" for a line, an invocation or a customer. */
static const std::size_t LN_NONE = static_cast<std::size_t>(-1);
 
/**
 * One LINE of execution: a point in the activity graph of one invocation.
 *
 * A sequential invocation has one line, its root. An AND fork gives it one line
 * per branch, all holding the SAME thread of the invocation, and the join
 * resumes the root; that is `LNRequest.copyWith` of `Solver_ssj_ln`.
 */
struct LnLine {
    std::size_t inv = LN_NONE;  ///< the invocation this line belongs to
    std::size_t act = 0;        ///< current activity, 0 for an entry with none
    /**
     * When the current activity began, -1 before it is reached. Its span runs
     * to its host completion, so it CARRIES THE NESTED CALLS it makes.
     */
    double act_t0 = -1.0;
    std::size_t call_pos = 0;           ///< next call index of the current activity
    /**
     * Repetitions of the current call still to be made; `NOTDRAWN` until drawn.
     * `calls-mean` is drawn once when the call is first reached, so a call of
     * mean 3 blocks three times per execution of the activity.
     */
    std::size_t calls_left = static_cast<std::size_t>(-1);
    double host_t0 = 0.0;   ///< when the line asked its processor, queueing included
    bool branch = false;    ///< a line an AND fork created (never the root)
    int hcol = -1;          ///< host admission column held, -1 = none
    int pcol = -1;          ///< operand column at a pooled host
    double remaining = 0.0; ///< work left at a pooled or PS host
    double siro_key = 0.0;  ///< uniform key in a SIRO host queue; the smallest is served next
    /**
     * How far the current activity is before its calls: 0 nothing started, 1 host demand started
     * or none, 2 think time started or none.
     */
    int dem_done = 0;
    /**
     * Host service in progress at a queueing host. The completion event carries `host_gen`, drawn
     * from one engine-wide counter because lines are recycled, so a preemption cancels it by
     * drawing a new one. `host_rem` is the work left when service (re)starts, negative meaning
     * "draw the demand": a fresh request or a PI victim.
     */
    std::uint64_t host_gen = 0;
    double host_rem = -1.0;
    double host_start = 0.0;
    double host_end = 0.0;
};
 
/** One invocation of an entry: what a call or a reference cycle starts. */
struct LnInv {
    std::size_t entry = 0, task = 0;
    std::size_t trep = 0;             ///< replica of the task serving it
    std::size_t ts = 0;               ///< task slot of that replica
    std::size_t th = static_cast<std::size_t>(-1);  ///< thread held, NOTHR on an INF task
    bool has_thread = false;
    bool started = false;             ///< the entry's service has begun (a thread is held)
    double entry_t0 = 0.0;            ///< when the entry took its thread, the start of its service
    std::size_t caller = LN_NONE;     ///< the line blocked on this call, if synchronous
    std::size_t cust = LN_NONE;       ///< the reference customer, at the top level only
    std::size_t root = LN_NONE;       ///< root line
    bool replied = false;
    bool early_reply = false;         ///< replied before completing (phase 2 follows)
    long pending = -1;                ///< live lines after a fork, -1 = never forked
    std::map<std::size_t, std::set<std::size_t> > joins;  ///< join activity -> inputs arrived
    int tcol = -1;                    ///< task admission column held, -1 = none
    bool fetching = false;            ///< this read is fetching an item (retrieval)
    std::size_t fetch_cache = 0, fetch_item = 0;
    double siro_key = 0.0;            ///< uniform key in a SIRO thread queue
};
 
/** A reference-task customer: it thinks, then starts an invocation. */
struct LnCustomer {
    std::size_t ref_task = 0;
    std::size_t repl = 0;   ///< replica of the reference task it belongs to
    double t_start = 0.0;   ///< when the current cycle began
};
 
/** Live content of one replica of a cache task, `LNCacheState` of the Java engine. */
struct LnCacheState {
    std::size_t nitems = 0;
    std::vector<int> cap;                       ///< per level list
    lang::ReplacementStrategy rs = lang::ReplacementStrategy::LRU;
    std::vector<std::list<std::size_t> > levels;  ///< per level, most recent first
    bool retrieval = false;
    std::vector<char> inflight;                 ///< [item] a fetch is in progress
    std::vector<std::vector<std::pair<std::size_t, std::size_t> > > held;  ///< (line, access)
    long long reads = 0, hits = 0, misses = 0, delayed = 0;
};
 
/** One ItemEntry read: its driver activity, its item law and its two branches. */
struct LnCacheAccess {
    std::vector<std::size_t> caches;  ///< cache state per replica of the cache task
    std::vector<double> cdf;
    std::size_t hit = 0, miss = 0, entry = 0;
};
 
/** Thread index meaning "no thread held", also used for an infinite-thread task. */
static const std::size_t NOTHR = static_cast<std::size_t>(-1);
 
/** `LnJob::calls_left` before the repetition count of a call has been drawn. */
static const std::size_t NOTDRAWN = static_cast<std::size_t>(-1);
 
/** Power state of one task thread, mirroring Solver_ssj_ln's THREAD_* constants. */
enum class ThreadState { ACTIVE, DELAYOFF, OFF, SETUP };
 
/** The layered result: per element, the mean measures. */
struct LnResult {
    std::size_t nidx = 0;
    Matrix<double> QLN, ULN, RLN, TLN;
    /**
     * Residence time per element, the ResidT column: the time an activity holds
     * ITS HOST PROCESSOR per visit of the request stream driving its task, and
     * for a task the sum over its activities.
     *
     * NOT the response time RLN, which also carries the nested synchronous
     * calls the activity makes. NaN at the element kinds that have no residence
     * -- processors and entries -- because a zero residence is a claim (nothing
     * is held there) and this is the absence of one, the same convention
     * SolverLN and LQNS report.
     */
    Matrix<double> WLN;
    /**
     * Every per-request ENTRY response time observed, one vector per entry in
     * LOCAL index space (0..nentries-1). `RLN` at the entry indices is their
     * mean; these are the observations themselves, which is what an empirical
     * CDF has to be built from -- a law fitted to the mean says nothing about
     * the tail, and the tail is the reason to ask a simulator at all.
     */
    std::vector<std::vector<double> > entry_resp_samples;
    double simulated_time = 0.0;
    long long completions = 0;
    /**
     * Cache measures, one row per ItemEntry that was read, pooled over the
     * replicas of its cache task as `Solver_ssj_ln` pools them: element indices
     * of the task and of the entry, the hit, miss and delayed-hit probabilities
     * per read, and the read rate. Hit + miss + delayed = 1.
     */
    std::vector<std::size_t> cache_task, cache_entry;
    std::vector<double> cache_hit, cache_miss, cache_delayed, cache_read_rate;
};
 
/** One entry's empirical response-time CDF: F[j] = P(R <= t[j]). */
struct LnEntryCdf {
    std::vector<double> t;
    std::vector<double> F;
};
 
/**
 * The empirical response time CDF of every ENTRY, the `getCdfRespTLN` of the
 * other codebases: one [F(t), t] table per entry in the LOCAL index space,
 * empty where the run observed nothing.
 *
 * The ecdf of each entry's observations: sort, step by 1/n, and collapse ties
 * keeping the LARGEST value at each distinct time, so an interpolating reader
 * cannot land on a multivalued point.
 */
inline std::vector<LnEntryCdf> ldes_ln_cdf_respt(const LnResult& r, std::size_t nentries) {
    std::vector<LnEntryCdf> out(nentries);
    for (std::size_t e = 0; e < nentries; ++e) {
        if (e >= r.entry_resp_samples.size()) continue;
        std::vector<double> x = r.entry_resp_samples[e];
        if (x.empty()) continue;
        std::sort(x.begin(), x.end());
        const double n = static_cast<double>(x.size());
        for (std::size_t i = 0; i < x.size(); ++i) {
            if (i + 1 < x.size() && x[i + 1] == x[i]) continue;
            out[e].t.push_back(x[i]);
            out[e].F.push_back(static_cast<double>(i + 1) / n);
        }
    }
    return out;
}
 
/** A column of the shared layered average table, as `line-cli` prints it. */
enum class LnColumn { QLen, Util, RespT, ResidT, Tput };
 
/**
 * Does the element at `i` HAVE the quantity in column `c`?
 *
 * THE MASK BELONGS TO THE TABLE, not to whichever engine filled it. A NaN there
 * is not a failed computation: it says the quantity is not defined for that
 * element kind, and SolverLN and LQNS report the same one --
 *
 *   Processor   QLen no,  Util yes, RespT no,  ResidT no,  Tput no
 *   Task        QLen yes, Util yes, RespT no,  ResidT yes, Tput yes
 *   Entry       QLen yes, Util yes, RespT yes, ResidT no,  Tput yes
 *   Activity    all five
 *
 * -- so a reader can diff the arms row by row. This engine MEASURES more than
 * that: a processor's completion rate is sitting in `TLN` and an entry's
 * occupancy in `QLN`, both perfectly well defined as sample-path quantities.
 * Reporting them where every other solver says NaN would make one column mean
 * different things depending on who filled it, which is the divergence the JAR
 * removed from `getLNAvgTable` on 2026-08-21. Nothing is discarded: the caller
 * still has the whole `LnResult`; only the SHARED TABLE is masked.
 *
 * ResidT needs no rule of its own here because the engine already applies it:
 * `WLN` starts at NaN and only activities and their tasks are written, so a
 * task with no activity at all keeps its NaN rather than claiming a residence
 * of zero (which would say something is held there for no time, not that
 * nothing is held).
 *
 * Twin of the mask `jline.solvers.ldes.SolverLDES.getLNAvgTable` applies and of
 * the `defined_*` flags of `ln::LnSolution`; the three are pinned against each
 * other by `cpp/tests/test_lqn_nan_mask_parity.cpp`.
 */
template <class T>
inline bool ln_defined(const lqn::LqnStruct<T>& lsn, const LnResult& r, std::size_t i,
                       LnColumn c) {
    if (i < 1 || i > r.nidx || i >= lsn.type.size()) return false;
    const bool host = lsn.type[i] == line::lang::LqnElement::HOST;
    switch (c) {
        case LnColumn::QLen: return !host;
        case LnColumn::Util: return true;
        case LnColumn::RespT:
            return lsn.type[i] == line::lang::LqnElement::ENTRY ||
                   lsn.type[i] == line::lang::LqnElement::ACTIVITY;
        case LnColumn::ResidT: return !std::isnan(r.WLN(i, 0));
        case LnColumn::Tput: return !host;
    }
    return false;
}
 
}  // namespace engine
 
/**
 * Simulate a layered model in process.
 *
 * The budget is ENTRY COMPLETIONS, mirroring the flat engine's service
 * completions: a layered model has no single notion of "a job leaving", so the
 * count of entry invocations that returned is what a horizon can be set on.
 */
template <class T>
engine::LnResult ldes_ln_engine_solve(const lqn::LqnStruct<T>& lsn, const LdesOptions& o) {
    using namespace engine;
    using lang::SchedStrategy;
 
    const std::size_t nidx = lsn.nidx;
    if (nidx == 0) throw InputError("SolverLDES (native LN engine): empty layered model");
 
    // ---- refuse what the engine does not simulate ---------------------------
    // The disciplines Solver_ssj_ln.assertSchedulingIsSupported admits, refused by name otherwise:
    // a discipline quietly served as another returns plausible numbers for a model nobody described.
    for (std::size_t k = 1; k <= nidx && k < lsn.sched.size(); ++k) {
        const bool is_host = k > lsn.hshift && k <= lsn.hshift + lsn.nhosts;
        const bool is_task = k > lsn.tshift && k <= lsn.tshift + lsn.ntasks;
        if (!is_host && !is_task) continue;
        const SchedStrategy sc = lsn.sched[k];
        const bool ok = sc == SchedStrategy::FCFS || sc == SchedStrategy::INF ||
                        sc == SchedStrategy::LCFS || sc == SchedStrategy::SIRO ||
                        sc == SchedStrategy::HOL || sc == SchedStrategy::REF ||
                        sc == SchedStrategy::LCFSPRIO ||
                        (is_host && (sc == SchedStrategy::PS || sc == SchedStrategy::PSPRIO ||
                                     sc == SchedStrategy::FCFSPRPRIO ||
                                     sc == SchedStrategy::FCFSPIPRIO ||
                                     sc == SchedStrategy::LCFSPRPRIO ||
                                     sc == SchedStrategy::LCFSPIPRIO));
        if (ok) continue;
        std::string scname = lang::sched_to_text(sc);
        for (std::size_t c = 0; c < scname.size(); ++c)
            scname[c] = static_cast<char>(std::toupper(static_cast<unsigned char>(scname[c])));
        throw UnsupportedError(
            "SolverLDES (native LN engine) does not simulate " + scname + " scheduling at " +
            (is_host ? "processor" : "task") + " '" + lsn.names[k] +
            "'. The layered engine implements " +
            (is_host ? "FCFS, LCFS, SIRO, HOL, FCFSPRIO, LCFSPRIO, FCFSPRPRIO, FCFSPIPRIO, "
                       "LCFSPRPRIO, LCFSPIPRIO, PS, PSPRIO and INF"
                     : "FCFS, LCFS, SIRO, HOL, FCFSPRIO, LCFSPRIO and INF (a task holds threads, "
                       "it does not divide or preempt them)") +
            "; a weighted discipline would additionally need the per-task shares the layered "
            "struct does not carry.");
    }
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        // A setup IS simulated, except on an infinite-thread task, which holds
        // no thread to power down. Same rule as the Java engine.
        if (t < lsn.hassetup.size() && lsn.hassetup[t]) {
            const double m = (t < lsn.mult.size()) ? lsn.mult[t] : 1.0;
            if (!std::isfinite(m) || (t < lsn.sched.size() && lsn.sched[t] == lang::SchedStrategy::INF))
                throw UnsupportedError("SolverLDES (native LN engine): task '" + lsn.names[t] +
                                       "' declares a setup time on an infinite-server task, which "
                                       "holds no thread to power down; give it a finite "
                                       "multiplicity");
        }
        // Server pools on a TASK: the Java engine holds pools on processors only
        // and throws here, so this does too rather than pick a reading of its own.
        if (t < lsn.pools.size() && lsn.pools[t].npools() > 0)
            throw UnsupportedError("SolverLDES (native LN engine): task '" + lsn.names[t] +
                                   "' declares heterogeneous server pools; LDES holds server "
                                   "pools on processors only");
    }
 
    const std::uint64_t max_events =
        static_cast<std::uint64_t>(o.events > 0 ? o.events : o.samples);
    const std::uint64_t base =
        (o.seed >= 0) ? static_cast<std::uint64_t>(o.seed) : std::random_device{}();
    /**
     * ONE STREAM PER ACTIVITY, as the flat engine now does per (node, class).
     *
     * The layered reference indexes its generators by the LQN entity rather
     * than by a station, so the offsets here are built the same way from the
     * activity index: host demands in the service band (+1000) and think times
     * in the arrival band, which is what those two are. Sharing one stream
     * across every activity interleaves draws that the reference keeps apart,
     * and the divergence is immediate rather than statistical.
     */
    const long long ln_seed = static_cast<long long>(base);
    std::vector<Rng> g_host, g_think;
    g_host.reserve(nidx + 1);
    g_think.reserve(nidx + 1);
    for (std::size_t k = 0; k <= nidx; ++k) {
        g_host.push_back(Rng(ln_seed, static_cast<long long>(k) * 10 + 1000));
        g_think.push_back(Rng(ln_seed, static_cast<long long>(k) * 10));
    }
    Rng g_call(ln_seed, 900000);
 
    // ---- resolved model -----------------------------------------------------
    std::vector<Sampler> hostdem(nidx + 1), think(nidx + 1), actthink(nidx + 1);
    std::vector<bool> has_hostdem(nidx + 1, false), has_think(nidx + 1, false),
        has_actthink(nidx + 1, false);
    // ACTIVITY think time (`Activity.setThinkTime`, .lqnx `think-time` on an activity): an
    // off-processor delay after the demand and before the calls, on its own offset block (6000).
    std::vector<Rng> g_actthink;
    g_actthink.reserve(nidx + 1);
    for (std::size_t k = 0; k <= nidx; ++k)
        g_actthink.push_back(Rng(ln_seed, static_cast<long long>(k) * 10 + 6000));
    for (std::size_t a = lsn.ashift + 1; a <= lsn.ashift + lsn.nacts && a < lsn.actthink.size(); ++a) {
        if (!lsn.actthink[a].disabled && num_traits<T>::to_double(lsn.actthink[a].mean) > 0.0) {
            actthink[a] = Sampler(lsn.actthink[a], "the think time of activity '" + lsn.names[a] + "'");
            has_actthink[a] = true;
        }
    }
    for (std::size_t k = 1; k <= nidx; ++k) {
        if (k < lsn.hostdem.size() && !lsn.hostdem[k].disabled &&
            num_traits<T>::to_double(lsn.hostdem[k].mean) > 0.0) {
            hostdem[k] = Sampler(lsn.hostdem[k], "the host demand of '" + lsn.names[k] + "'");
            has_hostdem[k] = true;
        }
        if (k < lsn.think.size() && !lsn.think[k].disabled &&
            num_traits<T>::to_double(lsn.think[k].mean) > 0.0) {
            think[k] = Sampler(lsn.think[k], "the think time of '" + lsn.names[k] + "'");
            has_think[k] = true;
        }
    }
 
    // Host servers and task threads: the two resources an activity holds.
    std::vector<std::size_t> host_servers(nidx + 1, 1), task_threads(nidx + 1, 1);
    std::vector<SchedStrategy> host_sched(nidx + 1, SchedStrategy::FCFS);
    std::vector<SchedStrategy> task_sched(nidx + 1, SchedStrategy::FCFS);
    // An INF host or task never queues, whatever multiplicity it declares, which is how
    // Solver_ssj_ln serves it (startHostService / enterTask branch on the discipline alone).
    for (std::size_t h = lsn.hshift + 1; h <= lsn.hshift + lsn.nhosts; ++h) {
        const double m = (h < lsn.mult.size()) ? lsn.mult[h] : 1.0;
        if (h < lsn.sched.size()) host_sched[h] = lsn.sched[h];
        host_servers[h] = (std::isfinite(m) && host_sched[h] != SchedStrategy::INF)
                              ? std::max<std::size_t>(1, static_cast<std::size_t>(m + 0.5))
                              : std::numeric_limits<std::size_t>::max();
    }
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        const double m = (t < lsn.mult.size()) ? lsn.mult[t] : 1.0;
        if (t < lsn.sched.size()) task_sched[t] = lsn.sched[t];
        task_threads[t] = (std::isfinite(m) && task_sched[t] != SchedStrategy::INF)
                              ? static_cast<std::size_t>(m + 0.5)
                              : std::numeric_limits<std::size_t>::max();
    }
 
    /**
     * REPLICATION. A processor or task declared with replication r is r
     * identical copies of itself, and a copy is a server of its own: pooling
     * them into one server of r times the capacity would let one queue absorb
     * what r separate queues cannot. So the host state below is indexed by a
     * SLOT -- the processor plus its replica -- while what the copies share by
     * definition (servers per copy, scheduling) stays indexed by the element.
     * The reported measures are summed over the copies, which keeps throughput
     * conserved across a call and utilization a fraction of the total capacity.
     */
    std::vector<std::size_t> repl(nidx + 1, 1);
    for (std::size_t k = 1; k < lsn.repl.size() && k <= nidx; ++k) {
        const double r = lsn.repl[k];
        repl[k] = (r > 1.0) ? static_cast<std::size_t>(r + 0.5) : 1;
    }
    std::vector<std::size_t> host_slot0(nidx + 1, 0);
    std::size_t nhost_slots = 0;
    for (std::size_t h = lsn.hshift + 1; h <= lsn.hshift + lsn.nhosts; ++h) {
        host_slot0[h] = nhost_slots;
        nhost_slots += repl[h];
    }
    if (nhost_slots == 0) nhost_slots = 1;
 
    /** Slot of the processor replica running replica trep of one of its tasks. */
    auto host_slot = [&](std::size_t host, std::size_t trep) {
        return host_slot0[host] + (trep % repl[host]);
    };
 
    // Task replicas get their own slots for exactly the reason processors do: r
    // copies are r thread pools, not one pool of r times the size.
    std::vector<std::size_t> task_slot0(nidx + 1, 0);
    std::size_t ntask_slots = 0;
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        task_slot0[t] = ntask_slots;
        ntask_slots += repl[t];
    }
    if (ntask_slots == 0) ntask_slots = 1;
    auto task_slot = [&](std::size_t task, std::size_t trep) {
        return task_slot0[task] + (trep % repl[task]);
    };
 
    /**
     * Replica of the callee reached by one call of replica `crep` of the caller.
     * An unset fan-out is the smallest value consistent with
     * repl(caller)*fanout = repl(callee)*fanin, and the caller reaches the block
     * {(i*f+k) mod r}. A call is one indivisible unit of work, so it goes to one
     * member of that block drawn uniformly: over many calls each member carries
     * the 1/f share that LQN2QN splits the call mean into. This is deterministic
     * pairing at f=1 and a uniform spread over every replica at f=r.
     */
    auto callee_replica = [&](std::size_t caller_task, std::size_t crep,
                              std::size_t callee_task) -> std::size_t {
        const std::size_t rb = repl[callee_task];
        if (rb <= 1) return 0;
        std::size_t f = static_cast<std::size_t>(lsn.fanout_at(caller_task, callee_task) + 0.5);
        if (f == 0) {
            const std::size_t ra = repl[caller_task];
            f = (rb > ra) ? std::max<std::size_t>(1, rb / ra) : 1;
        }
        f = std::min(std::max<std::size_t>(1, f), rb);
        const std::size_t k =
            (f > 1) ? static_cast<std::size_t>(uniform01(g_call) * static_cast<double>(f)) : 0;
        return ((crep * f) + std::min(k, f - 1)) % rb;
    };
 
    // ---- the activity graph, as buildActivityGraph builds it -----------------
    const std::size_t a_lo = lsn.ashift + 1, a_hi = lsn.ashift + lsn.nacts;
    auto is_act = [&](std::size_t k) { return k >= a_lo && k <= a_hi; };
    auto wgt = [&](std::size_t i, std::size_t j) {
        return num_traits<T>::to_double(lsn.graph.get(i, j));
    };
    auto pretype = [&](std::size_t a) {
        return a < lsn.actpretype.size() ? lsn.actpretype[a] : lang::PrecedenceType::NONE;
    };
    auto posttype = [&](std::size_t a) {
        return a < lsn.actposttype.size() ? lsn.actposttype[a] : lang::PrecedenceType::NONE;
    };
    auto phase_of = [&](std::size_t a) -> int {
        const std::size_t k = a - lsn.ashift;
        return (k < lsn.actphase.size() && lsn.actphase[k] > 0) ? lsn.actphase[k] : 1;
    };
    // Successors of an activity are the ACTIVITIES it precedes, ascending: the
    // element graph also carries its call edges, which are not precedence.
    std::vector<std::vector<std::size_t> > succ_act(nidx + 1);
    for (std::size_t a = a_lo; a <= a_hi; ++a)
        for (std::size_t s : lsn.graph.succ(a))
            if (is_act(s)) succ_act[a].push_back(s);
    /**
     * A multi-successor activity CHOOSES one branch when it is an OR fork or a
     * loop decision (its own post type, the first test of `buildActivityGraph`)
     * or when any branch weight is a probability (w != 1, the lazy test of
     * `completeActivity`); otherwise it FORKS every branch. The weights are
     * renormalised when they do not sum to one, and the draw is `u <= cum`.
     */
    std::vector<std::vector<double> > or_prob(nidx + 1);
    for (std::size_t a = a_lo; a <= a_hi; ++a) {
        const std::vector<std::size_t>& sc = succ_act[a];
        if (sc.size() < 2) continue;
        bool prob = posttype(a) == lang::PrecedenceType::POST_OR ||
                    posttype(a) == lang::PrecedenceType::POST_LOOP;
        for (std::size_t s : sc) {
            const double w = wgt(a, s);
            if (w != 1.0 && w > 0.0) prob = true;
        }
        if (!prob) continue;
        double tot = 0.0;
        for (std::size_t s : sc) tot += wgt(a, s);
        for (std::size_t s : sc) {
            const double w = wgt(a, s);
            or_prob[a].push_back((tot > 0.0 && tot != 1.0) ? w / tot : w);
        }
    }
    // An AND join waits for as many inputs as there are PRE_AND activities whose
    // first successor it is, joinRequiredCount of the Java engine.
    std::vector<std::size_t> join_required(nidx + 1, 0);
    for (std::size_t a = a_lo; a <= a_hi; ++a)
        if (pretype(a) == lang::PrecedenceType::PRE_AND && !succ_act[a].empty())
            ++join_required[succ_act[a][0]];
    // The Java engine ignores a quorum and waits for every input. Waiting for
    // all of them here would run a different system from the declared one, and
    // running a quorum would diverge from the reference, so it is refused.
    for (std::size_t a = a_lo; a <= a_hi; ++a) {
        const std::size_t q = a < lsn.actquorum.size() ? lsn.actquorum[a] : 0;
        if (q > 0 && join_required[a] > 0 && q < join_required[a])
            throw UnsupportedError("SolverLDES (native LN engine): the AND join into '" +
                                   lsn.names[a] + "' declares a quorum of " + std::to_string(q) +
                                   " of " + std::to_string(join_required[a]) +
                                   " inputs; LDES joins wait for every input");
    }
    /**
     * The activity an entry starts at: the one the entry's own graph edge names,
     * else the first predecessor-free sequential activity of its task, else the
     * first activity, which is `findBoundActivityForEntry`. 0 means none.
     */
    std::vector<std::size_t> bound(nidx + 1, 0);
    for (std::size_t e = lsn.eshift + 1; e <= lsn.eshift + lsn.nentries; ++e) {
        const std::size_t t = lsn.parent[e];
        const std::vector<std::size_t>& acts =
            (t < lsn.actsof.size()) ? lsn.actsof[t] : std::vector<std::size_t>();
        for (std::size_t a : acts)
            if (wgt(e, a) > 0.0) {
                bound[e] = a;
                break;
            }
        if (bound[e] != 0 || acts.empty()) continue;
        for (std::size_t a : acts) {
            const lang::PrecedenceType pt = pretype(a);
            if (pt != lang::PrecedenceType::PRE_SEQ && pt != lang::PrecedenceType::NONE) continue;
            bool haspred = false;
            for (std::size_t o2 : acts)
                if (o2 != a && wgt(o2, a) > 0.0) {
                    haspred = true;
                    break;
                }
            if (!haspred) {
                bound[e] = a;
                break;
            }
        }
        if (bound[e] == 0) bound[e] = acts[0];
    }
 
    // ---- cache tasks: one content per replica, one access per ItemEntry -------
    std::vector<LnCacheState> caches;
    std::vector<std::size_t> cache_of_slot(ntask_slots, LN_NONE);
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        if (t >= lsn.iscache.size() || !lsn.iscache[t]) continue;
        const std::size_t n = t < lsn.nitems.size() ? lsn.nitems[t] : 0;
        if (n == 0) continue;
        std::vector<int> cap = (t < lsn.itemcap.size()) ? lsn.itemcap[t] : std::vector<int>();
        if (cap.empty()) cap.assign(1, 1);
        for (std::size_t m = 0; m < repl[t]; ++m) {
            LnCacheState cs;
            cs.nitems = n;
            cs.cap = cap;
            cs.rs = (t < lsn.replacestrat.size()) ? lsn.replacestrat[t]
                                                  : lang::ReplacementStrategy::LRU;
            cs.levels.assign(cap.size(), std::list<std::size_t>());
            cs.retrieval = t < lsn.hasretrieval.size() && lsn.hasretrieval[t];
            if (cs.retrieval) {
                cs.inflight.assign(n, 0);
                cs.held.assign(n, {});
            }
            cache_of_slot[task_slot(t, m)] = caches.size();
            caches.push_back(cs);
        }
    }
    std::vector<LnCacheAccess> accesses;
    std::vector<std::size_t> access_of_driver(nidx + 1, LN_NONE);
    std::vector<std::size_t> branch_reply(nidx + 1, 0);  ///< hit/miss activity -> its ItemEntry
    for (std::size_t e = lsn.eshift + 1; e <= lsn.eshift + lsn.nentries; ++e) {
        const std::size_t t = lsn.parent[e];
        if (t <= lsn.tshift || t > lsn.tshift + lsn.ntasks) continue;
        if (cache_of_slot[task_slot(t, 0)] == LN_NONE) continue;
        const std::size_t drv = bound[e];
        const bool is_item = e < lsn.itemproc.size() && !lsn.itemproc[e].empty();
        // An entry of a cache task that is not an ItemEntry runs as an ordinary
        // entry, as it does in the Java engine.
        if (!is_item) continue;
        if (drv == 0 || succ_act[drv].size() < 2)
            throw UnsupportedError("SolverLDES (native LN engine): ItemEntry '" + lsn.names[e] +
                                   "' has no cache access with a hit and a miss branch");
        LnCacheAccess ca;
        for (std::size_t m = 0; m < repl[t]; ++m) ca.caches.push_back(cache_of_slot[task_slot(t, m)]);
        ca.hit = succ_act[drv][0];
        ca.miss = succ_act[drv][1];
        ca.entry = e;
        std::size_t card = e < lsn.nitems.size() ? lsn.nitems[e] : 0;
        if (card == 0) card = caches[ca.caches[0]].nitems;
        double sum = 0.0;
        for (std::size_t i = 0; i < card; ++i) {
            const double v = (i < lsn.itemproc[e].size())
                                 ? num_traits<T>::to_double(lsn.itemproc[e][i]) : 0.0;
            if (std::isnan(v) || v < 0.0)
                throw InputError("SolverLDES (native LN engine): ItemEntry '" + lsn.names[e] +
                                 "' has a popularity that is not a pmf");
            sum += v;
            ca.cdf.push_back(sum);
        }
        if (!(sum > 0.0))
            throw InputError("SolverLDES (native LN engine): ItemEntry '" + lsn.names[e] +
                             "' has a popularity of zero total mass");
        for (double& c : ca.cdf) c /= sum;
        access_of_driver[drv] = accesses.size();
        branch_reply[ca.hit] = e;
        branch_reply[ca.miss] = e;
        accesses.push_back(ca);
    }
    // A POST_CACHE branch with no resolvable read would otherwise run as an AND
    // fork of the hit AND the miss branch, which is not a cache at all.
    for (std::size_t a = a_lo; a <= a_hi; ++a) {
        if (access_of_driver[a] != LN_NONE) continue;
        for (std::size_t s : succ_act[a])
            if (posttype(s) == lang::PrecedenceType::POST_CACHE)
                throw UnsupportedError("SolverLDES (native LN engine): activity '" + lsn.names[a] +
                                       "' reads a cache, but no ItemEntry of a cache task with "
                                       "a popularity is bound to it");
    }
 
    // ---- admission constraints A n <= b --------------------------------------
    /**
     * The rows are declared once per element and hold separately inside each
     * replica, so A and b are per element and the occupancy per slot. A host's
     * columns are its tasks (`tasksof`), a task's its entries (`entriesof`).
     */
    struct Lincon {
        std::vector<std::vector<double> > A;
        std::vector<double> b;
        std::map<std::size_t, int> col;
    };
    std::vector<Lincon> lincon(nidx + 1);
    std::vector<bool> has_lincon(nidx + 1, false);
    for (std::size_t k = 1; k < lsn.lincon_A.size() && k <= nidx; ++k) {
        const Matrix<T>& A = lsn.lincon_A[k];
        if (A.rows() == 0) continue;
        const bool host = k > lsn.hshift && k <= lsn.hshift + lsn.nhosts;
        const std::vector<std::size_t>* ops = nullptr;
        if (host && k < lsn.tasksof.size()) ops = &lsn.tasksof[k];
        if (!host && k < lsn.entriesof.size()) ops = &lsn.entriesof[k];
        const std::size_t ncols = ops ? ops->size() : 0;
        if (A.cols() != ncols)
            throw InputError("SolverLDES (native LN engine): the admission constraint on '" +
                             lsn.names[k] + "' has " + std::to_string(A.cols()) +
                             " columns but the server has " + std::to_string(ncols) + " operands");
        if (ncols == 0) continue;
        Lincon lc;
        lc.A.assign(A.rows(), std::vector<double>(ncols, 0.0));
        for (std::size_t r = 0; r < A.rows(); ++r)
            for (std::size_t c = 0; c < ncols; ++c) lc.A[r][c] = num_traits<T>::to_double(A(r, c));
        for (std::size_t r = 0; r < A.rows(); ++r)
            lc.b.push_back((k < lsn.lincon_b.size() && r < lsn.lincon_b[k].size())
                               ? num_traits<T>::to_double(lsn.lincon_b[k][r])
                               : std::numeric_limits<double>::infinity());
        for (std::size_t c = 0; c < ncols; ++c) lc.col[(*ops)[c]] = static_cast<int>(c);
        lincon[k] = lc;
        has_lincon[k] = true;
    }
    auto lincon_col = [&](std::size_t elem, std::size_t operand) -> int {
        if (!has_lincon[elem]) return -1;
        const std::map<std::size_t, int>::const_iterator it = lincon[elem].col.find(operand);
        return it == lincon[elem].col.end() ? -1 : it->second;
    };
    /** True when one more job in column `col` still satisfies every row. */
    auto admits = [&](std::size_t elem, const std::vector<int>& occ, int col) {
        const Lincon& lc = lincon[elem];
        for (std::size_t r = 0; r < lc.A.size(); ++r) {
            double lhs = lc.A[r][static_cast<std::size_t>(col)];
            for (std::size_t k = 0; k < occ.size(); ++k) lhs += lc.A[r][k] * occ[k];
            if (lhs > lc.b[r] + 1e-9) return false;
        }
        return true;
    };
 
    // ---- heterogeneous server pools on processors ----------------------------
    /**
     * The pools ARE the servers of the host, so their counts must add up to its
     * multiplicity. An INF host ignores them, as the Java engine's service path
     * does. The operand column of a job is its task's position in `tasksof`.
     */
    std::vector<bool> pooled(nidx + 1, false);
    std::vector<std::vector<double> > pool_capacity(nidx + 1);
    std::vector<std::vector<std::vector<bool> > > pool_serves(nidx + 1);
    std::vector<std::map<std::size_t, int> > pool_col(nidx + 1);
    for (std::size_t h = lsn.hshift + 1; h <= lsn.hshift + lsn.nhosts; ++h) {
        if (h >= lsn.pools.size() || lsn.pools[h].npools() == 0) continue;
        if (host_sched[h] == SchedStrategy::INF) continue;  // the INF path ignores pools
        const lqn::ServerPools<T>& P = lsn.pools[h];
        const std::size_t np = P.npools(), nc = P.compat.cols();
        std::size_t declared = 0;
        for (std::size_t p = 0; p < np; ++p) {
            const double cnt = (p < P.counts.size()) ? P.counts[p] : 0.0;
            const std::size_t c = cnt > 0.0 ? static_cast<std::size_t>(cnt) : 0;
            const double r = (p < P.rates.size()) ? num_traits<T>::to_double(P.rates[p]) : 1.0;
            pool_capacity[h].push_back(static_cast<double>(c) * (r > 0.0 ? r : 1.0));
            declared += c;
            std::vector<bool> sv(nc, false);
            for (std::size_t j = 0; j < nc; ++j)
                sv[j] = num_traits<T>::to_double(P.compat(p, j)) != 0.0;
            pool_serves[h].push_back(sv);
        }
        if (declared != host_servers[h])
            throw InputError("SolverLDES (native LN engine): processor '" + lsn.names[h] +
                             "' declares server pools holding " + std::to_string(declared) +
                             " servers but has multiplicity " +
                             (host_servers[h] == std::numeric_limits<std::size_t>::max()
                                  ? std::string("Inf") : std::to_string(host_servers[h])) +
                             "; the pools must partition the servers");
        const std::vector<std::size_t>& ops =
            (h < lsn.tasksof.size()) ? lsn.tasksof[h] : std::vector<std::size_t>();
        for (std::size_t j = 0; j < ops.size() && j < nc; ++j) pool_col[h][ops[j]] = static_cast<int>(j);
        pooled[h] = true;
    }
 
    // ---- live state ---------------------------------------------------------
    double now = 0.0;
    std::priority_queue<LnEvent, std::vector<LnEvent>, LnEventLater> evq;
    std::uint64_t seq = 0;
    // Lines and invocations live in deques, whose elements do not move when
    // another is added, and are recycled through free lists.
    std::deque<LnLine> lines;
    std::deque<LnInv> invs;
    std::vector<std::size_t> free_lines, free_invs;
    std::vector<LnCustomer> custs;
 
    // Busy servers and waiting lines of ONE processor replica, addressed by slot.
    std::vector<std::size_t> host_busy(nhost_slots, 0);
    std::vector<std::deque<std::size_t>> host_queue(nhost_slots);
    // A pooled replica holds its lines in arrival order, the order the pools scan.
    std::vector<std::vector<std::size_t> > pool_jobs(nhost_slots);
    std::vector<double> pool_last(nhost_slots, 0.0), pool_rate(nhost_slots, 0.0);
    std::vector<std::uint64_t> pool_gen(nhost_slots, 0);
    std::vector<std::size_t> host_of_slot(nhost_slots, 0);
    for (std::size_t h = lsn.hshift + 1; h <= lsn.hshift + lsn.nhosts; ++h)
        for (std::size_t m = 0; m < repl[h]; ++m) host_of_slot[host_slot0[h] + m] = h;
    // Admission occupancy and blocked requests, per slot.
    std::vector<std::vector<int> > hocc(nhost_slots), tocc(ntask_slots);
    std::vector<std::deque<std::size_t> > hadm(nhost_slots), tadm(ntask_slots);
    for (std::size_t h = lsn.hshift + 1; h <= lsn.hshift + lsn.nhosts; ++h)
        if (has_lincon[h])
            for (std::size_t m = 0; m < repl[h]; ++m)
                hocc[host_slot0[h] + m].assign(lincon[h].A[0].size(), 0);
 
    // Task threads of ONE task replica, addressed by slot. An invocation that
    // finds every thread taken waits in thr_queue, holding no resource meanwhile.
    std::vector<std::vector<bool>> thr_held(ntask_slots);
    std::vector<std::vector<ThreadState>> thr_state(ntask_slots);
    std::vector<std::vector<double>> thr_setup_t0(ntask_slots);
    std::vector<std::vector<std::uint64_t>> thr_gen(ntask_slots);
    std::vector<std::deque<std::size_t>> thr_queue(ntask_slots);
    std::vector<std::size_t> thr_owner(ntask_slots, 0);   ///< task index of each slot
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        const std::size_t n = (task_threads[t] == std::numeric_limits<std::size_t>::max())
                                  ? 0 : task_threads[t];
        for (std::size_t m = 0; m < repl[t]; ++m) {
            const std::size_t ts = task_slot(t, m);
            thr_owner[ts] = t;
            thr_held[ts].assign(n, false);
            thr_state[ts].assign(n, ThreadState::ACTIVE);
            thr_setup_t0[ts].assign(n, 0.0);
            thr_gen[ts].assign(n, 0);
            if (has_lincon[t]) tocc[ts].assign(lincon[t].A[0].size(), 0);
        }
    }
 
    // A SetupTask's two clocks, per task. Both must be positive for the power
    // cycle to exist: a setup with no delay-off would never fire, because nothing
    // would ever power a thread down.
    std::vector<Sampler> setupd(nidx + 1), delayoffd(nidx + 1);
    std::vector<bool> has_setup(nidx + 1, false);
    std::vector<Rng> g_setup, g_doff;
    g_setup.reserve(nidx + 1);
    g_doff.reserve(nidx + 1);
    for (std::size_t k = 0; k <= nidx; ++k) {
        g_setup.push_back(Rng(ln_seed, static_cast<long long>(k) * 10 + 2000));
        g_doff.push_back(Rng(ln_seed, static_cast<long long>(k) * 10 + 3000));
    }
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        if (t >= lsn.hassetup.size() || !lsn.hassetup[t]) continue;
        const bool su = t < lsn.setuptime.size() && !lsn.setuptime[t].disabled &&
                        num_traits<T>::to_double(lsn.setuptime[t].mean) > 0.0;
        const bool df = t < lsn.delayofftime.size() && !lsn.delayofftime[t].disabled &&
                        num_traits<T>::to_double(lsn.delayofftime[t].mean) > 0.0;
        if (!su || !df) continue;
        setupd[t] = Sampler(lsn.setuptime[t], "the setup time of '" + lsn.names[t] + "'");
        delayoffd[t] = Sampler(lsn.delayofftime[t], "the delay-off time of '" + lsn.names[t] + "'");
        has_setup[t] = true;
    }
 
    /**
     * The new features draw on streams of their own, in blocks above every
     * existing offset, so that a model without them keeps its seeded sample
     * path: branch choices per activity (940000), item draws per cache replica
     * (950000) and random replacement per cache replica (960000). 910000 and
     * 920000 are the SIRO keys of hosts and tasks.
     */
    std::vector<Rng> g_route, g_item, g_rr;
    g_route.reserve(nidx + 1);
    for (std::size_t k = 0; k <= nidx; ++k)
        g_route.push_back(Rng(ln_seed, 940000 + static_cast<long long>(k) * 10));
    for (std::size_t c = 0; c < caches.size(); ++c) {
        g_item.push_back(Rng(ln_seed, 950000 + static_cast<long long>(c) * 10));
        g_rr.push_back(Rng(ln_seed, 960000 + static_cast<long long>(c) * 10));
    }
 
    /// Lines woken by a release, continued after the current walk returns.
    /// Waking them inside the walk would re-enter it on another line.
    std::deque<std::size_t> runnable;
 
    /**
     * PRIORITY ORDER, as Solver_ssj_ln orders it: a LARGER `lsn.prio` first (the lqns convention),
     * FIFO within one level under HOL/FCFSPRIO and the FCFS preemptive variants, LIFO under the
     * LCFS ones. A host request carries the priority of the task whose activity it is; a task
     * request that of its CALLER, since every request at one task belongs to that task, and one
     * with no caller (a reference task's own cycle, an open arrival, an asynchronous call) takes
     * 0, the lowest.
     */
    auto task_prio = [&](std::size_t task) -> int {
        return (task < lsn.prio.size()) ? lsn.prio[task] : 0;
    };
    auto host_req_prio = [&](std::size_t lid) -> int { return task_prio(lsn.parent[lines[lid].act]); };
    auto task_req_prio = [&](std::size_t vi) -> int {
        const std::size_t c = invs[vi].caller;
        return c == LN_NONE ? 0 : task_prio(invs[lines[c].inv].task);
    };
    // SIRO keys, one stream per resource kind on the engine's own offsets, so a SIRO queue does
    // not shift the demand, think or call variates of any other element.
    Rng g_siro_host(ln_seed, 910000), g_siro_task(ln_seed, 920000);
    auto fcfs_prio = [](SchedStrategy sc) {
        // FCFSPRIO has no enumerator of its own here: it is MATLAB's alias for HOL.
        return sc == SchedStrategy::HOL || sc == SchedStrategy::FCFSPRPRIO ||
               sc == SchedStrategy::FCFSPIPRIO;
    };
    auto lcfs_prio = [](SchedStrategy sc) {
        return sc == SchedStrategy::LCFSPRIO || sc == SchedStrategy::LCFSPRPRIO ||
               sc == SchedStrategy::LCFSPIPRIO;
    };
    auto preemptive = [](SchedStrategy sc) {
        return sc == SchedStrategy::FCFSPRPRIO || sc == SchedStrategy::FCFSPIPRIO ||
               sc == SchedStrategy::LCFSPRPRIO || sc == SchedStrategy::LCFSPIPRIO;
    };
    auto is_ps = [](SchedStrategy sc) {
        return sc == SchedStrategy::PS || sc == SchedStrategy::PSPRIO;
    };
    /**
     * Queue `id` behind every request of a higher priority and, within its own level, behind the
     * ones FIFO (or, `lifo`, LIFO) puts first. `key_of` is the arrival instant a level is ordered
     * by; a task queue passes none and relies on requests joining it in arrival order.
     */
    auto prio_enqueue = [](std::deque<std::size_t>& q, std::size_t id,
                           const std::function<int(std::size_t)>& prio_of, bool lifo,
                           const std::function<double(std::size_t)>* key_of) {
        const int prio = prio_of(id);
        const double key = key_of ? (*key_of)(id) : 0.0;
        std::deque<std::size_t>::iterator it = q.end();
        while (it != q.begin()) {
            const std::size_t o = *(it - 1);
            const int po = prio_of(o);
            if (po > prio) break;
            if (po == prio) {
                if (!lifo && (!key_of || (*key_of)(o) <= key)) break;
                if (lifo && key_of && (*key_of)(o) > key) break;
            }
            --it;
        }
        q.insert(it, id);
    };
    /**
     * Queue `id` under `sc`, the order Solver_ssj_ln's comparators impose: FCFS by arrival, LCFS
     * newest first, SIRO by a uniform key drawn now, the priority disciplines by `prio` then
     * arrival (`key_of`, the host residence start, when given).
     */
    auto enqueue = [&](std::deque<std::size_t>& q, std::size_t id, SchedStrategy sc,
                       const std::function<int(std::size_t)>& prio_of, double& key, Rng& g_siro,
                       const std::function<double(std::size_t)>* key_of) {
        if (fcfs_prio(sc) || lcfs_prio(sc)) {
            prio_enqueue(q, id, prio_of, lcfs_prio(sc), key_of);
        } else if (sc == SchedStrategy::LCFS) {
            q.push_front(id);
        } else {
            if (sc == SchedStrategy::SIRO) key = uniform01(g_siro);
            q.push_back(id);
        }
    };
    /** Take the request `sc` serves next: the head, or under SIRO the smallest key. */
    auto dequeue = [&](std::deque<std::size_t>& q, SchedStrategy sc,
                       const std::function<double(std::size_t)>& key_of) -> std::size_t {
        std::deque<std::size_t>::iterator best = q.begin();
        if (sc == SchedStrategy::SIRO)
            for (std::deque<std::size_t>::iterator it = q.begin(); it != q.end(); ++it)
                if (key_of(*it) < key_of(*best)) best = it;
        const std::size_t id = *best;
        q.erase(best);
        return id;
    };
    const std::function<double(std::size_t)> line_key = [&](std::size_t lid) {
        return lines[lid].siro_key;
    };
    const std::function<double(std::size_t)> inv_key = [&](std::size_t vi) {
        return invs[vi].siro_key;
    };
    const std::function<double(std::size_t)> host_key = [&](std::size_t lid) {
        return lines[lid].host_t0;
    };
 
    /**
     * PROCESSOR SHARING on a host replica, Solver_ssj_ln's startPSJob / advancePSService: every line
     * present is served at min(1, c/n) of a server, and the next completion is re-read whenever the
     * set changes. The event queue has no removal, so a stale completion is recognised by its
     * generation (event kind 5).
     */
    std::vector<std::vector<std::size_t> > ps_jobs(nhost_slots);
    std::vector<double> ps_last(nhost_slots, 0.0);
    std::vector<std::uint64_t> ps_gen(nhost_slots, 0);
    /**
     * The share of each line of slot `hs`: min(1, c/n) under PS; under PSPRIO, Solver_ssj_ln's
     * psprioJobRates, the levels served from the largest task priority down, each taking
     * min(servers left, its lines) servers split evenly. Either way the shares sum to min(n, c).
     */
    auto ps_rates = [&](std::size_t hs) -> std::vector<double> {
        const std::vector<std::size_t>& pj = ps_jobs[hs];
        const std::size_t n = pj.size();
        const std::size_t host = host_of_slot[hs];
        const double c = static_cast<double>(host_servers[host]);
        if (host_sched[host] != SchedStrategy::PSPRIO || static_cast<double>(n) <= c)
            return std::vector<double>(n, n == 0 ? 0.0 : std::min(1.0, c / static_cast<double>(n)));
        std::vector<double> rate(n, 0.0);
        std::vector<int> pr(n);
        for (std::size_t i = 0; i < n; ++i) pr[i] = task_prio(lsn.parent[lines[pj[i]].act]);
        std::vector<int> levels(pr);
        std::sort(levels.begin(), levels.end(), std::greater<int>());
        levels.erase(std::unique(levels.begin(), levels.end()), levels.end());
        double left = c;
        for (std::size_t l = 0; l < levels.size() && left > 0.0; ++l) {
            std::size_t cnt = 0;
            for (std::size_t i = 0; i < n; ++i) cnt += (pr[i] == levels[l]) ? 1 : 0;
            const double alloc = std::min(left, static_cast<double>(cnt));
            for (std::size_t i = 0; i < n; ++i)
                if (pr[i] == levels[l]) rate[i] = alloc / static_cast<double>(cnt);
            left -= alloc;
        }
        return rate;
    };
    /** Serve every line of slot `hs` up to now at the shares in force since the last change. */
    auto ps_advance = [&](std::size_t hs) {
        const double dt = now - ps_last[hs];
        if (dt > 0.0) {
            const std::vector<double> rate = ps_rates(hs);
            for (std::size_t i = 0; i < ps_jobs[hs].size(); ++i) {
                LnLine& L = lines[ps_jobs[hs][i]];
                L.remaining = std::max(0.0, L.remaining - dt * rate[i]);
            }
        }
        ps_last[hs] = now;
    };
    /** Busy servers of a PS slot holding n lines, the min(n, c) Solver_ssj_ln.updateHostStats counts. */
    auto ps_busy = [&](std::size_t hs, std::size_t n) -> double {
        return static_cast<double>(std::min(n, host_servers[host_of_slot[hs]]));
    };
 
    // Time integrals per element.
    std::vector<double> tot_q(nidx + 1, 0.0), tot_u(nidx + 1, 0.0);
    std::vector<double> cur_q(nidx + 1, 0.0), cur_u(nidx + 1, 0.0);
    std::vector<double> last_upd(nidx + 1, 0.0);
    std::vector<double> completions(nidx + 1, 0.0), resp_sum(nidx + 1, 0.0),
        resp_cnt(nidx + 1, 0.0);
    // The observations behind `resp_sum` at the ENTRY indices, kept for
    // get_cdf_respt. Only entries: an activity's and a task's response times are
    // means the layered table reports, with no distributional getter on them.
    std::vector<std::vector<double> > entry_resp_samples(lsn.nentries);
    // Total time the executions of an activity spend AT ITS HOST PROCESSOR --
    // queueing there plus running. Over the elapsed time this is the mean
    // occupancy at the host, which is what ResidT is built from below.
    // MEASURED rather than derived from the nominal `hostdem` mean: that mean is
    // the SERVICE time, and on a contended processor the residence is that plus
    // the wait, so a derived figure would report an uncontended residence at
    // exactly the models where it matters.
    std::vector<double> act_host_resid(nidx + 1, 0.0);
 
    // ONE STREAM PER CALL for the repetition count, on the engine's own offset
    // block (4000), so that adding a call does not shift the variates of the
    // think times or the host demands and re-baseline every seeded result.
    std::vector<Rng> g_callmult;
    g_callmult.reserve(lsn.ncalls + 1);
    for (std::size_t c = 0; c <= lsn.ncalls; ++c)
        g_callmult.push_back(Rng(ln_seed, static_cast<long long>(c) * 10 + 4000));
 
    /**
     * How many times one execution of the caller makes call `cidx`.
     *
     * `floor(mean)` calls plus one more with probability equal to the fraction,
     * which is `Solver_ssj_ln.sampleCallCount` exactly. It is DETERMINISTIC at
     * an integer mean, which is what `calls-mean="3"` asks for: three calls
     * every time, not three on average.
     */
    auto sample_call_count = [&](std::size_t cidx) -> std::size_t {
        const double m = (cidx < lsn.callproc_mean.size())
                             ? num_traits<T>::to_double(lsn.callproc_mean[cidx])
                             : 1.0;
        if (!(m > 0.0)) return 0;
        std::size_t n = static_cast<std::size_t>(std::floor(m));
        const double frac = m - static_cast<double>(n);
        if (frac > 0.0 && cidx < g_callmult.size() && uniform01(g_callmult[cidx]) < frac) ++n;
        return n;
    };
 
    auto touch = [&](std::size_t k) {
        const double dt = now - last_upd[k];
        if (dt > 0.0) {
            tot_q[k] += cur_q[k] * dt;
            tot_u[k] += cur_u[k] * dt;
        }
        last_upd[k] = now;
    };
    auto push_ev = [&](LnEvent e) {
        e.seq = seq++;
        evq.push(e);
    };
 
    std::uint64_t done = 0;
 
    // ---- lines and invocations ------------------------------------------------
    auto new_line = [&](std::size_t inv, std::size_t act, bool branch) -> std::size_t {
        std::size_t id;
        if (!free_lines.empty()) {
            id = free_lines.back();
            free_lines.pop_back();
            lines[id] = LnLine();
        } else {
            id = lines.size();
            lines.push_back(LnLine());
        }
        lines[id].inv = inv;
        lines[id].act = act;
        lines[id].calls_left = NOTDRAWN;
        lines[id].branch = branch;
        return id;
    };
    auto free_line = [&](std::size_t id) { free_lines.push_back(id); };
    /** Move a line to activity `a`, not yet begun. */
    auto set_act = [&](std::size_t lid, std::size_t a) {
        LnLine& L = lines[lid];
        L.act = a;
        L.act_t0 = -1.0;
        L.call_pos = 0;
        L.calls_left = NOTDRAWN;
        L.dem_done = 0;
    };
    auto new_inv =[&](std::size_t entry, std::size_t trep, std::size_t caller,
                       std::size_t cust) -> std::size_t {
        std::size_t id;
        if (!free_invs.empty()) {
            id = free_invs.back();
            free_invs.pop_back();
            invs[id] = LnInv();
        } else {
            id = invs.size();
            invs.push_back(LnInv());
        }
        LnInv& v = invs[id];
        v.entry = entry;
        v.task = lsn.parent[entry];
        v.trep = trep;
        v.ts = task_slot(v.task, trep);
        v.th = NOTHR;
        v.caller = caller;
        v.cust = cust;
        const std::size_t root = new_line(id, bound[entry], false);
        invs[id].root = root;
        return id;
    };
 
    /**
     * Enter the task: the invocation now counts at the task. It counts at its ENTRY only once it
     * holds a thread (see `advance`), since the entry's service starts there.
     */
    auto enter_task = [&](std::size_t vi) {
        LnInv& v = invs[vi];
        touch(v.task);
        cur_q[v.task] += 1.0;
    };
    /**
     * Arrive at the task, through its admission constraint when it has one. A
     * blocked request holds neither a thread nor a place in the task's queue,
     * so it is counted nowhere while it waits, as `arriveAtTask` has it.
     */
    auto arrive_at_task = [&](std::size_t vi) -> bool {
        LnInv& v = invs[vi];
        const int col = lincon_col(v.task, v.entry);
        if (col >= 0) {
            v.tcol = col;
            if (!admits(v.task, tocc[v.ts], col)) {
                tadm[v.ts].push_back(vi);
                return false;
            }
            ++tocc[v.ts][static_cast<std::size_t>(col)];
        }
        enter_task(vi);
        return true;
    };
 
    // ---- the activity walk --------------------------------------------------
    std::function<void(std::size_t)> advance;
    /** Continue every line a release woke, and every line those wake. */
    auto drain = [&]() {
        while (!runnable.empty()) {
            const std::size_t id = runnable.front();
            runnable.pop_front();
            advance(id);
        }
    };
 
    /** Put line `lid` into service at its host replica. */
    std::function<void(std::size_t)> advance_pool;
    auto pool_rates = [&](std::size_t hs) {
        const std::size_t h = host_of_slot[hs];
        const std::vector<std::size_t>& jobs = pool_jobs[hs];
        std::vector<double> rate(jobs.size(), 0.0);
        for (std::size_t p = 0; p < pool_capacity[h].size(); ++p)
            for (std::size_t i = 0; i < jobs.size(); ++i) {
                const int j = lines[jobs[i]].pcol;
                if (j >= 0 && static_cast<std::size_t>(j) < pool_serves[h][p].size() &&
                    pool_serves[h][p][static_cast<std::size_t>(j)]) {
                    rate[i] += pool_capacity[h][p];
                    break;
                }
            }
        return rate;
    };
    /** Serve every job of a pooled replica up to now at the rates in force. */
    advance_pool = [&](std::size_t hs) {
        const std::vector<double> r = pool_rates(hs);
        const double dt = now - pool_last[hs];
        for (std::size_t i = 0; i < pool_jobs[hs].size(); ++i) {
            LnLine& L = lines[pool_jobs[hs][i]];
            L.remaining = std::max(0.0, L.remaining - dt * r[i]);
        }
        pool_last[hs] = now;
    };
    /**
     * Re-read the delivered rate after the job set changed, and arm the next
     * completion. A pooled host's Util integrates the DELIVERED rate, which is
     * `hostBusyTime += delivered*dt` of `advancePSService`.
     */
    auto rearm_pool = [&](std::size_t hs) {
        const std::size_t h = host_of_slot[hs];
        const std::vector<double> r = pool_rates(hs);
        double tot = 0.0, tmin = std::numeric_limits<double>::infinity();
        for (std::size_t i = 0; i < r.size(); ++i) {
            tot += r[i];
            if (r[i] > 0.0) tmin = std::min(tmin, lines[pool_jobs[hs][i]].remaining / r[i]);
        }
        touch(h);
        pool_rate[hs] = tot;
        double u = 0.0;
        for (std::size_t m = 0; m < repl[h]; ++m) u += pool_rate[host_slot0[h] + m];
        cur_u[h] = u;
        const std::uint64_t g = ++pool_gen[hs];
        if (std::isfinite(tmin)) {
            LnEvent e;
            e.t = now + tmin;
            e.kind = 4;
            e.who = hs;
            e.gen = g;
            push_ev(e);
        }
    };
    /** Arm the next completion of PS slot `hs` under the shares that hold from now on. */
    auto ps_schedule = [&](std::size_t hs) {
        ++ps_gen[hs];
        if (ps_jobs[hs].empty()) return;
        const std::vector<double> rate = ps_rates(hs);
        double tmin = std::numeric_limits<double>::infinity();
        for (std::size_t i = 0; i < ps_jobs[hs].size(); ++i)
            if (rate[i] > 0.0) tmin = std::min(tmin, lines[ps_jobs[hs][i]].remaining / rate[i]);
        LnEvent e;
        e.t = now + tmin;
        e.kind = 5;
        e.who = hs;
        e.gen = ps_gen[hs];
        push_ev(e);
    };
 
    /// Lines in service on a host slot under a preemptive discipline, the candidates to preempt.
    std::vector<std::vector<std::size_t> > in_service(nhost_slots);
    std::uint64_t host_gen_ctr = 0;
    /**
     * Put line `lid` on a server of queueing host slot `hs` for the work it has left, drawing the
     * demand when it has none yet (a fresh request, or a PI victim resuming).
     */
    auto begin_host = [&](std::size_t lid, std::size_t hs) {
        LnLine& L = lines[lid];
        if (L.host_rem < 0.0) L.host_rem = hostdem[L.act].next(g_host[L.act]);
        L.host_start = now;
        L.host_end = now + L.host_rem;
        if (preemptive(host_sched[host_of_slot[hs]])) in_service[hs].push_back(lid);
        L.host_gen = ++host_gen_ctr;
        LnEvent e;
        e.t = L.host_end;
        e.kind = 0;
        e.who = lid;
        e.gen = L.host_gen;
        push_ev(e);
    };
    /**
     * Preempt a line in service on full slot `hs` for `lid`, Solver_ssj_ln.preemptFor: FCFSPR/PIPRIO
     * take only a STRICTLY lower priority, the lowest in service and the latest started among ties;
     * LCFSPR/PIPRIO take the lowest strictly lower one (earliest arrival among ties), else the equal
     * one that started earliest. The victim rejoins the queue at its original arrival, keeping its
     * residual work under PR and drawing afresh under PI. Returns whether `lid` took a server.
     */
    auto preempt_for = [&](std::size_t lid, std::size_t hs) -> bool {
        const SchedStrategy sc = host_sched[host_of_slot[hs]];
        const bool lifo = lcfs_prio(sc);
        const int pa = host_req_prio(lid);
        std::vector<std::size_t>& ins = in_service[hs];
        std::size_t vi = ins.size();
        int pv = 0;
        for (std::size_t i = 0; i < ins.size(); ++i) {
            const int pr = host_req_prio(ins[i]);
            if (pr >= pa) continue;
            bool better = vi == ins.size() || pr < pv;
            if (!better && pr == pv) {
                const LnLine& r = lines[ins[i]];
                const LnLine& v = lines[ins[vi]];
                better = lifo ? r.host_t0 < v.host_t0 : r.host_start > v.host_start;
            }
            if (better) {
                vi = i;
                pv = pr;
            }
        }
        if (vi == ins.size() && lifo)
            for (std::size_t i = 0; i < ins.size(); ++i)
                if (host_req_prio(ins[i]) == pa &&
                    (vi == ins.size() || lines[ins[i]].host_start < lines[ins[vi]].host_start))
                    vi = i;
        if (vi == ins.size()) return false;
        const std::size_t victim = ins[vi];
        ins.erase(ins.begin() + static_cast<std::ptrdiff_t>(vi));
        LnLine& v = lines[victim];
        v.host_gen = ++host_gen_ctr;  // its pending completion is now stale
        const bool resume = sc == SchedStrategy::FCFSPRPRIO || sc == SchedStrategy::LCFSPRPRIO;
        v.host_rem = resume ? std::max(0.0, v.host_end - now) : -1.0;
        enqueue(host_queue[hs], victim, sc, host_req_prio, v.siro_key, g_siro_host, &host_key);
        begin_host(lid, hs);
        return true;
    };
 
    auto start_host_service = [&](std::size_t lid) {
        LnLine& L = lines[lid];
        const std::size_t act = L.act;
        const std::size_t host = lsn.host_of(act);
        const std::size_t hs = host_slot(host, invs[L.inv].trep);
        touch(host);
        cur_q[host] += 1.0;
        if (pooled[host]) {
            advance_pool(hs);
            const std::size_t t = lsn.parent[act];
            const std::map<std::size_t, int>::const_iterator it = pool_col[host].find(t);
            if (it == pool_col[host].end())
                throw InputError("SolverLDES (native LN engine): activity '" + lsn.names[act] +
                                 "' runs on a host whose pool declaration does not name its "
                                 "task, so no server can be assigned");
            lines[lid].pcol = it->second;
            lines[lid].remaining = hostdem[act].next(g_host[act]);
            pool_jobs[hs].push_back(lid);
            rearm_pool(hs);
            return;
        }
        lines[lid].host_rem = -1.0;
        if (is_ps(host_sched[host])) {
            // Charge the interval that just ended at the share in force during it, then admit.
            ps_advance(hs);
            const std::size_t n0 = ps_jobs[hs].size();
            lines[lid].remaining = hostdem[act].next(g_host[act]);
            ps_jobs[hs].push_back(lid);
            ++host_busy[hs];
            cur_u[host] += ps_busy(hs, n0 + 1) - ps_busy(hs, n0);
            ps_schedule(hs);
            return;
        }
        if (host_busy[hs] < host_servers[host]) {
            ++host_busy[hs];
            cur_u[host] += 1.0;
            begin_host(lid, hs);
        } else if (!(preemptive(host_sched[host]) && preempt_for(lid, hs))) {
            enqueue(host_queue[hs], lid, host_sched[host], host_req_prio, lines[lid].siro_key,
                    g_siro_host, &host_key);
        }
    };
    /**
     * Ask the host for service. The residence starts HERE, when the line asks,
     * whether it then queues at the processor or waits outside it for a host
     * admission token; a blocked line keeps its thread, being still inside
     * its activity, but holds no server.
     */
    auto request_host = [&](std::size_t lid) {
        LnLine& L = lines[lid];
        const std::size_t act = L.act;
        const std::size_t host = lsn.host_of(act);
        const std::size_t hs = host_slot(host, invs[L.inv].trep);
        L.host_t0 = now;
        const int col = lincon_col(host, lsn.parent[act]);
        if (col >= 0) {
            L.hcol = col;
            if (!admits(host, hocc[hs], col)) {
                hadm[hs].push_back(lid);
                return;
            }
            ++hocc[hs][static_cast<std::size_t>(col)];
        }
        start_host_service(lid);
    };
    /** Give back a host admission token and admit the blocked head while it fits. */
    auto release_host_admission = [&](std::size_t lid) {
        LnLine& L = lines[lid];
        if (L.hcol < 0) return;
        const std::size_t host = lsn.host_of(L.act);
        const std::size_t hs = host_slot(host, invs[L.inv].trep);
        --hocc[hs][static_cast<std::size_t>(L.hcol)];
        L.hcol = -1;
        while (!hadm[hs].empty()) {
            const std::size_t head = hadm[hs].front();
            const int c = lines[head].hcol;
            if (!admits(host, hocc[hs], c)) return;
            hadm[hs].pop_front();
            ++hocc[hs][static_cast<std::size_t>(c)];
            start_host_service(head);
        }
    };
 
    /** Arm THREAD's cold start; the request that woke it waits in thr_queue. */
    auto start_setup = [&](std::size_t ts, std::size_t th) {
        const std::size_t task = thr_owner[ts];
        thr_state[ts][th] = ThreadState::SETUP;
        thr_setup_t0[ts][th] = now;
        LnEvent e;
        e.t = now + setupd[task].next(g_setup[task]);
        e.kind = 2;
        e.who = ts;
        e.aux = th;
        e.gen = ++thr_gen[ts][th];
        push_ev(e);
    };
 
    /** Arm THREAD's idle countdown; when it expires the thread is OFF. */
    auto start_delayoff = [&](std::size_t ts, std::size_t th) {
        const std::size_t task = thr_owner[ts];
        thr_state[ts][th] = ThreadState::DELAYOFF;
        LnEvent e;
        e.t = now + delayoffd[task].next(g_doff[task]);
        e.kind = 3;
        e.who = ts;
        e.aux = th;
        e.gen = ++thr_gen[ts][th];
        push_ev(e);
    };
 
    /**
     * Give invocation `vi` a thread of its task.
     *
     * Returns false when it could not proceed: either every thread is taken, or
     * the only thread available was powered off and is now warming up. Either
     * way it is parked in thr_queue and continued by whoever frees or finishes
     * warming a thread. An infinite-thread task never blocks and holds no
     * identified thread.
     */
    auto acquire_thread = [&](std::size_t vi) -> bool {
        LnInv& v = invs[vi];
        const std::size_t task = v.task;
        const std::size_t ts = v.ts;
        touch(task);
        if (task_threads[task] == std::numeric_limits<std::size_t>::max()) {
            cur_u[task] += 1.0;
            v.has_thread = true;
            v.th = NOTHR;
            return true;
        }
        std::vector<bool>& held = thr_held[ts];
        std::vector<ThreadState>& st = thr_state[ts];
        // An ACTIVE idle thread first, then one still counting down -- caught
        // before it powers off, it serves having paid nothing -- and only then a
        // thread that is already OFF, which must warm up first.
        for (std::size_t th = 0; th < held.size(); ++th) {
            if (!held[th] && st[th] == ThreadState::ACTIVE) {
                held[th] = true;
                cur_u[task] += 1.0;
                v.has_thread = true;
                v.th = th;
                return true;
            }
        }
        for (std::size_t th = 0; th < held.size(); ++th) {
            if (!held[th] && st[th] == ThreadState::DELAYOFF) {
                ++thr_gen[ts][th];  // cancel the countdown
                st[th] = ThreadState::ACTIVE;
                held[th] = true;
                cur_u[task] += 1.0;
                v.has_thread = true;
                v.th = th;
                return true;
            }
        }
        enqueue(thr_queue[ts], vi, task_sched[task], task_req_prio, invs[vi].siro_key, g_siro_task,
                nullptr);
        for (std::size_t th = 0; th < held.size(); ++th) {
            if (!held[th] && st[th] == ThreadState::OFF) {
                start_setup(ts, th);
                break;
            }
        }
        return false;
    };
 
    /** Free the thread invocation `vi` holds and hand it to the next one waiting. */
    auto release_thread = [&](std::size_t vi) {
        LnInv& v = invs[vi];
        const std::size_t ts = v.ts;
        const std::size_t th = v.th;
        const std::size_t task = thr_owner[ts];
        touch(task);
        cur_u[task] -= 1.0;
        v.has_thread = false;
        v.th = NOTHR;
        if (th == NOTHR) return;  // an infinite-thread task holds nothing
        thr_held[ts][th] = false;
        if (!thr_queue[ts].empty()) {
            const std::size_t nxt = dequeue(thr_queue[ts], task_sched[task], inv_key);
            thr_held[ts][th] = true;
            cur_u[task] += 1.0;
            invs[nxt].has_thread = true;
            invs[nxt].th = th;
            runnable.push_back(invs[nxt].root);
        } else if (has_setup[task]) {
            start_delayoff(ts, th);
        }
    };
 
    /**
     * One execution of the line's activity has finished.
     *
     * ITS RESPONSE TIME CARRIES THE CALLS IT MADE, which is what separates it
     * from the residence `act_host_resid` accumulates: on an activity that calls
     * a server three times, the response is its own demand plus the three
     * replies and the residence is its own demand alone. The occupancy is closed
     * off through `cur_q` at the same instant, so QLen, RespT and Tput satisfy
     * Little's law by construction rather than by three separate estimators
     * happening to agree.
     */
    auto finish_act = [&](LnLine& L) {
        const std::size_t act = L.act;
        completions[act] += 1.0;
        if (L.act_t0 >= 0.0) {
            resp_sum[act] += now - L.act_t0;
            resp_cnt[act] += 1.0;
            touch(act);
            cur_q[act] -= 1.0;
        }
        L.act_t0 = -1.0;
    };
 
    /**
     * The entry replies: its throughput and response time are recorded, and the
     * synchronous caller is released. `direct` is a reply at completion, whose
     * caller the walk continues itself; an early reply (phase 2 follows, or a
     * cache branch replied) wakes the caller through `runnable`.
     *
     * AN ENTRY'S RESPONSE TIME IS THE LQN ENTRY SERVICE TIME, from the instant it
     * took a thread of its task to the reply, and so carries every nested call its
     * activities made. The wait for the thread belongs to the caller's call, as
     * lqns, lqsim, SolverLN and Solver_ssj_ln.sendEntryReply all count it.
     */
    auto send_reply = [&](std::size_t vi, bool direct) {
        LnInv& v = invs[vi];
        if (v.replied) return;
        v.replied = true;
        v.early_reply = !direct;
        completions[v.entry] += 1.0;
        const double r = now - v.entry_t0;
        resp_sum[v.entry] += r;
        resp_cnt[v.entry] += 1.0;
        // 1-BASED element indices in this engine, unlike the JAR and python
        // ports: entries are eshift+1 .. eshift+nentries, so the local index is
        // one less than the difference.
        if (v.entry > lsn.eshift && v.entry <= lsn.eshift + lsn.nentries)
            entry_resp_samples[v.entry - lsn.eshift - 1].push_back(r);
        if (!direct && v.caller != LN_NONE) runnable.push_back(v.caller);
    };
 
    // ---- cache content -------------------------------------------------------
    /** `cacheMiss`: insert into the entry list, evicting by the policy's rule. */
    auto cache_insert = [&](std::size_t ci, std::size_t item) {
        LnCacheState& cs = caches[ci];
        std::list<std::size_t>& lst = cs.levels[0];
        const std::size_t cap = static_cast<std::size_t>(std::max(0, cs.cap[0]));
        if (cs.rs == lang::ReplacementStrategy::RR) {
            if (lst.size() >= cap && !lst.empty()) {
                std::size_t r = static_cast<std::size_t>(uniform01(g_rr[ci]) *
                                                         static_cast<double>(lst.size()));
                r = std::min(r, lst.size() - 1);
                std::list<std::size_t>::iterator it = lst.begin();
                std::advance(it, static_cast<long>(r));
                *it = item;
            } else {
                lst.push_front(item);
            }
        } else {
            lst.push_front(item);
            if (lst.size() > cap) lst.pop_back();
        }
    };
    /** `cacheHit`: refresh or promote the item one list up, demoting the displaced one. */
    auto cache_promote = [&](std::size_t ci, std::size_t item, std::size_t level) {
        LnCacheState& cs = caches[ci];
        const lang::ReplacementStrategy rs = cs.rs;
        const std::size_t inew = std::min(level + 1, cs.levels.size() - 1);
        const bool front_rule = rs == lang::ReplacementStrategy::LRU ||
                                rs == lang::ReplacementStrategy::HLRU ||
                                rs == lang::ReplacementStrategy::QLRU;
        if (inew <= level) {
            if (front_rule) {
                cs.levels[level].remove(item);
                cs.levels[level].push_front(item);
            }
            return;
        }
        std::list<std::size_t>& lo = cs.levels[level];
        std::list<std::size_t>& hi = cs.levels[inew];
        std::size_t kpos = 0;
        for (std::list<std::size_t>::iterator it = lo.begin(); it != lo.end(); ++it, ++kpos)
            if (*it == item) break;
        lo.remove(item);
        if (hi.size() < static_cast<std::size_t>(std::max(0, cs.cap[inew]))) {
            hi.push_front(item);
            return;
        }
        std::size_t demoted;
        if (rs == lang::ReplacementStrategy::RR) {
            std::size_t r = static_cast<std::size_t>(uniform01(g_rr[ci]) *
                                                     static_cast<double>(hi.size()));
            r = std::min(r, hi.size() - 1);
            std::list<std::size_t>::iterator it = hi.begin();
            std::advance(it, static_cast<long>(r));
            demoted = *it;
            *it = item;
        } else {
            demoted = hi.back();
            hi.pop_back();
            hi.push_front(item);
        }
        if (front_rule || rs == lang::ReplacementStrategy::SFIFO) {
            lo.push_front(demoted);
        } else {
            std::list<std::size_t>::iterator it = lo.begin();
            std::advance(it, static_cast<long>(std::min(kpos, lo.size())));
            lo.insert(it, demoted);
        }
    };
    static const std::size_t HELD = static_cast<std::size_t>(-2);
    /**
     * One read of an ItemEntry, `accessCache`: draw the item, look it up in the
     * live content of the replica serving the request, update the content, and
     * return the branch to continue on -- or HELD, when retrieval parks a read
     * of an item that is already being fetched.
     */
    auto access_cache = [&](std::size_t ai, std::size_t lid) -> std::size_t {
        const LnCacheAccess& ca = accesses[ai];
        LnInv& v = invs[lines[lid].inv];
        const std::size_t ci = ca.caches[std::min(v.trep, ca.caches.size() - 1)];
        LnCacheState& cs = caches[ci];
        ++cs.reads;
        const double u = uniform01(g_item[ci]);
        std::size_t item = ca.cdf.size() - 1;
        for (std::size_t i = 0; i < ca.cdf.size(); ++i)
            if (u <= ca.cdf[i]) {
                item = i;
                break;
            }
        for (std::size_t l = 0; l < cs.levels.size(); ++l)
            if (std::find(cs.levels[l].begin(), cs.levels[l].end(), item) != cs.levels[l].end()) {
                ++cs.hits;
                cache_promote(ci, item, l);
                return ca.hit;
            }
        if (cs.retrieval) {
            if (item < cs.inflight.size() && cs.inflight[item]) {
                cs.held[item].push_back(std::make_pair(lid, ai));
                return HELD;
            }
            ++cs.misses;
            cs.inflight[item] = 1;
            v.fetching = true;
            v.fetch_cache = ci;
            v.fetch_item = item;
            return ca.miss;
        }
        ++cs.misses;
        cache_insert(ci, item);
        return ca.miss;
    };
    /** The fetch completed: insert the item and let every parked read hit it. */
    auto release_fetch = [&](std::size_t vi) {
        LnInv& v = invs[vi];
        if (!v.fetching) return;
        v.fetching = false;
        LnCacheState& cs = caches[v.fetch_cache];
        const std::size_t item = v.fetch_item;
        cs.inflight[item] = 0;
        cache_insert(v.fetch_cache, item);
        std::vector<std::pair<std::size_t, std::size_t> > held;
        held.swap(cs.held[item]);
        for (std::size_t k = 0; k < held.size(); ++k) {
            ++caches[v.fetch_cache].delayed;
            set_act(held[k].first, accesses[held[k].second].hit);
            runnable.push_back(held[k].first);
        }
    };
 
    /**
     * The invocation is complete, `completeRequest`: reply if it has not, leave
     * the entry and the task, free the thread and the task admission token, and
     * return the line to continue -- the synchronous caller, when the reply was
     * made just now -- or LN_NONE. A top-level invocation sends its customer
     * back to think.
     */
    auto complete_request = [&](std::size_t vi) -> std::size_t {
        send_reply(vi, true);
        release_fetch(vi);
        LnInv& v = invs[vi];
        touch(v.entry);
        cur_q[v.entry] -= 1.0;
        completions[v.task] += 1.0;
        touch(v.task);
        cur_q[v.task] -= 1.0;
        release_thread(vi);
        if (v.tcol >= 0) {
            const std::size_t ts = v.ts, task = v.task;
            --tocc[ts][static_cast<std::size_t>(v.tcol)];
            v.tcol = -1;
            while (!tadm[ts].empty()) {
                const std::size_t head = tadm[ts].front();
                const int c = invs[head].tcol;
                if (!admits(task, tocc[ts], c)) break;
                tadm[ts].pop_front();
                ++tocc[ts][static_cast<std::size_t>(c)];
                enter_task(head);
                runnable.push_back(invs[head].root);
            }
        }
        LnInv& w = invs[vi];
        const std::size_t caller = w.caller, cust = w.cust;
        const bool early = w.early_reply;
        free_line(w.root);
        free_invs.push_back(vi);
        if (cust != LN_NONE) {
            // THE COMPLETION IS NOT COUNTED HERE: the reference task's own entry
            // just replied and `completions[task]` above is this task already.
            // `resp_sum` is the CYCLE response time and is the task's own.
            const std::size_t rt = custs[cust].ref_task;
            resp_sum[rt] += now - custs[cust].t_start;
            resp_cnt[rt] += 1.0;
            ++done;
            LnEvent e;
            e.t = now + (has_think[rt] ? think[rt].next(g_think[rt]) : 0.0);
            e.kind = 1;
            e.who = cust;
            push_ev(e);
            return LN_NONE;
        }
        return (caller != LN_NONE && !early) ? caller : LN_NONE;
    };
 
    /** A line reached an activity with no successor, `finishLine`. */
    auto finish_line = [&](std::size_t lid) -> std::size_t {
        const std::size_t vi = lines[lid].inv;
        LnInv& v = invs[vi];
        const bool br = lines[lid].branch;
        if (v.pending < 0) return complete_request(vi);
        const long rem = v.pending - 1;
        if (br) free_line(lid);
        if (rem <= 0) {
            invs[vi].pending = -1;
            return complete_request(vi);
        }
        invs[vi].pending = rem;
        return LN_NONE;
    };
 
    /**
     * The line's activity is complete, `completeActivity`: a cache read picks
     * its branch; an AND-join input is recorded and the join fires on the ROOT
     * once all have arrived; a phase-1 activity entering phase 2, or a cache
     * branch of the entry, replies; then one successor continues, an OR or loop
     * decision draws one, and an AND fork starts every branch. Returns the line
     * to continue, or LN_NONE.
     */
    auto complete_activity = [&](std::size_t lid) -> std::size_t {
        finish_act(lines[lid]);
        const std::size_t act = lines[lid].act;
        const std::size_t vi = lines[lid].inv;
        if (access_of_driver[act] != LN_NONE) {
            const std::size_t nx = access_cache(access_of_driver[act], lid);
            if (nx == HELD) return LN_NONE;
            set_act(lid, nx);
            return lid;
        }
        const std::vector<std::size_t>& sc = succ_act[act];
        const bool br = lines[lid].branch;
        if ((br || invs[vi].pending >= 0) && pretype(act) == lang::PrecedenceType::PRE_AND &&
            !sc.empty()) {
            const std::size_t j = sc[0];
            std::set<std::size_t>& got = invs[vi].joins[j];
            got.insert(act);
            if (br) free_line(lid);
            const std::size_t req = join_required[j] > 0 ? join_required[j] : 1;
            if (got.size() >= req) {
                invs[vi].joins.erase(j);
                if (invs[vi].pending >= 0) invs[vi].pending -= static_cast<long>(req - 1);
                set_act(invs[vi].root, j);
                return invs[vi].root;
            }
            return LN_NONE;
        }
        if (!invs[vi].replied && !br) {
            bool reply = branch_reply[act] != 0 && branch_reply[act] == invs[vi].entry;
            if (!reply && phase_of(act) == 1)
                for (std::size_t s : sc)
                    if (phase_of(s) == 2) {
                        reply = true;
                        break;
                    }
            if (reply) send_reply(vi, false);
        }
        if (invs[vi].replied && invs[vi].fetching) release_fetch(vi);
        if (sc.empty()) return finish_line(lid);
        if (sc.size() == 1) {
            set_act(lid, sc[0]);
            return lid;
        }
        if (!or_prob[act].empty()) {
            const double u = uniform01(g_route[act]);
            double cum = 0.0;
            std::size_t pick = sc.back();
            for (std::size_t k = 0; k < sc.size(); ++k) {
                cum += or_prob[act][k];
                if (u <= cum) {
                    pick = sc[k];
                    break;
                }
            }
            set_act(lid, pick);
            return lid;
        }
        // AND fork: one line becomes |sc| concurrent lines. The forking line goes
        // dormant (the root) or retires (a branch); the join resumes the root.
        invs[vi].pending = (invs[vi].pending < 0 ? 1 : invs[vi].pending) +
                           static_cast<long>(sc.size()) - 1;
        if (br) free_line(lid);
        std::vector<std::size_t> kids;
        for (std::size_t k = 0; k < sc.size(); ++k) kids.push_back(new_line(vi, sc[k], true));
        for (std::size_t k = kids.size(); k-- > 1;) runnable.push_front(kids[k]);
        return kids[0];
    };
 
    /**
     * A request at `entry` on task replica `rep` that nobody waits on: an asynchronous callee or
     * an open arrival. It runs as an invocation with no caller and no customer, so it just ends.
     */
    auto spawn_request = [&](std::size_t entry, std::size_t rep) {
        const std::size_t vi = new_inv(entry, rep, LN_NONE, LN_NONE);
        if (arrive_at_task(vi)) runnable.push_back(invs[vi].root);
    };
 
    /**
     * Continue line `lid` from wherever it is: take a thread, begin the current
     * activity, ask its host, think, issue its next call, or complete it.
     */
    advance = [&](std::size_t lid) {
        while (lid != LN_NONE) {
            const std::size_t vi = lines[lid].inv;
            // An invocation holds a THREAD of its task for the whole entry,
            // across its calls, and waits outside when none is free.
            if (!invs[vi].has_thread) {
                if (!acquire_thread(vi)) return;
            }
            // The entry's service starts with its thread; its occupancy runs from here to the
            // end of its last phase, Solver_ssj_ln's entryResidenceTime from serviceStartTime.
            if (!invs[vi].started) {
                LnInv& v = invs[vi];
                v.started = true;
                v.entry_t0 = now;
                touch(v.entry);
                cur_q[v.entry] += 1.0;
            }
            LnLine& L = lines[lid];
            const std::size_t act = L.act;
            if (act == 0) {  // an entry with no activity replies at once
                lid = complete_request(vi);
                continue;
            }
            if (L.act_t0 < 0.0) {
                L.act_t0 = now;
                touch(act);
                cur_q[act] += 1.0;
                L.dem_done = 0;
            }
            // THE HOST DEMAND RUNS FIRST, THEN THE THINK TIME, THEN THE CALLS: the order LQN2QN
            // builds (the demand step, then one step per call stage) and Solver_ssj_ln runs.
            if (L.dem_done == 0) {
                L.dem_done = 1;
                if (has_hostdem[act]) {
                    request_host(lid);
                    return;
                }
            }
            // The activity's think time, Solver_ssj_ln.proceedAfterHost: it holds the thread but
            // not the processor, so it adds to the activity's response and not to host busy.
            if (L.dem_done == 1) {
                L.dem_done = 2;
                if (has_actthink[act]) {
                    LnEvent e;
                    e.t = now + actthink[act].next(g_actthink[act]);
                    e.kind = 6;
                    e.who = lid;
                    push_ev(e);
                    return;
                }
            }
            const std::vector<std::size_t>& calls =
                (act < lsn.callsof.size()) ? lsn.callsof[act] : std::vector<std::size_t>();
            if (L.call_pos < calls.size()) {
                // A pending call of the current activity, made `calls-mean` times
                // per execution, drawn once when it is first reached.
                const std::size_t cidx = calls[L.call_pos];
                if (L.calls_left == NOTDRAWN) L.calls_left = sample_call_count(cidx);
                if (L.calls_left == 0) {
                    ++L.call_pos;
                    L.calls_left = NOTDRAWN;
                    continue;
                }
                --L.calls_left;
                if (L.calls_left == 0) {
                    ++L.call_pos;
                    L.calls_left = NOTDRAWN;
                }
                const std::size_t callee = lsn.callpair_dst[cidx];
                const std::size_t rep =
                    callee_replica(invs[vi].task, invs[vi].trep, lsn.parent[callee]);
                if (lsn.calltype[cidx] == lang::CallType::SYNC) {
                    // THE CALLER KEEPS ITS THREAD across the call. That is the
                    // layered contention the model exists to represent.
                    const std::size_t child = new_inv(callee, rep, lid, LN_NONE);
                    if (!arrive_at_task(child)) return;
                    lid = invs[child].root;
                    continue;
                }
                // Asynchronous: the callee runs as a request of its own, the caller
                // does not wait, and the reply is never collected.
                spawn_request(callee, rep);
                continue;
            }
            // The demand, the think time and the calls are done: the activity
            // completes, its response time carrying all three.
            L.call_pos = 0;
            L.calls_left = NOTDRAWN;
            lid = complete_activity(lid);
        }
    };
 
    // ---- reference tasks ----------------------------------------------------
    // Each reference task runs `mult` independent customers, which is what makes
    // the closed population of a layered model.
    auto start_cycle = [&](std::size_t c) {
        const std::size_t t = custs[c].ref_task;
        custs[c].t_start = now;
        const std::size_t vi = new_inv(lsn.entriesof[t][0], custs[c].repl, LN_NONE, c);
        if (arrive_at_task(vi)) advance(invs[vi].root);
    };
    bool any_ref = false;
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        if (t >= lsn.isref.size() || !lsn.isref[t]) continue;
        any_ref = true;
        const std::size_t n = (task_threads[t] == std::numeric_limits<std::size_t>::max())
                                  ? 1
                                  : task_threads[t];
        const std::vector<std::size_t>& es =
            (t < lsn.entriesof.size()) ? lsn.entriesof[t] : std::vector<std::size_t>();
        if (es.empty())
            throw InputError("SolverLDES (native LN engine): reference task '" + lsn.names[t] +
                             "' declares no entry to invoke");
        // A replicated reference task is r separate populations of `mult`
        // customers, not one population of r*mult sharing a queue.
        for (std::size_t m = 0; m < repl[t]; ++m) {
            for (std::size_t k = 0; k < n; ++k) {
                LnCustomer cu;
                cu.ref_task = t;
                cu.repl = m;
                custs.push_back(cu);
                start_cycle(custs.size() - 1);
            }
        }
    }
    // ---- open arrivals --------------------------------------------------------
    // Each replica of the entry's task receives its own stream, as Solver_ssj_ln.startOpenArrivals
    // gives it: the declared rate is what one replica sees. Offset block 5000, per entry.
    std::vector<Sampler> arrival(nidx + 1);
    std::vector<Rng> g_arr;
    g_arr.reserve(nidx + 1);
    for (std::size_t k = 0; k <= nidx; ++k)
        g_arr.push_back(Rng(ln_seed, static_cast<long long>(k) * 10 + 5000));
    for (std::size_t e = lsn.eshift + 1; e <= lsn.eshift + lsn.nentries; ++e) {
        if (e >= lsn.has_arrival.size() || !lsn.has_arrival[e] || e >= lsn.arrival.size() ||
            lsn.arrival[e].disabled)
            continue;
        arrival[e] = Sampler(lsn.arrival[e], "the open arrivals at '" + lsn.names[e] + "'");
        for (std::size_t m = 0; m < repl[lsn.parent[e]]; ++m) {
            LnEvent ev;
            ev.t = now + arrival[e].next(g_arr[e]);
            ev.kind = 7;
            ev.who = e;
            ev.aux = m;
            push_ev(ev);
        }
    }
    drain();
    if (!any_ref)
        throw InputError("SolverLDES (native LN engine): the model has no reference task, so "
                         "nothing drives it");
 
    // ---- the event loop -----------------------------------------------------
    while (!evq.empty() && done < max_events) {
        const LnEvent ev = evq.top();
        evq.pop();
        now = ev.t;
 
        if (ev.kind == 1) {
            // A reference customer finished thinking: start the next cycle on the
            // replica it belongs to.
            start_cycle(ev.who);
            drain();
            continue;
        }
 
        if (ev.kind == 2 || ev.kind == 3) {
            const std::size_t ts = ev.who, th = ev.aux;
            if (ev.gen != thr_gen[ts][th]) continue;  // the timer was cancelled
            if (ev.kind == 3) {
                if (thr_state[ts][th] == ThreadState::DELAYOFF)
                    thr_state[ts][th] = ThreadState::OFF;
                continue;
            }
            // A cold start finished: the thread is up and takes the head of the
            // queue it was woken for.
            thr_state[ts][th] = ThreadState::ACTIVE;
            if (!thr_held[ts][th] && !thr_queue[ts].empty()) {
                const std::size_t nxt = dequeue(thr_queue[ts], task_sched[thr_owner[ts]], inv_key);
                thr_held[ts][th] = true;
                touch(thr_owner[ts]);
                cur_u[thr_owner[ts]] += 1.0;
                invs[nxt].has_thread = true;
                invs[nxt].th = th;
                advance(invs[nxt].root);
                drain();
            }
            continue;
        }
 
        if (ev.kind == 4) {
            // The next completion at a pooled replica, unless the job set moved it.
            const std::size_t hs = ev.who;
            if (ev.gen != pool_gen[hs]) continue;
            const std::size_t host = host_of_slot[hs];
            advance_pool(hs);
            std::vector<std::size_t> finished, keep;
            for (std::size_t lid : pool_jobs[hs])
                (lines[lid].remaining <= 1e-9 ? finished : keep).push_back(lid);
            pool_jobs[hs].swap(keep);
            rearm_pool(hs);
            for (std::size_t lid : finished) {
                touch(host);
                cur_q[host] -= 1.0;
                completions[host] += 1.0;
                act_host_resid[lines[lid].act] += now - lines[lid].host_t0;
                release_host_admission(lid);
                advance(lid);  // on to its think time and calls
                drain();
            }
            continue;
        }
 
        if (ev.kind == 5) {
            // A PS host replica's next completion: every line whose work is done leaves together.
            const std::size_t hs = ev.who;
            if (ev.gen != ps_gen[hs]) continue;  // the job set changed since this was armed
            const std::size_t host = host_of_slot[hs];
            ps_advance(hs);
            std::vector<std::size_t> finished, keep;
            for (std::size_t lid : ps_jobs[hs])
                (lines[lid].remaining <= 1e-9 ? finished : keep).push_back(lid);
            const std::size_t n0 = ps_jobs[hs].size();
            ps_jobs[hs].swap(keep);
            touch(host);
            cur_u[host] -= ps_busy(hs, n0) - ps_busy(hs, ps_jobs[hs].size());
            ps_schedule(hs);
            for (std::size_t lid : finished) {
                cur_q[host] -= 1.0;
                --host_busy[hs];
                completions[host] += 1.0;
                act_host_resid[lines[lid].act] += now - lines[lid].host_t0;
                release_host_admission(lid);
                advance(lid);  // on to its think time and calls
                drain();
            }
            continue;
        }
 
        if (ev.kind == 6) {
            // A line finished its activity's think time and goes on to its calls.
            advance(ev.who);
            drain();
            continue;
        }
 
        if (ev.kind == 7) {
            // An open arrival: a request nobody waits on, then the next arrival of the stream.
            spawn_request(ev.who, ev.aux);
            LnEvent nx = ev;
            nx.t = now + arrival[ev.who].next(g_arr[ev.who]);
            push_ev(nx);
            drain();
            continue;
        }
 
        // A host completion at a c-server (or INF) processor, unless a preemption made it stale.
        const std::size_t lid = ev.who;
        if (ev.gen != lines[lid].host_gen) continue;
        const std::size_t act = lines[lid].act;
        const std::size_t host = lsn.host_of(act);
        // The demand ran on the processor replica paired with the task replica
        // that holds the request, and that pairing is fixed for the invocation.
        const std::size_t hs = host_slot(host, invs[lines[lid].inv].trep);
        touch(host);
        cur_q[host] -= 1.0;
        cur_u[host] -= 1.0;
        --host_busy[hs];
        completions[host] += 1.0;
        act_host_resid[act] += now - lines[lid].host_t0;
        if (preemptive(host_sched[host]))
            in_service[hs].erase(std::find(in_service[hs].begin(), in_service[hs].end(), lid));
        if (!host_queue[hs].empty()) {
            const std::size_t nxt = dequeue(host_queue[hs], host_sched[host], line_key);
            ++host_busy[hs];
            cur_u[host] += 1.0;
            begin_host(nxt, hs);
        }
        release_host_admission(lid);
        // The demand is the FIRST thing an activity does: it goes on to its think
        // time and its calls, and advance() completes it once they are done.
        advance(lid);
        drain();
    }
 
    // A run that stops early with a request still blocked by an admission
    // constraint did not finish its budget: the constraint deadlocked it, and
    // the measures so far describe a model that has stopped moving.
    if (done < max_events && evq.empty()) {
        for (std::size_t s = 0; s < nhost_slots; ++s)
            if (!hadm[s].empty())
                throw InputError("SolverLDES (native LN engine): the admission constraint on "
                                 "host '" + lsn.names[host_of_slot[s]] + "' deadlocked the model: " +
                                 std::to_string(hadm[s].size()) +
                                 " request(s) are blocked with no event left");
        for (std::size_t s = 0; s < ntask_slots; ++s)
            if (!tadm[s].empty())
                throw InputError("SolverLDES (native LN engine): the admission constraint on "
                                 "task '" + lsn.names[thr_owner[s]] + "' deadlocked the model: " +
                                 std::to_string(tadm[s].size()) +
                                 " request(s) are blocked with no event left");
    }
 
    // ---- result -------------------------------------------------------------
    for (std::size_t k = 1; k <= nidx; ++k) touch(k);
    LnResult res;
    res.nidx = nidx;
    res.QLN = Matrix<double>(nidx + 1, 1, 0.0);
    res.ULN = Matrix<double>(nidx + 1, 1, 0.0);
    res.RLN = Matrix<double>(nidx + 1, 1, 0.0);
    res.TLN = Matrix<double>(nidx + 1, 1, 0.0);
    // Residence starts at NaN everywhere and is written only where it exists:
    // activities and their tasks. See LnResult::WLN.
    res.WLN = Matrix<double>(nidx + 1, 1, std::numeric_limits<double>::quiet_NaN());
    for (std::size_t k = 1; k <= nidx; ++k) {
        if (now > 0.0) {
            res.QLN(k, 0) = tot_q[k] / now;
            // A host's utilization is its busy servers over its capacity; an
            // infinite server has no capacity to be busy against, so the
            // integral itself is the answer.
            // Replicas add capacity: r copies of a c-server processor hold r*c
            // servers, and the integral above already runs over all of them.
            // A task's capacity is its threads, the same way a processor's is its
            // servers: both are now integrated above, so both normalise the same.
            const bool is_task = (k > lsn.tshift && k <= lsn.tshift + lsn.ntasks);
            const std::size_t cap_k = is_task ? task_threads[k] : host_servers[k];
            const double c = (cap_k == std::numeric_limits<std::size_t>::max())
                                 ? 1.0
                                 : static_cast<double>(cap_k * repl[k]);
            res.ULN(k, 0) = tot_u[k] / (now * c);
            res.TLN(k, 0) = completions[k] / now;
        }
        if (resp_cnt[k] > 0.0) res.RLN(k, 0) = resp_sum[k] / resp_cnt[k];
    }
    res.entry_resp_samples.swap(entry_resp_samples);
 
    // ResidT, per activity and summed onto its task. SolverLN builds
    // residt(aidx) as QN(host)/TN_ref, and TN_ref -- the reference of the
    // activity's layer -- is the completion rate of the activity's TASK: the
    // cycle rate of a reference task, the request rate into a called one,
    // TLN(task) either way. QN(host) is the mean occupancy at the host, which
    // is the total host time over the run divided by it, so
    //
    //   ResidT = QN(host) / X_task = (total host time) / (elapsed * X_task)
    //
    // i.e. host time accumulated per driving visit. An activity executing once
    // per request reports its host residence, one inside a loop the loop's
    // total, and neither reading needs a visit count of its own. Accumulated
    // aside so a task with no activity at all keeps its NaN instead of 0.
    if (now > 0.0) {
        std::vector<double> task_resid(nidx + 1, 0.0);
        std::vector<bool> task_has_act(nidx + 1, false);
        for (std::size_t a = lsn.ashift + 1; a <= lsn.ashift + lsn.nacts && a <= nidx; ++a) {
            const std::size_t task = lsn.parent[a];
            const double x_task = (task >= 1 && task <= nidx) ? res.TLN(task, 0) : 0.0;
            const double resid = (x_task > 0.0) ? act_host_resid[a] / (now * x_task) : 0.0;
            res.WLN(a, 0) = resid;
            if (task >= 1 && task <= nidx) {
                task_resid[task] += resid;
                task_has_act[task] = true;
            }
        }
        for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks && t <= nidx; ++t)
            if (task_has_act[t]) res.WLN(t, 0) = task_resid[t];
    }
 
    // ENTRY AND ACTIVITY UTILIZATION ARE DERIVED, the way `Solver_ssj_ln`
    // derives them: the utilization an activity causes at its host is its
    // throughput times its host demand, an entry's is the sum over the
    // activities bound under it, and an activity's own figure is that share of
    // its task's total. There is nothing to integrate for either kind -- an
    // entry is not a server and holds none -- so the `cur_u` integral that
    // answers for hosts and tasks leaves both at zero, which read as idle.
    //
    // The HOST AND TASK columns are NOT touched here: this port measures them as
    // the integral of busy servers and of busy threads, which is a different
    // quantity from the JAR's derived task utilization (sum of X*D over the
    // task's activities) and the one its own tests pin. The two ports therefore
    // still disagree on the task Util column; see `_kb/09-ldes-and-cache.md`.
    if (now > 0.0) {
        std::vector<double> proc_util(nidx + 1, 0.0), task_proc_util(nidx + 1, 0.0);
        for (std::size_t a = lsn.ashift + 1; a <= lsn.ashift + lsn.nacts && a <= nidx; ++a) {
            const double d = (a < lsn.hostdem.size() && !lsn.hostdem[a].disabled)
                                 ? num_traits<T>::to_double(lsn.hostdem[a].mean)
                                 : 0.0;
            proc_util[a] = res.TLN(a, 0) * d;
            const std::size_t task = lsn.parent[a];
            if (task >= 1 && task <= nidx) task_proc_util[task] += proc_util[a];
        }
        for (std::size_t e = lsn.eshift + 1; e <= lsn.eshift + lsn.nentries && e <= nidx; ++e) {
            double u = 0.0;
            const std::vector<std::size_t>& acts =
                (e < lsn.actsof.size()) ? lsn.actsof[e] : std::vector<std::size_t>();
            for (std::size_t i = 0; i < acts.size(); ++i)
                if (acts[i] <= nidx) u += proc_util[acts[i]];
            res.ULN(e, 0) = u;
        }
        for (std::size_t a = lsn.ashift + 1; a <= lsn.ashift + lsn.nacts && a <= nidx; ++a) {
            const std::size_t task = lsn.parent[a];
            const double tu = (task >= 1 && task <= nidx) ? task_proc_util[task] : 0.0;
            if (tu > 0.0 && tu < 1.0)
                res.ULN(a, 0) = std::min(1.0, proc_util[a] / tu);
            else
                res.ULN(a, 0) = (proc_util[a] > 0.0) ? 1.0 : 0.0;
        }
    }
 
    // Cache measures per ItemEntry, pooled over the replicas of its cache task.
    for (std::size_t ai = 0; ai < accesses.size(); ++ai) {
        long long reads = 0, hits = 0, misses = 0, delayed = 0;
        for (std::size_t ci : accesses[ai].caches) {
            reads += caches[ci].reads;
            hits += caches[ci].hits;
            misses += caches[ci].misses;
            delayed += caches[ci].delayed;
        }
        if (reads <= 0 || !(now > 0.0)) continue;
        const double n = static_cast<double>(reads);
        res.cache_task.push_back(lsn.parent[accesses[ai].entry]);
        res.cache_entry.push_back(accesses[ai].entry);
        res.cache_hit.push_back(static_cast<double>(hits) / n);
        res.cache_miss.push_back(static_cast<double>(misses) / n);
        res.cache_delayed.push_back(static_cast<double>(delayed) / n);
        res.cache_read_rate.push_back(n / now);
    }
 
    res.simulated_time = now;
    res.completions = static_cast<long long>(done);
    return res;
}
 
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_LN_ENGINE_H