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 THIS INCREMENT COVERS: reference tasks with think time and
 * multiplicity, entries with sequential activity graphs, host demand under
 * FCFS / PS / INF host scheduling, task multiplicity as a thread semaphore,
 * synchronous calls (blocking, with a mean call multiplicity) and asynchronous
 * calls (fire and forget), open arrivals at an entry, and REPLICATION of a
 * processor or a task, materialised one queue per copy.
 *
 * WHAT IT REFUSES by name: activity graphs with AND/OR forks or loops, cache
 * tasks, setup/delay-off on a task, and admission constraints.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <functional>
#include <limits>
#include <map>
#include <queue>
#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, 1 = think completion, 2 = a task thread finished its
     * cold start, 3 = an idle task thread's delay-off countdown expired.
     *
     * 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. 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;
    }
};
 
/** One in-flight request: a reference-task job somewhere in its call tree. */
struct LnJob {
    std::size_t id = 0;
    std::size_t ref_task = 0;
    /** The activity stack: what this job is executing and what it will return to. */
    std::vector<std::size_t> stack;      ///< entry indices, innermost last
    std::vector<std::size_t> act_pos;    ///< next activity of each stacked entry
    /**
     * Repetitions of the CURRENT call still to be made, per stack level.
     *
     * `NOTDRAWN` means the count has not been sampled yet for this execution of
     * the activity. An LQN call carries a MEAN NUMBER OF CALLS (`calls-mean`),
     * not a single dispatch, and the count is drawn once when the call is first
     * reached -- so `A.synchCall(E, 3)` blocks on E three times per execution of
     * A, which is three times the callee's demand in the caller's cycle.
     */
    std::vector<std::size_t> calls_left;
    std::vector<std::size_t> call_pos;   ///< next call index of the current activity
    /** When each stacked entry was entered, for its response time and occupancy. */
    std::vector<double> entry_t0;
    /**
     * When the current activity of each level began; -1 between activities.
     *
     * An activity's span runs from the moment it becomes current to the moment
     * its host demand completes, so it CARRIES THE NESTED CALLS it makes -- that
     * is what makes it a response time and not a residence.
     */
    std::vector<double> act_t0;
    /**
     * Replica of the task owning each stacked entry. A call picks the callee
     * replica once and the request stays on it until it replies, so the whole
     * of one invocation runs on one copy of the callee and its processor.
     */
    std::vector<std::size_t> repl_of;
    /**
     * The task thread each stacked entry holds. A caller keeps its thread across
     * a synchronous call, so every level of the stack holds one of its own task's
     * threads at once, and that is the layered contention the model represents.
     * `NOTHR` marks a level whose thread is not yet held (it is queueing, or its
     * thread is powering up) or a level on an infinite-thread task.
     */
    std::vector<std::size_t> thr_slot;
    std::vector<std::size_t> thr_id;
    double t_start = 0.0;                ///< when the current top-level invocation began
    /**
     * When this job last asked its host processor for service -- the instant
     * the residence at the host starts, queueing included. A job is at one host
     * at a time, so one scalar carries it.
     */
    double host_t0 = 0.0;
    bool holding_thread = false;
    std::size_t thread_task = 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;
};
 
/** 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 this increment does not simulate -----------------------
    for (std::size_t t = lsn.tshift + 1; t <= lsn.tshift + lsn.ntasks; ++t) {
        if (t < lsn.iscache.size() && lsn.iscache[t])
            throw UnsupportedError("SolverLDES (native LN engine): cache task '" + lsn.names[t] +
                                   "' is not ported yet");
        // A setup IS simulated now, 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");
        }
    }
    // The Java engine (Solver_ssj_ln) admits requests against A n <= b and parks
    // the ones that would breach it; this engine does not, so it must say so
    // rather than return the unconstrained sample path under a constrained model.
    for (std::size_t k = 1; k < lsn.lincon_A.size(); ++k) {
        if (lsn.lincon_A[k].rows() > 0)
            throw UnsupportedError("SolverLDES (native LN engine): the admission constraint on '" +
                                   lsn.names[k] + "' is not ported yet; solve this model with the "
                                   "Java LDES engine");
    }
    // Heterogeneous server pools on a layer server. The Java engine holds them
    // concretely -- a request occupies one server of one compatible pool, and
    // under PS the per-job rates are the max-min fair allocation of
    // SnCompatShare -- so a model that declares them is a DIFFERENT system from
    // the multiserver it would otherwise look like. Flattening them into the
    // multiplicity here would finish and report plausible numbers for that other
    // system, which is the failure this validation exists to prevent.
    for (std::size_t k = 0; k < lsn.pools.size(); ++k) {
        if (lsn.pools[k].npools() > 0)
            throw UnsupportedError("SolverLDES (native LN engine): '" + lsn.names[k] +
                                   "' declares heterogeneous server pools, which are not ported "
                                   "yet; solve this model with the Java LDES engine");
    }
    for (std::size_t a = lsn.ashift + 1; a <= lsn.ashift + lsn.nacts; ++a) {
        if (a < lsn.actpretype.size() && lsn.actpretype[a] != lang::PrecedenceType::NONE &&
            lsn.actpretype[a] != lang::PrecedenceType::PRE_SEQ)
            throw UnsupportedError("SolverLDES (native LN engine): activity '" + lsn.names[a] +
                                   "' has a non-sequential precedence, which is not ported yet");
        if (a < lsn.actposttype.size() && lsn.actposttype[a] != lang::PrecedenceType::NONE &&
            lsn.actposttype[a] != lang::PrecedenceType::POST_SEQ)
            throw UnsupportedError("SolverLDES (native LN engine): activity '" + lsn.names[a] +
                                   "' has a non-sequential post-precedence, which is not ported "
                                   "yet");
    }
 
    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);
    std::vector<bool> has_hostdem(nidx + 1, false), has_think(nidx + 1, false);
    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);
    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;
        host_servers[h] = std::isfinite(m) ? static_cast<std::size_t>(m + 0.5)
                                           : std::numeric_limits<std::size_t>::max();
        if (h < lsn.sched.size()) host_sched[h] = lsn.sched[h];
    }
    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;
        task_threads[t] = std::isfinite(m) ? 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;
    };
 
    // ---- live state ---------------------------------------------------------
    double now = 0.0;
    std::priority_queue<LnEvent, std::vector<LnEvent>, LnEventLater> evq;
    std::uint64_t seq = 0;
    std::map<std::size_t, LnJob> jobs;
    std::size_t next_job = 0;
 
    // Busy servers and waiting jobs 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);
    std::vector<std::size_t> running_on(nhost_slots, 0);  ///< job currently on a host slot
 
    // Task threads of ONE task replica, addressed by slot. A request 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);
        }
    }
 
    // 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;
    }
 
    /// Jobs woken by a thread release, continued after the current walk returns.
    /// Waking them inside the walk would re-enter it on another job while this
    /// one still holds a reference into `jobs`.
    std::deque<std::size_t> runnable;
 
    // 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;
 
    // ---- the activity walk --------------------------------------------------
    std::function<void(std::size_t)> advance;
    /** Continue every job a thread release woke, and every job those wake. */
    auto drain = [&]() {
        while (!runnable.empty()) {
            const std::size_t id = runnable.front();
            runnable.pop_front();
            advance(id);
        }
    };
 
    /** Put `job` on its host replica, queueing when every server of it is taken. */
    auto start_host = [&](std::size_t jid, std::size_t act, std::size_t trep) {
        const std::size_t task = lsn.parent[act];
        const std::size_t host = lsn.host_of(act);
        const std::size_t hs = host_slot(host, trep);
        touch(host);
        cur_q[host] += 1.0;
        // Stamped BEFORE the branch: a request that has to queue starts its
        // residence when it asks, not when a server frees up.
        jobs[jid].host_t0 = now;
        if (host_busy[hs] < host_servers[host]) {
            ++host_busy[hs];
            cur_u[host] += 1.0;
            running_on[hs] = jid;
            LnEvent e;
            e.t = now + hostdem[act].next(g_host[act]);
            e.kind = 0;
            e.who = jid;
            push_ev(e);
        } else {
            host_queue[hs].push_back(jid);
        }
        (void)task;
    };
 
    /** 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 `jid` a thread of the task owning the entry on top of its stack.
     *
     * Returns false when the job could not proceed: either every thread is taken,
     * or the only thread available was powered off and is now warming up. Either
     * way the job is parked in thr_queue and will be 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 jid, std::size_t entry) -> bool {
        LnJob& j = jobs[jid];
        const std::size_t task = lsn.parent[entry];
        const std::size_t trep = j.repl_of.back();
        const std::size_t ts = task_slot(task, trep);
        touch(task);
        if (task_threads[task] == std::numeric_limits<std::size_t>::max()) {
            cur_u[task] += 1.0;
            j.thr_slot.back() = ts;
            j.thr_id.back() = 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;
                j.thr_slot.back() = ts;
                j.thr_id.back() = 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;
                j.thr_slot.back() = ts;
                j.thr_id.back() = th;
                return true;
            }
        }
        thr_queue[ts].push_back(jid);
        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 the innermost entry of `jid` holds and wake the next job. */
    auto release_thread = [&](std::size_t jid) {
        LnJob& j = jobs[jid];
        const std::size_t ts = j.thr_slot.back();
        const std::size_t th = j.thr_id.back();
        const std::size_t task = thr_owner[ts];
        touch(task);
        cur_u[task] -= 1.0;
        j.thr_slot.back() = NOTHR;
        j.thr_id.back() = 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 = thr_queue[ts].front();
            thr_queue[ts].pop_front();
            LnJob& nj = jobs[nxt];
            thr_held[ts][th] = true;
            cur_u[task] += 1.0;
            nj.thr_slot.back() = ts;
            nj.thr_id.back() = th;
            runnable.push_back(nxt);
        } else if (has_setup[task]) {
            start_delayoff(ts, th);
        }
    };
 
    /**
     * Continue `jid` from wherever it is: run the next activity of the innermost
     * entry, issue its next call, or return to the caller.
     */
    /**
     * One execution of `act` 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 = [&](LnJob& j, std::size_t act) {
        const double t0 = j.act_t0.empty() ? -1.0 : j.act_t0.back();
        completions[act] += 1.0;
        if (t0 >= 0.0) {
            resp_sum[act] += now - t0;
            resp_cnt[act] += 1.0;
            touch(act);
            cur_q[act] -= 1.0;
        }
        if (!j.act_t0.empty()) j.act_t0.back() = -1.0;
    };
 
    advance = [&](std::size_t jid) {
        LnJob& j = jobs[jid];
        while (true) {
            if (j.stack.empty()) {
                // The invocation is finished: the reference task thinks again.
                //
                // THE COMPLETION IS NOT COUNTED HERE. The reference task's own
                // entry just replied, and the reply branch below already charged
                // `completions[parent(entry)]`, which for a reference task IS
                // this task -- so counting it again reported TWICE the cycle
                // rate for every reference task, and with it half the residence
                // of every activity on one (ResidT divides by TLN(task)).
                // `resp_sum` is the CYCLE response time and is the task's own,
                // so it does belong here.
                const std::size_t rt = j.ref_task;
                resp_sum[rt] += now - j.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 = jid;
                push_ev(e);
                return;
            }
            const std::size_t entry = j.stack.back();
            // A request holds a THREAD of the entry's task for the whole entry,
            // across its calls, and waits outside when none is free.
            if (j.thr_slot.back() == NOTHR && j.thr_id.back() == NOTHR) {
                if (!acquire_thread(jid, entry)) return;
            }
            const std::vector<std::size_t>& acts =
                (entry < lsn.actsof.size()) ? lsn.actsof[entry] : std::vector<std::size_t>();
            if (j.act_pos.back() >= acts.size()) {
                // Every activity of this entry has run: the entry replies.
                //
                // AN ENTRY'S RESPONSE TIME IS THE WHOLE INVOCATION, from the
                // instant the request reached it -- queueing for a thread of its
                // task included -- to this reply, and so carries every nested
                // call its activities made. Its occupancy follows from the same
                // interval through `cur_q`, which is Little's law rather than a
                // second measurement of the same thing.
                completions[entry] += 1.0;
                const double e0 = j.entry_t0.back();
                resp_sum[entry] += now - e0;
                resp_cnt[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 (entry > lsn.eshift && entry <= lsn.eshift + lsn.nentries) {
                    entry_resp_samples[entry - lsn.eshift - 1].push_back(now - e0);
                }
                touch(entry);
                cur_q[entry] -= 1.0;
                const std::size_t task = lsn.parent[entry];
                completions[task] += 1.0;
                touch(task);
                cur_q[task] -= 1.0;
                release_thread(jid);
                j.stack.pop_back();
                j.act_pos.pop_back();
                j.call_pos.pop_back();
                j.calls_left.pop_back();
                j.entry_t0.pop_back();
                j.act_t0.pop_back();
                j.repl_of.pop_back();
                j.thr_slot.pop_back();
                j.thr_id.pop_back();
                continue;
            }
            const std::size_t act = acts[j.act_pos.back()];
            // The activity begins the first time it is reached and stays open
            // across its calls, so `act_t0` is set once and cleared when it
            // completes. -1 marks the gap between two activities of one entry.
            if (j.act_t0.back() < 0.0) {
                j.act_t0.back() = now;
                touch(act);
                cur_q[act] += 1.0;
            }
            const std::vector<std::size_t>& calls =
                (act < lsn.callsof.size()) ? lsn.callsof[act] : std::vector<std::size_t>();
            if (j.call_pos.back() < calls.size()) {
                // A pending call of the current activity. HOW MANY TIMES it is
                // made is `calls-mean`, drawn once per execution: walking the
                // call LIST and issuing each entry once ignored the multiplicity
                // entirely, so `synchCall(E, 3)` blocked on E once and every
                // number on such a model was computed on the wrong cycle.
                const std::size_t cidx = calls[j.call_pos.back()];
                if (j.calls_left.back() == NOTDRAWN)
                    j.calls_left.back() = sample_call_count(cidx);
                if (j.calls_left.back() == 0) {
                    // A mean below 1 can draw zero: this call is not made at all
                    // on this execution.
                    ++j.call_pos.back();
                    j.calls_left.back() = NOTDRAWN;
                    continue;
                }
                --j.calls_left.back();
                if (j.calls_left.back() == 0) {
                    ++j.call_pos.back();
                    j.calls_left.back() = NOTDRAWN;
                }
                const std::size_t callee = lsn.callpair_dst[cidx];
                const bool sync = (lsn.calltype[cidx] == lang::CallType::SYNC);
                if (sync) {
                    // THE CALLER KEEPS ITS THREAD across the call. That is the
                    // layered contention the model exists to represent.
                    j.stack.push_back(callee);
                    j.act_pos.push_back(0);
                    j.call_pos.push_back(0);
                    j.calls_left.push_back(NOTDRAWN);
                    j.entry_t0.push_back(now);
                    j.act_t0.push_back(-1.0);
                    j.repl_of.push_back(callee_replica(lsn.parent[entry], j.repl_of.back(),
                                                       lsn.parent[callee]));
                    j.thr_slot.push_back(NOTHR);
                    j.thr_id.push_back(NOTHR);
                    touch(lsn.parent[callee]);
                    cur_q[lsn.parent[callee]] += 1.0;
                    touch(callee);
                    cur_q[callee] += 1.0;
                    continue;
                }
                // Asynchronous: the callee runs on its own, the caller does not
                // wait, and the reply is never collected.
                completions[callee] += 1.0;
                continue;
            }
            // The activity's calls are done; run its host demand, if any.
            ++j.act_pos.back();
            j.call_pos.back() = 0;
            j.calls_left.back() = NOTDRAWN;
            if (has_hostdem[act]) {
                start_host(jid, act, j.repl_of.back());
                return;
            }
            // No host demand: the activity ends here rather than at a host
            // completion, and still has a response time -- the calls it made.
            finish_act(j, act);
        }
    };
 
    // ---- reference tasks ----------------------------------------------------
    // Each reference task runs `mult` independent jobs, which is what makes the
    // closed population of a layered model.
    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) {
                LnJob j;
                j.id = ++next_job;
                j.ref_task = t;
                j.t_start = 0.0;
                j.stack.push_back(es[0]);
                j.act_pos.push_back(0);
                j.call_pos.push_back(0);
                j.calls_left.push_back(NOTDRAWN);
                j.entry_t0.push_back(now);
                j.act_t0.push_back(-1.0);
                j.repl_of.push_back(m);
                j.thr_slot.push_back(NOTHR);
                j.thr_id.push_back(NOTHR);
                jobs[j.id] = j;
                touch(t);
                cur_q[t] += 1.0;
                touch(es[0]);
                cur_q[es[0]] += 1.0;
                advance(j.id);
            }
        }
    }
    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 task finished thinking: start the next invocation.
            LnJob& j = jobs[ev.who];
            const std::vector<std::size_t>& es = lsn.entriesof[j.ref_task];
            j.t_start = now;
            const std::size_t rep = j.repl_of.empty() ? 0 : j.repl_of.front();
            j.stack.assign(1, es[0]);
            j.act_pos.assign(1, 0);
            j.call_pos.assign(1, 0);
            j.calls_left.assign(1, NOTDRAWN);
            j.entry_t0.assign(1, now);
            j.act_t0.assign(1, -1.0);
            // The customer belongs to one replica of its reference task for the
            // whole run: it returns to the copy it came from.
            j.repl_of.assign(1, rep);
            j.thr_slot.assign(1, NOTHR);
            j.thr_id.assign(1, NOTHR);
            touch(j.ref_task);
            cur_q[j.ref_task] += 1.0;
            touch(es[0]);
            cur_q[es[0]] += 1.0;
            advance(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 = thr_queue[ts].front();
                thr_queue[ts].pop_front();
                LnJob& nj = jobs[nxt];
                thr_held[ts][th] = true;
                touch(thr_owner[ts]);
                cur_u[thr_owner[ts]] += 1.0;
                nj.thr_slot.back() = ts;
                nj.thr_id.back() = th;
                advance(nxt);
                drain();
            }
            continue;
        }
 
        // A host completion.
        LnJob& j = jobs[ev.who];
        const std::size_t entry = j.stack.empty() ? 0 : j.stack.back();
        const std::vector<std::size_t>& acts =
            (entry != 0 && entry < lsn.actsof.size()) ? lsn.actsof[entry]
                                                      : std::vector<std::size_t>();
        // act_pos was advanced before the demand was started, so the activity
        // that just finished is the one before it.
        const std::size_t ai = j.act_pos.empty() ? 0 : j.act_pos.back();
        const std::size_t act = (ai > 0 && ai - 1 < acts.size()) ? acts[ai - 1] : 0;
        if (act == 0) continue;
        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, j.repl_of.empty() ? 0 : j.repl_of.back());
        touch(host);
        cur_q[host] -= 1.0;
        cur_u[host] -= 1.0;
        --host_busy[hs];
        completions[host] += 1.0;
        act_host_resid[act] += now - j.host_t0;
        // The host demand is the LAST thing an activity does, so this is where
        // the activity itself completes.
        finish_act(j, act);
 
        if (!host_queue[hs].empty()) {
            const std::size_t nxt = host_queue[hs].front();
            host_queue[hs].pop_front();
            LnJob& nj = jobs[nxt];
            const std::size_t ne = nj.stack.back();
            const std::size_t nai = nj.act_pos.back();
            const std::size_t nact = lsn.actsof[ne][nai - 1];
            ++host_busy[hs];
            cur_u[host] += 1.0;
            LnEvent e;
            e.t = now + hostdem[nact].next(g_host[nact]);
            e.kind = 0;
            e.who = nxt;
            push_ev(e);
        }
        advance(ev.who);
        drain();
    }
 
    // ---- 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;
        }
    }
 
    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