ldes_sampler.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_SAMPLER_H
#define LINE_SOLVERS_LDES_LDES_SAMPLER_H
 
/**
 * @file
 * @ingroup line_solvers
 * The variate generators of the native LDES engine.
 *
 * ONE SAMPLER PER (station, class) PAIR, built once and then advanced, which is
 * what `Solver_ssj.initializeGenerators` does and is not an optimization: a
 * MAP, a RAP and an ME carry a PHASE across successive samples, and rebuilding
 * the generator per variate would restart that phase every time and silently
 * turn a correlated process into a renewal one with the same marginal. The
 * autocorrelation is the reason those processes are in the model at all, so
 * losing it produces a run that looks converged and answers a different
 * question.
 *
 * WHERE THE PARAMETERS COME FROM, transcribed from
 * `createNonMarkovianArrivalGen` and `firingGenFromMeanScv`:
 *
 *  - ERLANG, HYPEREXP, PH, APH, COXIAN, COX2, MAP, MMPP2, ME, RAP, DMAP read
 *    the (D0,D1) REPRESENTATION, so their higher moments are exactly the
 *    model's;
 *  - UNIFORM, GAMMA, WEIBULL, LOGNORMAL, PARETO read the MOMENT PAIR (mean,
 *    SCV) and invert it, because `sn.proc` carries an Erlang FIT for them
 *    rather than their own parameters (see _kb/09-ldes-and-cache.md);
 *  - DET, IMMEDIATE, REPLAYER, and the counting families read their own
 *    parameters directly.
 *
 * A moment pair that no member of the family can realise is REFUSED by name.
 * A non-negative uniform needs SCV <= 1/3; papering over that with a clamp
 * would run a model with a different variance and report it as the user's.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <memory>
#include <random>
#include <string>
#include <vector>
 
#include "line/api/mam/map_sample.h"
#include "line/api/mam/me_sample.h"
#include "line/solvers/ldes/ldes_ssj_variates.h"
#include "line/util/rng_ssj.h"
#include "line/lang/lang_types.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace ldes {
namespace engine {
 
/**
 * The engine's randomness, in the SHAPE the Java engine uses it.
 *
 * `Solver_ssj` carries TWO generators wherever it samples, and which one serves
 * depends on the family rather than on the call site:
 *
 *     MRG32k3a stream = new MRG32k3a();
 *     stream.setSeed(new long[]{ seed+offset, ..., seed+offset+5 });
 *     this.svcRng[i][k] = new java.util.Random(seed + offset);
 *
 * The SSJ `randvar` generators (Exponential, Erlang, Uniform, Weibull, Pareto,
 * Lognormal, Gamma, Poisson, Binomial, Bernoulli) are built on the MRG STREAM;
 * the Markovian families do not go through SSJ at all -- `MapSampleGen` calls
 * `Map_sample.map_sample(D0, D1, 1, rng)` with the java.util.Random. So a C++
 * engine that means to walk the same sample path needs both, seeded from the
 * same `seed + offset`, and must send each family to the same one.
 *
 * `mc` is the third generator and the one still out of step: `rap_sample` and
 * `me_sample` take `pfqn::McRng` by reference, and the RAP/ME path has no SSJ
 * or java.util.Random counterpart to align to yet.
 */
struct Rng {
    rng::Mrg32k3a stream;  ///< SSJ MRG32k3a: every renewal family
    rng::JavaRandom aux;   ///< java.util.Random: the Markovian samplers
    pfqn::McRng mc;        ///< the residual RAP/ME path, not yet aligned
 
    /**
     * A run seed and the TWO offsets the reference derives, one per generator.
     *
     * They are not always equal, which is why this takes both. At a SERVICE
     * site `Solver_ssj` seeds the stream and the java.util.Random from the same
     * `((numSources + svcIdx) * numClasses + k) * 10 + 1000`; at an ARRIVAL site
     * the stream gets `(srcIdx * numClasses + k) * 10` and the Random the same
     * expression PLUS 2000. Seeding both from one offset would put the
     * Markovian arrival samplers on a stream the reference never uses.
     */
    Rng(long long seed, long long stream_offset, long long aux_offset)
        : aux(seed + aux_offset), mc(static_cast<std::uint64_t>(seed + aux_offset)) {
        stream.set_seed_offset(seed, stream_offset);
    }
 
    /** The common case, where the reference uses one offset for both. */
    Rng(long long seed, long long offset) : Rng(seed, offset, offset) {}
};
 
/** Uniform on (0,1) off the MRG stream; `next_double` never returns 0. */
inline double uniform01(Rng& g) { return g.stream.next_double(); }
 
/**
 * One variate generator, holding whatever state its family needs.
 *
 * Deliberately a tagged struct rather than a class hierarchy: the sampler is
 * called once per event on the engine's hot path, and a virtual dispatch there
 * costs more than the switch. The tag IS `ProcessType`, so a family added to
 * the language shows up here as a missing case rather than as silence.
 */
class Sampler {
public:
    Sampler() {}
 
    /** Build the generator of `d`; the (station, class) names are for diagnostics. */
    template <class T>
    Sampler(const lang::Distrib<T>& d, const std::string& where) : where_(where) {
        type_ = d.type;
        mean_ = num_traits<T>::to_double(d.mean);
        scv_ = num_traits<T>::to_double(d.scv);
        switch (type_) {
            case lang::ProcessType::DISABLED:
                break;
            case lang::ProcessType::EXP:
                require_mean();
                rate_ = 1.0 / mean_;
                break;
            case lang::ProcessType::DET:
                require_mean();
                break;
            case lang::ProcessType::IMMEDIATE:
                // The reference serves an Immediate in 1e-8 time units, not in
                // zero: a zero service time is an infinite rate the estimators
                // cannot carry. See Distrib::immediate().
                mean_ = 1.0 / lang::GlobalConstants::Immediate;
                break;
            case lang::ProcessType::REPLAYER:
                for (const T& v : d.trace) trace_.push_back(num_traits<T>::to_double(v));
                if (trace_.empty())
                    throw InputError("SolverLDES (native engine): the Replayer at " + where_ +
                                     " carries no samples");
                break;
            case lang::ProcessType::UNIFORM: {
                require_moments();
                const double half = mean_ * std::sqrt(3.0 * scv_);
                a_ = mean_ - half;
                b_ = mean_ + half;
                if (!(a_ >= 0.0) || !(b_ > a_))
                    throw InputError("SolverLDES (native engine): the Uniform at " + where_ +
                                     " implies a support reaching below zero; a non-negative "
                                     "uniform requires an SCV of at most 1/3");
                break;
            }
            case lang::ProcessType::GAMMA:
                require_moments();
                a_ = 1.0 / scv_;          // shape
                b_ = mean_ * scv_;        // SCALE, not the reference's rate
                break;
            case lang::ProcessType::WEIBULL: {
                require_moments();
                const double c = std::sqrt(scv_);
                a_ = std::pow(c, -1.086);                      // shape
                b_ = mean_ / std::tgamma(1.0 + 1.0 / a_);      // scale
                if (!(a_ > 0.0) || !(b_ > 0.0))
                    throw InputError("SolverLDES (native engine): the Weibull at " + where_ +
                                     " has unusable moments");
                break;
            }
            case lang::ProcessType::LOGNORMAL: {
                require_moments();
                const double c2p1 = scv_ + 1.0;
                a_ = std::log(mean_ / std::sqrt(c2p1));  // mu
                b_ = std::sqrt(std::log(c2p1));          // sigma
                if (!(b_ > 0.0))
                    throw InputError("SolverLDES (native engine): the Lognormal at " + where_ +
                                     " has a zero log-variance");
                break;
            }
            case lang::ProcessType::PARETO:
                require_moments();
                // SCV = 1/(a(a-2)) for a > 2, so a = 1 + sqrt(1 + 1/SCV), and the
                // scale follows from mean = a*m/(a-1).
                a_ = std::sqrt(1.0 + 1.0 / scv_) + 1.0;  // shape
                b_ = mean_ * (a_ - 1.0) / a_;            // scale
                break;
            case lang::ProcessType::GEOMETRIC:
                require_mean();
                // Supported on {1,2,...}: the trial index of the first success,
                // mean 1/p and SCV 1-p. SSJ's own Geometric counts FAILURES and
                // would admit a zero-length interval.
                a_ = 1.0 / mean_;
                if (!(a_ > 0.0) || a_ > 1.0)
                    throw InputError("SolverLDES (native engine): the Geometric at " + where_ +
                                     " has a success probability outside (0,1]");
                break;
            case lang::ProcessType::POISSON:
                require_mean();
                break;
            case lang::ProcessType::BERNOULLI:
                if (!(mean_ >= 0.0) || !(mean_ <= 1.0))
                    throw InputError("SolverLDES (native engine): the Bernoulli at " + where_ +
                                     " has a mean outside [0,1]");
                break;
            case lang::ProcessType::BINOMIAL: {
                require_moments();
                const double p = 1.0 - scv_ * mean_;
                if (!(p > 0.0) || p > 1.0)
                    throw InputError("SolverLDES (native engine): the Binomial at " + where_ +
                                     " has moments no (n,p) pair realises");
                a_ = p;
                b_ = std::floor(mean_ / p + 0.5);
                break;
            }
            case lang::ProcessType::MAPT:
            case lang::ProcessType::PHT:
            case lang::ProcessType::NHPP:
                load_schedule(d);
                break;
            case lang::ProcessType::ERLANG:
            case lang::ProcessType::HYPEREXP:
            case lang::ProcessType::PH:
            case lang::ProcessType::APH:
            case lang::ProcessType::COXIAN:
            case lang::ProcessType::COX2:
            case lang::ProcessType::MAP:
            case lang::ProcessType::MMPP2:
            case lang::ProcessType::ME:
            case lang::ProcessType::RAP:
            case lang::ProcessType::DMAP:
            case lang::ProcessType::MMAP:
            case lang::ProcessType::BMAP:
                load_map(d);
                break;
            default:
                throw UnsupportedError("SolverLDES (native engine): the process at " + where_ +
                                       " is of a family this engine does not sample");
        }
    }
 
    bool disabled() const { return type_ == lang::ProcessType::DISABLED; }
    lang::ProcessType type() const { return type_; }
    double mean() const { return mean_; }
 
    /**
     * One variate from a TIME-INHOMOGENEOUS process, given the instant it
     * starts at.
     *
     * The families whose parameters depend on absolute time cannot be sampled
     * from a duration alone, so this overload takes `from` and every other
     * family ignores it. The engine calls it wherever a duration is drawn; a
     * homogeneous process is unaffected.
     */
    double next_at(Rng& g, double from) {
        if (!has_schedule_) return next(g);
        return schedule_sample(g, from);
    }
 
    bool time_varying() const { return has_schedule_; }
 
    /**
     * One variate. Advances whatever phase the family carries.
     *
     * EVERY RENEWAL FAMILY BUT THE ERLANG COSTS EXACTLY ONE UNIFORM off the MRG
     * stream, which is not an optimization but the reference's behaviour: SSJ's
     * instance generators invert, and a family that drew twice would shift
     * every later event in the run. The quantiles live in
     * `ldes_ssj_variates.h` and are pinned there against SSJ's own output. The
     * Erlang is the exception, and deliberately so: it costs k uniforms in both
     * engines because the Gamma inversion SSJ used is not reproducible across
     * JVM releases.
     */
    double next(Rng& g) {
        switch (type_) {
            case lang::ProcessType::EXP:
                // SSJ's ExponentialDist.inverseF is -log1p(-u)/lambda, which is
                // NOT -log(u)/lambda: they differ in the last bits for small u,
                // exactly where an interarrival time matters most.
                return ssj::exponential_inverse(rate_, uniform01(g));
            case lang::ProcessType::DET:
            case lang::ProcessType::IMMEDIATE:
                return mean_;
            case lang::ProcessType::REPLAYER: {
                const double v = trace_[trace_idx_];
                trace_idx_ = (trace_idx_ + 1) % trace_.size();
                // The reference floors a non-positive trace entry rather than
                // scheduling an event in the past.
                return (v <= 0.0) ? 1e-9 : v;
            }
            case lang::ProcessType::UNIFORM:
                return ssj::uniform_inverse(a_, b_, uniform01(g));
            case lang::ProcessType::GAMMA:
                // b_ is the SCALE here; SSJ's GammaGen takes the RATE.
                return ssj::gamma_inverse(a_, 1.0 / b_, uniform01(g));
            case lang::ProcessType::WEIBULL:
                // b_ is the SCALE; SSJ's WeibullGen takes lambda = 1/scale.
                return ssj::weibull_inverse(a_, 1.0 / b_, 0.0, uniform01(g));
            case lang::ProcessType::LOGNORMAL:
                return ssj::lognormal_inverse(a_, b_, uniform01(g));
            case lang::ProcessType::PARETO:
                return ssj::pareto_inverse(a_, b_, uniform01(g));
            case lang::ProcessType::GEOMETRIC: {
                if (a_ >= 1.0) {
                    uniform01(g);  // keep the stream consumption independent of p
                    return 1.0;
                }
                return std::ceil(std::log(1.0 - uniform01(g)) / std::log(1.0 - a_));
            }
            case lang::ProcessType::POISSON:
                return ssj::poisson_inverse(mean_, uniform01(g));
            case lang::ProcessType::BERNOULLI:
                return ssj::bernoulli_inverse(mean_, uniform01(g));
            case lang::ProcessType::BINOMIAL:
                return ssj::binomial_inverse(static_cast<int>(b_), a_, uniform01(g));
            case lang::ProcessType::ERLANG: {
                // An Erlang costs k uniforms, one per phase, which is the one
                // renewal family that does not invert in a single draw. The
                // reference used to invert the Gamma here, but that inversion
                // is iterative and its last bits move between JVM releases, so
                // a seeded run was not reproducible across JVMs; both engines
                // now convolve exponentials instead.
                if (erlang_k_ <= 0) return next_markovian(g);
                double sum = 0.0;
                for (int i = 0; i < erlang_k_; ++i)
                    sum += ssj::exponential_inverse(erlang_rate_, uniform01(g));
                return sum;
            }
            default:
                return next_markovian(g);
        }
    }
 
    /**
     * The phase the process is in, for a caller that must carry it across a
     * preemption or across successive executions of an LQN entry. -1 when the
     * family has no phase.
     */
    int phase() const { return has_phase_ ? static_cast<int>(phase_) : -1; }
    void set_phase(int p) {
        if (has_phase_ && p >= 0 && static_cast<std::size_t>(p) < map_.order()) {
            phase_ = static_cast<std::size_t>(p);
            phase_known_ = true;
        }
    }
 
private:
    void require_mean() const {
        if (!(mean_ > 0.0) || !std::isfinite(mean_))
            throw InputError("SolverLDES (native engine): the process at " + where_ +
                             " has a mean that is neither finite nor positive");
    }
    void require_moments() const {
        require_mean();
        if (!(scv_ > 0.0) || !std::isfinite(scv_))
            throw InputError("SolverLDES (native engine): the process at " + where_ +
                             " has an SCV that is neither finite nor positive");
    }
 
    template <class T>
    void load_map(const lang::Distrib<T>& d) {
        if (d.D0.rows() == 0 || d.D0.rows() != d.D1.rows())
            throw InputError("SolverLDES (native engine): the process at " + where_ +
                             " declares no (D0,D1) representation to sample");
        const std::size_t K = d.D0.rows();
        map_.D0 = Matrix<double>(K, K, 0.0);
        map_.D1 = Matrix<double>(K, K, 0.0);
        for (std::size_t i = 0; i < K; ++i)
            for (std::size_t j = 0; j < K; ++j) {
                map_.D0(i, j) = num_traits<T>::to_double(d.D0(i, j));
                map_.D1(i, j) = num_traits<T>::to_double(d.D1(i, j));
            }
        // THE MOMENTS OF A (D0,D1) PROCESS COME FROM THE PAIR, not from the
        // Distrib's own fields: `Distrib::map_dist` leaves mean 0 and SCV 1 as
        // placeholders, so reading them would give this station a zero mean
        // service time and a utilization computed against nothing.
        if (!(mean_ > 0.0)) {
            mean_ = num_traits<T>::to_double(mam::map_mean(map_));
            if (!(mean_ > 0.0))
                throw InputError("SolverLDES (native engine): the process at " + where_ +
                                 " has a non-positive mean under its (D0,D1) representation");
        }
        // ONLY the genuinely correlated families carry a phase between samples.
        // A PH, an Erlang, a HyperExp and a Coxian are RENEWAL: restarting them
        // from the entry law each time is what makes successive services
        // independent, and carrying the phase instead would introduce a
        // correlation the model does not have.
        has_phase_ = (type_ == lang::ProcessType::MAP || type_ == lang::ProcessType::MMPP2 ||
                      type_ == lang::ProcessType::RAP || type_ == lang::ProcessType::MMAP ||
                      type_ == lang::ProcessType::BMAP || type_ == lang::ProcessType::DMAP);
        // ONCE PER STATION, not once per event: the table costs a matrix
        // exponential and a thousand row-vector products, and every service at
        // this station then inverts by binary search.
        if (type_ == lang::ProcessType::ME)
            me_.reset(new mam::MeSampler<double>(map_));
        if (type_ == lang::ProcessType::ERLANG) recover_erlang();
    }
 
    /**
     * Recovers (k, lambda) from an Erlang's (D0,D1) so `next` can convolve k
     * exponentials off the MRG stream instead of walking the phases off the
     * auxiliary java.util.Random. The reference engine samples the same
     * convolution, so the two engines then agree sample for sample on an
     * Erlang model; walking the phases does not, because it draws 2k uniforms
     * from a different stream. A cell that does not carry the bidiagonal
     * equal-rate structure of an Erlang is left to the MAP walk rather than
     * approximated.
     */
    void recover_erlang() {
        erlang_k_ = 0;
        const std::size_t n = map_.order();
        if (n == 0) return;
        const double lambda = -map_.D0(0, 0);
        if (!(lambda > 0.0)) return;
        const double tol = 1e-12 * lambda;
        for (std::size_t i = 0; i < n; ++i) {
            if (std::fabs(-map_.D0(i, i) - lambda) > tol) return;
            for (std::size_t j = 0; j < n; ++j) {
                if (j == i) continue;
                const double want = (j == i + 1) ? lambda : 0.0;
                if (std::fabs(map_.D0(i, j) - want) > tol) return;
            }
            // Only the last phase completes, and it completes into phase 0.
            for (std::size_t j = 0; j < n; ++j) {
                const double want = (i + 1 == n && j == 0) ? lambda : 0.0;
                if (std::fabs(map_.D1(i, j) - want) > tol) return;
            }
        }
        erlang_k_ = static_cast<int>(n);
        erlang_rate_ = lambda;
    }
 
    /**
     * One interval from a MAP, transcribed from `Map_sample.MapSampler.next`.
     *
     * THE GENERATOR IS java.util.Random, NOT the MRG stream. `MapSampleGen`
     * extends SSJ's `RandomVariateGen` and is handed the stream like every other
     * generator, but its `nextDouble` ignores it and calls
     * `Map_sample.map_sample(D0, D1, 1, rng)` with the java.util.Random built
     * beside it. Drawing this family off the MRG stream would both use the wrong
     * numbers and desynchronize every renewal draw that follows.
     *
     * The walk itself: hold the phase across calls (it is only advanced here,
     * so an idle server freezes its service process), draw the initial phase
     * from map_pie on the first call, then alternate an exponential holding
     * time at rate -D0(i,i) with a choice among the 2n competing transitions,
     * ending when the chosen one is a D1 (an event) rather than a D0 (a hidden
     * phase change). One uniform for the holding time and one for the choice,
     * per step, in that order.
     */
    double next_map_java(Rng& g) {
        const std::size_t n = map_.order();
        if (n == 1) {
            // The reference's exponential shortcut: -log(u)/lambda, with plain
            // log and not log1p, because Map_sample spells it that way.
            const double lambda = map_.D1(0, 0);
            return -std::log(g.aux.next_double()) / lambda;
        }
        if (!phase_known_) {
            const std::vector<double> pie = mam::map_pie(map_);
            double sum = 0.0;
            const double r = g.aux.next_double();
            phase_ = n - 1;
            for (std::size_t i = 0; i < n; ++i) {
                sum += pie[i];
                if (r < sum) {
                    phase_ = i;
                    break;
                }
            }
            phase_known_ = true;
        }
        double sample = 0.0;
        std::vector<double> row(2 * n, 0.0);
        bool go = true;
        while (go) {
            const double rate = -map_.D0(phase_, phase_);
            sample += -std::log(g.aux.next_double()) / rate;
            for (std::size_t k = 0; k < n; ++k) {
                row[k] = map_.D0(phase_, k);
                row[n + k] = map_.D1(phase_, k);
            }
            row[phase_] = 0.0;  // the diagonal is the exit rate, not a target
            std::size_t next_state = 2 * n - 1;
            double sum = 0.0;
            const double r = g.aux.next_double();
            for (std::size_t j = 0; j < 2 * n; ++j) {
                sum += row[j] / rate;
                if (r < sum) {
                    if (j >= n) {
                        next_state = j - n;
                        go = false;
                    } else {
                        next_state = j;
                        go = true;
                    }
                    break;
                }
            }
            if (next_state >= 2 * n - 1 && go) {
                // The loop above fell through: the reference leaves nextState at
                // its initialized 2n-1 and continues, which is a D1 transition.
                next_state = n - 1;
                go = false;
            }
            phase_ = (next_state < n) ? next_state : next_state - n;
        }
        return sample;
    }
 
    double next_markovian(Rng& g) {
        std::vector<double> out;
        if (type_ == lang::ProcessType::MAP || type_ == lang::ProcessType::MMPP2 ||
            type_ == lang::ProcessType::PH || type_ == lang::ProcessType::APH ||
            type_ == lang::ProcessType::COXIAN || type_ == lang::ProcessType::COX2 ||
            type_ == lang::ProcessType::ERLANG || type_ == lang::ProcessType::HYPEREXP) {
            return next_map_java(g);
        }
        if (type_ == lang::ProcessType::ME) {
            // An ME is a RENEWAL process: the entry law is the same at every
            // sample, so the inversion table built in `load_map` serves them
            // all. Routing it through `rap_sample` instead cost a matrix
            // exponential per bisection step and ran 140x slower than the Java
            // engine, which is a wall-clock timeout rather than a slow answer.
            return me_->next(g.mc);
        }
        if (type_ == lang::ProcessType::RAP) {
            // A RAP's state is the real-valued ENTRY LAW, not a discrete phase,
            // so it is chained through `a_out` rather than recovered from a
            // trace.
            std::vector<double> a_next;
            out = mam::rap_sample(map_, 1, g.mc, entry_, &a_next);
            entry_ = a_next;
        } else {
            mam::SampleTrace tr;
            std::vector<double> start;
            if (has_phase_ && phase_known_) {
                start.assign(map_.order(), 0.0);
                start[phase_] = 1.0;
            }
            out = mam::map_sample(map_, 1, g.mc, start, &tr);
            if (has_phase_ && !tr.last.empty()) {
                phase_ = tr.last[0];
                phase_known_ = true;
            }
        }
        if (out.empty()) return 0.0;
        // A discrete-time MAP samples a SLOT COUNT, which is already the
        // interval on the lattice; the continuous families sample the interval
        // directly. Neither needs rescaling here.
        return out[0];
    }
 
    template <class T>
    void load_schedule(const lang::Distrib<T>& d) {
        if (!d.has_schedule())
            throw InputError("SolverLDES (native engine): the time-inhomogeneous process at " +
                             where_ + " carries no segment schedule");
        for (const T& b : d.sched_bp) bp_.push_back(num_traits<T>::to_double(b));
        if (bp_.size() != d.sched_D0.size() + 1)
            throw InputError("SolverLDES (native engine): the schedule at " + where_ +
                             " has a boundary vector that does not bound its segments");
        for (std::size_t k = 0; k < d.sched_D0.size(); ++k) {
            const std::size_t H = d.sched_D0[k].rows();
            mam::Map<double> seg;
            seg.D0 = Matrix<double>(H, H, 0.0);
            seg.D1 = Matrix<double>(H, H, 0.0);
            for (std::size_t a = 0; a < H; ++a)
                for (std::size_t b2 = 0; b2 < H; ++b2) {
                    seg.D0(a, b2) = num_traits<T>::to_double(d.sched_D0[k](a, b2));
                    seg.D1(a, b2) = num_traits<T>::to_double(d.sched_D1[k](a, b2));
                }
            segs_.push_back(seg);
        }
        cyclic_ = d.sched_cyclic;
        has_schedule_ = true;
        phase_ = 0;
        if (!(mean_ > 0.0)) mean_ = num_traits<T>::to_double(d.mean);
        if (!(mean_ > 0.0) && !segs_.empty()) mean_ = num_traits<double>::to_double(
            mam::map_mean(segs_[0]));
    }
 
    double period() const { return bp_.back() - bp_.front(); }
 
    /** Index of the segment in force at t, or -1 past a non-cyclic horizon. */
    int segment_at(double t) const {
        double offset = t - bp_.front();
        const double per = period();
        if (cyclic_) {
            offset = std::fmod(offset, per);
            if (offset < 0.0) offset += per;
        } else if (offset < 0.0 || offset >= per) {
            return -1;
        }
        const double pos = bp_.front() + offset;
        for (std::size_t k = 0; k + 1 < bp_.size(); ++k)
            if (pos < bp_[k + 1]) return static_cast<int>(k);
        return static_cast<int>(segs_.size()) - 1;
    }
 
    /**
     * One interval of a MAPt / PHt / NHPP, started at `from`.
     *
     * The walk advances through the phase process and STOPS AT EVERY SEGMENT
     * BOUNDARY, resampling the holding time under the new generator. That is
     * what makes the schedule piecewise-constant rather than merely
     * time-averaged: an exponential holding time drawn under one segment's rate
     * and allowed to run past the boundary would carry the old rate into the
     * new segment, which is exactly the approximation the family exists to
     * avoid. Transcribes `sampleMAPtInterarrival`.
     */
    double schedule_sample(Rng& g, double from) {
        const std::size_t H = segs_[0].order();
        double elapsed = 0.0, pos = from;
        for (int guard = 0; guard < 1000000; ++guard) {
            const int idx = segment_at(pos);
            if (idx < 0) return 0.0;  // past a non-cyclic horizon: no more arrivals
            double offset = pos - bp_.front();
            if (cyclic_) {
                offset = std::fmod(offset, period());
                if (offset < 0.0) offset += period();
            }
            const double to_boundary = (bp_[static_cast<std::size_t>(idx) + 1] - bp_.front()) - offset;
            const mam::Map<double>& seg = segs_[static_cast<std::size_t>(idx)];
            const double total = -seg.D0(phase_, phase_);
            const bool last_segment =
                (!cyclic_ && static_cast<std::size_t>(idx) + 1 == segs_.size());
            if (!(total > 0.0)) {
                if (last_segment) return 0.0;
                elapsed += to_boundary;
                pos += to_boundary;
                continue;
            }
            const double holding = -std::log(uniform01(g)) / total;
            if (holding >= to_boundary) {
                if (last_segment) return 0.0;
                elapsed += to_boundary;
                pos += to_boundary;
                continue;
            }
            elapsed += holding;
            pos += holding;
            // The transitions competing out of the current phase, ARRIVALS
            // FIRST: a D1 entry ends the interval, a D0 entry only moves phase.
            const double u = uniform01(g) * total;
            double cum = 0.0;
            int chosen = -1;
            for (std::size_t j = 0; j < 2 * H; ++j) {
                const double w = (j < H) ? seg.D1(phase_, j)
                                         : ((j - H == phase_) ? 0.0 : seg.D0(phase_, j - H));
                cum += w;
                if (u < cum) {
                    chosen = static_cast<int>(j);
                    break;
                }
            }
            if (chosen < 0)
                for (std::size_t j = 2 * H; j-- > 0;) {
                    const double w = (j < H) ? seg.D1(phase_, j)
                                             : ((j - H == phase_) ? 0.0 : seg.D0(phase_, j - H));
                    if (w > 0.0) {
                        chosen = static_cast<int>(j);
                        break;
                    }
                }
            if (chosen < static_cast<int>(H)) {
                phase_ = static_cast<std::size_t>(chosen);
                return elapsed;
            }
            phase_ = static_cast<std::size_t>(chosen) - H;
        }
        return 0.0;
    }
 
    lang::ProcessType type_ = lang::ProcessType::DISABLED;
    std::string where_;
    double mean_ = 0.0, scv_ = 1.0, rate_ = 0.0;
    double a_ = 0.0, b_ = 0.0;  ///< the family's two shape/scale slots
    int erlang_k_ = 0;          ///< phases of an Erlang, 0 when not recovered
    double erlang_rate_ = 0.0;  ///< per-phase rate of an Erlang
    std::vector<double> trace_;
    std::size_t trace_idx_ = 0;
    mam::Map<double> map_;
    std::vector<double> entry_;  ///< RAP entry law, chained across samples
    std::shared_ptr<mam::MeSampler<double>> me_;  ///< ME inversion table, built once
    bool has_phase_ = false, phase_known_ = false;
    // The time-inhomogeneous schedule: boundaries, per-segment (D0,D1), and
    // whether the schedule repeats past its last boundary.
    bool has_schedule_ = false, cyclic_ = false;
    std::vector<double> bp_;
    std::vector<mam::Map<double>> segs_;
    std::size_t phase_ = 0;
};
 
}  // namespace engine
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_SAMPLER_H