ldes_samplepath.h

Go to the documentation of this file.

// Copyright (c) 2012-2026, Imperial College London
// All rights reserved.
 
/**
 * @file ldes_samplepath.h
 * @brief The exact sample-path recorder of the native LDES engine (`--sample-path N`).
 *
 * Transcribes the reference's `SamplePathRecorder` and the `sp*` block of
 * `Solver_ssj`. The engine observes itself at every event boundary: it takes a
 * snapshot of the jobs held at each station, diffs it against the previous one
 * by job identity, and turns the difference into SolverSSA-style events plus one
 * row of LINE State per stateful node. Phase changes that fall between two real
 * events (the stages of an Erlang, the phases of a walked PH or MAP) are replayed
 * from the phase SCRIPT of each drawn interval and recorded as virtual PHASE rows.
 *
 * NO EXTRA DRAWS: recording consumes nothing from any simulation stream. A walked
 * family reports the phases it visited through `Sampler::set_trace`, which reads
 * the draws the walk takes anyway; an Erlang (drawn by convolution, which gives
 * no stage epochs) is reconstructed from k-1 sorted uniforms of the recorder's
 * own java.util.Random(seed + 900001).
 *
 * Row convention (= SolverSSA): row k is the state OCCUPIED UNTIL t[k] (the left
 * limit); events tagged step k happen at t[k].
 */
 
#ifndef LINE_SOLVERS_LDES_SAMPLEPATH_H
#define LINE_SOLVERS_LDES_SAMPLEPATH_H
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <deque>
#include <limits>
#include <map>
#include <memory>
#include <random>
#include <string>
#include <vector>
 
#include "line/solvers/ldes/ldes_sampler.h"
#include "line/util/rng_ssj.h"
 
namespace line {
namespace ldes {
namespace engine {
namespace sp {
 
/** Seed offset of the recorder's own generator (Erlang stage reconstruction). */
constexpr long long SEED_OFFSET = 900001LL;
 
/** Event wire codes (= MATLAB EventType values). */
constexpr double EV_ARV = 1.0;
constexpr double EV_DEP = 2.0;
constexpr double EV_PHASE = 3.0;
constexpr double EV_START = 15.0;
constexpr double EV_DROP = 17.0;
 
/** Java's Math.ulp for a double. */
inline double ulp(double v) {
    v = std::fabs(v);
    if (std::isnan(v)) return v;
    if (std::isinf(v)) return std::numeric_limits<double>::infinity();
    if (v == std::numeric_limits<double>::max()) return v - std::nextafter(v, 0.0);
    return std::nextafter(v, std::numeric_limits<double>::infinity()) - v;
}
 
/**
 * True for a family whose draw is a phase WALK the sampler can trace: the
 * reference's MapSampler / map_sample families, plus at a Source the BMAP and
 * marked-stationary walks. An Erlang the sampler recovered for convolution
 * (`erlang_k() > 0`) gives no stage epochs and is reconstructed instead.
 */
inline bool walked(const Sampler& s, bool arrival) {
    using lang::ProcessType;
    switch (s.type()) {
        case ProcessType::MAP:
        case ProcessType::MMPP2:
        case ProcessType::PH:
        case ProcessType::APH:
        case ProcessType::HYPEREXP:
        case ProcessType::COXIAN:
        case ProcessType::COX2:
            return true;
        case ProcessType::ERLANG:
            return s.erlang_k() == 0;
        case ProcessType::BMAP:
            return arrival;
        default:
            return arrival && lang::process_is_marked_stationary(s.type());
    }
}
 
/**
 * Phases the export reports for a law (`spPhaseCount`), or -1 when the law is
 * time-inhomogeneous and has no fixed phase space.
 */
inline int phase_count(const Sampler& s, bool disabled, bool arrival) {
    using lang::ProcessType;
    const ProcessType pt = s.type();
    if (disabled || pt == ProcessType::DISABLED) return 1;
    if (pt == ProcessType::MAPT || pt == ProcessType::PHT || lang::process_is_marked_schedule(pt)) return -1;
    if (pt == ProcessType::REPLAYER) return 1;
    if (pt == ProcessType::ERLANG && s.erlang_k() > 0) return std::max(1, s.erlang_k());
    const bool w = pt == ProcessType::PH || pt == ProcessType::APH || pt == ProcessType::HYPEREXP ||
                   pt == ProcessType::COXIAN || pt == ProcessType::COX2 || pt == ProcessType::MAP ||
                   pt == ProcessType::MMPP2 ||
                   (arrival && (pt == ProcessType::BMAP || lang::process_is_marked_stationary(pt)));
    if (w && s.map_order() >= 1) return std::max(1, static_cast<int>(s.map_order()));
    return 1;
}
 
/** Phase script of one drawn interval: segment j ends at cum[j] and is spent in ph[j] (0-based). */
struct Script {
    std::vector<double> cum;
    std::vector<int> ph;
 
    static std::shared_ptr<Script> single(double total) {
        std::shared_ptr<Script> s = std::make_shared<Script>();
        s->cum.assign(1, total);
        s->ph.assign(1, 0);
        return s;
    }
 
    /** Segment holding attained amount `a`, never below `from`. */
    int seg_at(double a, int from) const {
        int s = std::max(0, from);
        const int n = static_cast<int>(cum.size());
        while (s < n - 1 && cum[static_cast<std::size_t>(s)] <= a) ++s;
        return s;
    }
};
typedef std::shared_ptr<Script> ScriptP;
 
/** One job held at a detailed station in a snapshot. */
struct Job {
    std::uint64_t key = 0;
    std::size_t cls = 0;
    bool in_service = false;
    /** Attained amount at time t0 and its rate of growth (1 in wall time, the PS share otherwise). */
    double t0 = 0.0;
    double att0 = 0.0;
    double rate = 0.0;
    ScriptP script;
    int seg = 0;
 
    int phase() const { return script ? script->ph[static_cast<std::size_t>(seg)] : 0; }
 
    /** Epoch at which the current segment ends, +Inf if it is the last. */
    double next_boundary() const {
        if (!in_service || !script || seg >= static_cast<int>(script->cum.size()) - 1 ||
            !(rate > 0.0))
            return std::numeric_limits<double>::infinity();
        return t0 + (script->cum[static_cast<std::size_t>(seg)] - att0) / rate;
    }
};
 
/** Jobs held everywhere at one instant, indexed by station. */
struct Snap {
    /** Per station: waiting jobs head first, then the jobs in service; meaningful when detail[i]. */
    std::vector<std::vector<Job>> jobs;
    /** Per station and class: job count. */
    std::vector<std::vector<int>> counts;
    /** Per station (sources only) and class: interarrival script in progress, start, segment. */
    std::vector<std::vector<ScriptP>> src;
    std::vector<std::vector<double>> src_start;
    std::vector<std::vector<int>> src_seg;
};
 
/** One recorded event without its step and time: {type, node, class, jobId, phase}. */
struct Ev {
    double v[5];
};
 
/** Transcribes `SamplePathRecorder`: the rows and events of the path. */
class Recorder {
public:
    Recorder() = default;
    Recorder(std::size_t limit, const std::vector<std::size_t>& stateful_nodes, std::size_t nclasses,
             const std::vector<std::string>& refusal, const std::vector<bool>& refused)
        : limit_(limit), nodes_(stateful_nodes), R_(nclasses), refusal_(refusal),
          refused_(refused), state_rows_(stateful_nodes.size()), aggr_rows_(stateful_nodes.size()) {}
 
    bool is_full() const { return t_.size() >= limit_; }
    std::size_t size() const { return t_.size(); }
 
    void add_row(double time, const std::vector<std::vector<double>>& state,
                 const std::vector<std::vector<double>>& aggr, const std::vector<Ev>& evs) {
        const double step = static_cast<double>(t_.size());
        t_.push_back(time);
        for (std::size_t i = 0; i < nodes_.size(); ++i) {
            state_rows_[i].push_back(refused_[i] ? std::vector<double>() : state[i]);
            aggr_rows_[i].push_back(aggr[i]);
        }
        for (const Ev& e : evs) {
            std::vector<double> row(7);
            row[0] = step;
            row[1] = time;
            for (int c = 0; c < 5; ++c) row[static_cast<std::size_t>(c) + 2] = e.v[c];
            events_.push_back(row);
        }
    }
 
    void count_unbound_script() { ++unbound_; }
    long long unbound_scripts() const { return unbound_; }
 
    const std::vector<double>& times() const { return t_; }
    const std::vector<std::vector<double>>& events() const { return events_; }
    std::size_t num_stateful() const { return nodes_.size(); }
    std::size_t stateful_node(std::size_t isf) const { return nodes_[isf]; }
    bool refused(std::size_t isf) const { return refused_[isf]; }
    const std::string& refusal(std::size_t isf) const { return refusal_[isf]; }
 
    /** The detailed rows of node isf, LEFT-padded with zeros to the widest row. */
    std::vector<std::vector<double>> state(std::size_t isf) const {
        const std::vector<std::vector<double>>& rows = state_rows_[isf];
        std::size_t w = 0;
        for (const std::vector<double>& r : rows) w = std::max(w, r.size());
        std::vector<std::vector<double>> out(rows.size(), std::vector<double>(w, 0.0));
        for (std::size_t k = 0; k < rows.size(); ++k)
            std::copy(rows[k].begin(), rows[k].end(), out[k].begin() + static_cast<std::ptrdiff_t>(w - rows[k].size()));
        return out;
    }
    std::vector<std::vector<double>> aggr(std::size_t isf) const {
        std::vector<std::vector<double>> out(aggr_rows_[isf].size(), std::vector<double>(R_, 0.0));
        for (std::size_t k = 0; k < out.size(); ++k)
            for (std::size_t r = 0; r < R_ && r < aggr_rows_[isf][k].size(); ++r)
                out[k][r] = aggr_rows_[isf][k][r];
        return out;
    }
 
private:
    std::size_t limit_ = 0;
    std::vector<std::size_t> nodes_;
    std::size_t R_ = 0;
    std::vector<std::string> refusal_;
    std::vector<bool> refused_;
    std::vector<double> t_;
    std::vector<std::vector<double>> events_;
    std::vector<std::vector<std::vector<double>>> state_rows_, aggr_rows_;
    long long unbound_ = 0;
};
 
/**
 * The recorder's bookkeeping between two boundaries: scripts drawn and not yet
 * bound to a job, source departures and drops, the persistent job ids, and the
 * snapshot held since the last boundary. The engine fills the snapshots (it owns
 * the stations); everything that does not read station internals lives here.
 */
struct Tracker {
    // ---- configuration, set by the engine -----------------------------------
    std::size_t M = 0, R = 0;
    std::vector<long long> node_of;            ///< per station: 0-based node index
    std::vector<bool> is_source;               ///< per station
    std::vector<bool> is_service;              ///< per station: a Queue or a Delay
    std::vector<bool> detail;                  ///< per station: detailed rows encoded
    std::vector<bool> fcfs_buffer;             ///< per station: the row carries a buffer
    std::vector<std::vector<int>> phases;      ///< per station and class
    std::vector<std::vector<int>> map_vars;    ///< per station and class
    std::vector<std::vector<double>> map_var_init;
    std::vector<long long> station_of_sf;      ///< per stateful node: station or -1
    std::vector<std::vector<double>> const_row; ///< per stateful node that is no station
 
    // ---- run state ------------------------------------------------------------
    Recorder rec;
    bool done = false;
    bool have_held = false;
    Snap held;
    rng::JavaRandom rng;
    PhaseTrace scratch;
    bool traced = false;
    struct Pending {
        std::size_t q, k;
        ScriptP s;
    };
    std::vector<Pending> pending;
    std::vector<std::vector<ScriptP>> src_script;   ///< [station][class]
    std::vector<std::vector<double>> src_start;
    std::vector<std::pair<std::size_t, std::size_t>> src_deps, drops;
    std::map<std::uint64_t, ScriptP> scripts;
    std::map<std::uint64_t, long long> ids;
    std::deque<long long> id_fifo;
    long long next_id = 0;
    /** True when the recorder stopped the run (stoppingReason "sample_path"). */
    bool stopped_run = false;
    double stop_time = 0.0;
 
    void seed(long long s) {
        if (s > 0) {
            rng.set_seed(s + SEED_OFFSET);
        } else {
            std::random_device rd;
            rng.set_seed((static_cast<long long>(rd()) << 32) ^ static_cast<long long>(rd()));
        }
    }
 
    /** Hands the scratch trace to a recording sampler call. */
    PhaseTrace* trace_fill() {
        scratch.clear();
        traced = true;
        return &scratch;
    }
 
    /** Script of the draw that just returned `value`: its walk, an Erlang reconstruction, or one segment. */
    ScriptP script_of(double value, int erlang_k) {
        const bool tr = traced;
        traced = false;
        if (tr && scratch.size() > 0) {
            const double tot = scratch.total();
            if (std::fabs(tot - value) <= 1e-9 * std::fabs(value) + 8.0 * ulp(value)) {
                const std::size_t n = scratch.size();
                ScriptP s = std::make_shared<Script>();
                s->cum.assign(n, 0.0);
                s->ph.assign(n, 0);
                double acc = 0.0;
                for (std::size_t i = 0; i < n; ++i) {
                    acc += scratch.dt[i];
                    s->cum[i] = acc;
                    s->ph[i] = scratch.phase[i];
                }
                s->cum[n - 1] = value;
                return s;
            }
            return Script::single(value);
        }
        if (erlang_k > 1) {
            const std::size_t k = static_cast<std::size_t>(erlang_k);
            std::vector<double> u(k - 1);
            for (std::size_t i = 0; i + 1 < k; ++i) u[i] = rng.next_double();
            std::sort(u.begin(), u.end());
            ScriptP s = std::make_shared<Script>();
            s->cum.assign(k, 0.0);
            s->ph.assign(k, 0);
            for (std::size_t i = 0; i + 1 < k; ++i) {
                s->cum[i] = u[i] * value;
                s->ph[i] = static_cast<int>(i);
            }
            s->cum[k - 1] = value;
            s->ph[k - 1] = static_cast<int>(k - 1);
            return s;
        }
        return Script::single(value);
    }
 
    void capture_service(std::size_t q, std::size_t k, double value, int erlang_k) {
        Pending p;
        p.q = q;
        p.k = k;
        p.s = script_of(value, erlang_k);
        pending.push_back(p);
    }
    void capture_arrival(std::size_t st, std::size_t k, double value, int erlang_k, double now) {
        src_script[st][k] = script_of(value, erlang_k);
        src_start[st][k] = now;
    }
    void source_departure(std::size_t st, std::size_t k) { src_deps.push_back(std::make_pair(st, k)); }
    void drop(std::size_t q, std::size_t k) { drops.push_back(std::make_pair(q, k)); }
 
    /** Binds a drawn script of (q, k) to a job whose drawn total is `expected`. */
    ScriptP bind(std::size_t q, std::size_t k, double expected, double now) {
        long best = -1;
        double best_err = std::numeric_limits<double>::infinity();
        for (std::size_t i = 0; i < pending.size(); ++i) {
            if (pending[i].q != q || pending[i].k != k) continue;
            const Script& s = *pending[i].s;
            const double err = std::fabs(s.cum.back() - expected);
            if (err < best_err) {
                best_err = err;
                best = static_cast<long>(i);
            }
        }
        const double tol = 1e-9 * std::fabs(expected) + 8.0 * ulp(std::fabs(now) + std::fabs(expected));
        if (best >= 0 && best_err <= tol) {
            ScriptP s = pending[static_cast<std::size_t>(best)].s;
            pending.erase(pending.begin() + best);
            const std::size_t n = s->cum.size();
            s->cum[n - 1] = std::max(expected, n > 1 ? s->cum[n - 2] : 0.0);
            return s;
        }
        rec.count_unbound_script();
        return Script::single(expected);
    }
 
    /** Script already bound to `key`, or a fresh binding (stored under `key`). */
    ScriptP script_for(std::uint64_t key, std::size_t q, std::size_t k, double expected, double now) {
        std::map<std::uint64_t, ScriptP>::const_iterator it = scripts.find(key);
        if (it != scripts.end()) return it->second;
        ScriptP s = bind(q, k, expected, now);
        scripts[key] = s;
        return s;
    }
 
    /** Source part of a snapshot. */
    void snap_sources(Snap& s, double now, const Snap* prev) const {
        s.src.assign(M, std::vector<ScriptP>(R));
        s.src_start.assign(M, std::vector<double>(R, 0.0));
        s.src_seg.assign(M, std::vector<int>(R, 0));
        for (std::size_t i = 0; i < M; ++i) {
            if (!is_source[i]) continue;
            for (std::size_t k = 0; k < R; ++k) {
                const ScriptP& sc = src_script[i][k];
                s.src[i][k] = sc;
                s.src_start[i][k] = src_start[i][k];
                if (sc) {
                    const int from = (prev && prev->src[i][k] == sc) ? prev->src_seg[i][k] : 0;
                    s.src_seg[i][k] = sc->seg_at(now - src_start[i][k], from);
                }
            }
        }
    }
 
    /**
     * Encodes snapshot `s` as LINE State rows and per-class counts. `sampler_phase(i, k)`
     * is the service sampler's walk phase (the reference's `getPhase()`), read for the
     * MAP local variable of an idle class.
     */
    template <class PhaseFn>
    void encode(const Snap& s, std::vector<std::vector<double>>& state,
                std::vector<std::vector<double>>& aggr, PhaseFn sampler_phase) const {
        const std::size_t nsf = station_of_sf.size();
        state.assign(nsf, std::vector<double>());
        aggr.assign(nsf, std::vector<double>());
        const double inf = std::numeric_limits<double>::infinity();
        for (std::size_t isf = 0; isf < nsf; ++isf) {
            const long long ist = station_of_sf[isf];
            std::vector<double> ag(R, 0.0);
            std::vector<double> row;
            if (ist >= 0 && is_source[static_cast<std::size_t>(ist)]) {
                const std::size_t i = static_cast<std::size_t>(ist);
                std::size_t w = 1;
                for (std::size_t k = 0; k < R; ++k) {
                    w += static_cast<std::size_t>(phases[i][k]);
                    ag[k] = inf;
                }
                row.assign(w, 0.0);
                row[0] = inf;
                std::size_t off = 1;
                for (std::size_t k = 0; k < R; ++k) {
                    const ScriptP& sc = s.src[i][k];
                    if (sc) {
                        const int p = std::min(sc->ph[static_cast<std::size_t>(s.src_seg[i][k])], phases[i][k] - 1);
                        row[off + static_cast<std::size_t>(p)] = 1.0;
                    }
                    off += static_cast<std::size_t>(phases[i][k]);
                }
            } else if (ist >= 0 && is_service[static_cast<std::size_t>(ist)]) {
                const std::size_t i = static_cast<std::size_t>(ist);
                for (std::size_t k = 0; k < R; ++k) ag[k] = s.counts[i][k];
                if (detail[i]) {
                    std::size_t nsrv = 0;
                    std::vector<std::size_t> base(R, 0);
                    for (std::size_t k = 0; k < R; ++k) {
                        base[k] = nsrv;
                        nsrv += static_cast<std::size_t>(phases[i][k]);
                    }
                    std::vector<double> srv(nsrv, 0.0);
                    std::vector<double> buf;
                    std::vector<double> map_var(R, 0.0);
                    std::vector<bool> map_busy(R, false);
                    for (const Job& j : s.jobs[i]) {
                        if (j.in_service) {
                            const int p = std::min(j.phase(), phases[i][j.cls] - 1);
                            srv[base[j.cls] + static_cast<std::size_t>(p)] += 1.0;
                            map_var[j.cls] = p + 1;
                            map_busy[j.cls] = true;
                        } else {
                            buf.push_back(static_cast<double>(j.cls + 1));
                        }
                    }
                    const std::size_t nbuf = fcfs_buffer[i] ? std::max<std::size_t>(1, buf.size()) : 0;
                    std::size_t nvar = 0;
                    for (std::size_t k = 0; k < R; ++k) nvar += static_cast<std::size_t>(map_vars[i][k]);
                    row.assign(nbuf + nsrv + nvar, 0.0);
                    for (std::size_t b = 0; b < buf.size(); ++b) row[nbuf - 1 - b] = buf[b];
                    std::copy(srv.begin(), srv.end(), row.begin() + static_cast<std::ptrdiff_t>(nbuf));
                    std::size_t off = nbuf + nsrv;
                    for (std::size_t k = 0; k < R; ++k) {
                        if (map_vars[i][k] > 0) {
                            double v = map_var[k];
                            if (!map_busy[k]) {
                                const int ph = sampler_phase(i, k);
                                v = (ph >= 0) ? ph + 1 : map_var_init[i][k];
                            }
                            row[off++] = v;
                        }
                    }
                }
            } else {
                row = const_row[isf];
            }
            state[isf] = row;
            aggr[isf] = ag;
        }
    }
 
    long long take_id() {
        if (id_fifo.empty()) return next_id++;
        const long long id = id_fifo.front();
        id_fifo.pop_front();
        return (id < 0) ? next_id++ : id;
    }
 
    /** Events turning snapshot `a` into `b`, in the order DEP, ARV, DROP, START, PHASE. */
    std::vector<Ev> diff(const Snap& a, const Snap& b) {
        std::vector<Ev> deps, arvs, drps, starts, phs;
        auto mk = [](double t, double n, double c, double id, double ph) {
            Ev e;
            e.v[0] = t;
            e.v[1] = n;
            e.v[2] = c;
            e.v[3] = id;
            e.v[4] = ph;
            return e;
        };
        for (std::size_t q = 0; q < M; ++q) {
            if (!is_service[q]) continue;
            const double node = static_cast<double>(node_of[q]);
            if (!detail[q]) {
                for (std::size_t k = 0; k < R; ++k)
                    for (int d = b.counts[q][k]; d < a.counts[q][k]; ++d) {
                        deps.push_back(mk(EV_DEP, node, static_cast<double>(k), -1.0, 0.0));
                        id_fifo.push_back(-1);
                    }
                continue;
            }
            std::map<std::uint64_t, const Job*> bm;
            for (const Job& j : b.jobs[q]) bm[j.key] = &j;
            for (const Job& j : a.jobs[q]) {
                if (bm.count(j.key)) continue;
                long long jid = -1;
                std::map<std::uint64_t, long long>::iterator it = ids.find(j.key);
                if (it != ids.end()) {
                    jid = it->second;
                    ids.erase(it);
                }
                scripts.erase(j.key);
                deps.push_back(mk(EV_DEP, node, static_cast<double>(j.cls), static_cast<double>(jid),
                                  j.in_service ? j.phase() + 1.0 : 0.0));
                id_fifo.push_back(jid);
            }
        }
        for (const std::pair<std::size_t, std::size_t>& sd : src_deps) {
            const long long id = next_id++;
            deps.push_back(mk(EV_DEP, static_cast<double>(node_of[sd.first]), static_cast<double>(sd.second),
                              static_cast<double>(id), 0.0));
            id_fifo.push_back(id);
        }
        for (std::size_t q = 0; q < M; ++q) {
            if (!is_service[q]) continue;
            const double node = static_cast<double>(node_of[q]);
            if (!detail[q]) {
                for (std::size_t k = 0; k < R; ++k)
                    for (int d = a.counts[q][k]; d < b.counts[q][k]; ++d) {
                        take_id();
                        arvs.push_back(mk(EV_ARV, node, static_cast<double>(k), -1.0, 0.0));
                    }
                continue;
            }
            std::map<std::uint64_t, const Job*> am;
            for (const Job& j : a.jobs[q]) am[j.key] = &j;
            for (const Job& j : b.jobs[q]) {
                std::map<std::uint64_t, const Job*>::const_iterator it = am.find(j.key);
                const Job* old = (it == am.end()) ? nullptr : it->second;
                if (!old) {
                    const long long id = take_id();
                    ids[j.key] = id;
                    arvs.push_back(mk(EV_ARV, node, static_cast<double>(j.cls), static_cast<double>(id), 0.0));
                }
                const double id = static_cast<double>(ids[j.key]);
                if (j.in_service && (!old || !old->in_service))
                    starts.push_back(mk(EV_START, node, static_cast<double>(j.cls), id, j.phase() + 1.0));
                else if (j.in_service && old->phase() != j.phase())
                    phs.push_back(mk(EV_PHASE, node, static_cast<double>(j.cls), id, j.phase() + 1.0));
            }
        }
        for (const std::pair<std::size_t, std::size_t>& dr : drops)
            drps.push_back(mk(EV_DROP, static_cast<double>(node_of[dr.first]), static_cast<double>(dr.second),
                              static_cast<double>(take_id()), 0.0));
        for (std::size_t i = 0; i < M; ++i) {
            if (!is_source[i]) continue;
            for (std::size_t k = 0; k < R; ++k) {
                const ScriptP& sc = b.src[i][k];
                if (sc && a.src[i][k] == sc && a.src_seg[i][k] != b.src_seg[i][k])
                    phs.push_back(mk(EV_PHASE, static_cast<double>(node_of[i]), static_cast<double>(k), -1.0,
                                     sc->ph[static_cast<std::size_t>(b.src_seg[i][k])] + 1.0));
            }
        }
        std::vector<Ev> evs(deps);
        evs.insert(evs.end(), arvs.begin(), arvs.end());
        evs.insert(evs.end(), drps.begin(), drps.end());
        evs.insert(evs.end(), starts.begin(), starts.end());
        evs.insert(evs.end(), phs.begin(), phs.end());
        return evs;
    }
 
    /**
     * Event-boundary hook (`spBoundary`): records what the event that just ran
     * changed, then, unless the run is stopping, the phase changes due before
     * `t_next`. `snapshot(now, prev)` builds the current snapshot.
     *
     * @return 0 to continue; 1 when the recorder filled at `now` and the run stops
     *         there; 2 when it filled at a virtual row and the run stops at
     *         `stop_time` (> now), before the next event.
     */
    template <class SnapFn, class PhaseFn>
    int boundary(bool stopped, double now, double t_next, SnapFn snapshot, PhaseFn sampler_phase) {
        if (!have_held) {
            held = snapshot(now, nullptr);
            have_held = true;
            for (std::size_t i = 0; i < M; ++i)
                if (is_service[i] && detail[i])
                    for (const Job& j : held.jobs[i]) ids[j.key] = next_id++;
        } else {
            Snap cur = snapshot(now, &held);
            const std::vector<Ev> evs = diff(held, cur);
            if (!evs.empty()) {
                std::vector<std::vector<double>> st, ag;
                encode(held, st, ag, sampler_phase);
                rec.add_row(now, st, ag, evs);
            }
            held = cur;
        }
        pending.clear();
        src_deps.clear();
        drops.clear();
        id_fifo.clear();
        if (rec.is_full()) {
            done = true;
            if (!stopped) {
                stopped_run = true;
                stop_time = now;
                return 1;
            }
            return 0;
        }
        if (stopped) {
            done = true;
            return 0;
        }
        // Virtual PHASE rows: scripted phase changes due before the next event.
        while (true) {
            double tau = std::numeric_limits<double>::infinity();
            for (std::size_t i = 0; i < M; ++i)
                if (is_service[i] && detail[i])
                    for (const Job& j : held.jobs[i]) tau = std::min(tau, j.next_boundary());
            for (std::size_t i = 0; i < M; ++i) {
                if (!is_source[i]) continue;
                for (std::size_t k = 0; k < R; ++k) {
                    const ScriptP& sc = held.src[i][k];
                    const int sg = held.src_seg[i][k];
                    if (sc && sg < static_cast<int>(sc->cum.size()) - 1)
                        tau = std::min(tau, held.src_start[i][k] + sc->cum[static_cast<std::size_t>(sg)]);
                }
            }
            if (!(tau < t_next)) break;
            tau = std::max(tau, now);
            std::vector<std::vector<double>> st, ag;
            encode(held, st, ag, sampler_phase);
            std::vector<Ev> evs;
            for (std::size_t i = 0; i < M; ++i) {
                if (!(is_service[i] && detail[i])) continue;
                const double node = static_cast<double>(node_of[i]);
                for (Job& j : held.jobs[i]) {
                    const double tb = j.next_boundary();
                    if (tb <= tau) {
                        const int before = j.phase();
                        ++j.seg;
                        while (j.next_boundary() <= tau) ++j.seg;
                        if (j.phase() != before) {
                            Ev e;
                            e.v[0] = EV_PHASE;
                            e.v[1] = node;
                            e.v[2] = static_cast<double>(j.cls);
                            e.v[3] = static_cast<double>(ids[j.key]);
                            e.v[4] = j.phase() + 1.0;
                            evs.push_back(e);
                        }
                    }
                }
            }
            for (std::size_t i = 0; i < M; ++i) {
                if (!is_source[i]) continue;
                for (std::size_t k = 0; k < R; ++k) {
                    const ScriptP& sc = held.src[i][k];
                    if (!sc) continue;
                    int& sg = held.src_seg[i][k];
                    const int before = sc->ph[static_cast<std::size_t>(sg)];
                    while (sg < static_cast<int>(sc->cum.size()) - 1 &&
                           held.src_start[i][k] + sc->cum[static_cast<std::size_t>(sg)] <= tau)
                        ++sg;
                    if (sc->ph[static_cast<std::size_t>(sg)] != before) {
                        Ev e;
                        e.v[0] = EV_PHASE;
                        e.v[1] = static_cast<double>(node_of[i]);
                        e.v[2] = static_cast<double>(k);
                        e.v[3] = -1.0;
                        e.v[4] = sc->ph[static_cast<std::size_t>(sg)] + 1.0;
                        evs.push_back(e);
                    }
                }
            }
            if (evs.empty()) continue;  // a segment boundary that kept the phase (a self-loop sojourn)
            rec.add_row(tau, st, ag, evs);
            if (rec.is_full()) {
                // The last row is stamped tau, later than the clock: stop the run AT
                // tau, so the statistics integrate the same interval the path covers.
                done = true;
                stopped_run = true;
                stop_time = tau;
                return 2;
            }
        }
        return 0;
    }
};
 
}  // namespace sp
}  // namespace engine
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_SAMPLEPATH_H