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(LdesStreamSeeds.mrg32k3a(seed, offset));
 *     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. The MRG
 * stream's six seeds are HASHED from (seed, offset) (`rng::SeedHash`, since
 * equally spaced seed vectors give linearly dependent streams); the
 * java.util.Random keeps `seed + offset`, whose seed scrambling and 48-bit
 * state show no such dependence (u1 - 2*u2 + u3 over three spaced seeds is
 * uniform to within its standard error).
 *
 * `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::DUNIFORM: {
                // The reference carries a DiscreteUniform as the moment pair
                // {mean, SCV}, not as its bounds, so the bounds are recovered from
                // the pair HERE TOO rather than read off `d.params`: the two
                // engines must place the lattice identically, and the Java one has
                // nothing but the moments. var = (w^2-1)/12 gives the integral
                // width, and the mean then places it, which leaves a lower bound
                // off the integers where the user put it.
                require_moments();
                const double var = scv_ * mean_ * mean_;
                b_ = std::floor(std::sqrt(12.0 * var + 1.0) + 0.5);  // width
                a_ = mean_ - (b_ - 1.0) / 2.0;                       // lower bound
                if (!(b_ >= 1.0))
                    throw InputError("SolverLDES (native engine): the DiscreteUniform at " +
                                     where_ + " has moments no (min,max) pair realises");
                if (!(a_ >= 0.0))
                    throw InputError("SolverLDES (native engine): the DiscreteUniform at " +
                                     where_ + " reaches below zero");
                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:
            case lang::ProcessType::MMAPT:
            case lang::ProcessType::MPHT:
            case lang::ProcessType::BMMAPT:
                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::MPH:
            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) { return resolve_zero_atom(draw(g)); }
 
    /**
     * The reference's `resolveZeroAtom`: the discrete laws supported on
     * {0,1,...} can return a zero-length interval, and a zero-length service is
     * an infinite rate no estimator carries. `Immediate` is exactly what a
     * zero-length interval denotes elsewhere in the engine, so the atom is
     * served at that scale. Under SLOTTED mode 1/Immediate is not a lattice
     * point, so `slot_snap` refuses the sample rather than inventing one, which
     * is the refusal the reference raises in its own words.
     */
    double resolve_zero_atom(double v) const {
        if (v != 0.0) return v;
        switch (type_) {
            case lang::ProcessType::BERNOULLI:
            case lang::ProcessType::BINOMIAL:
            case lang::ProcessType::POISSON:
            case lang::ProcessType::DUNIFORM:
                return 1.0 / lang::GlobalConstants::Immediate;
            default:
                return v;
        }
    }
 
private:
    double draw(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::DUNIFORM:
                return ssj::uniform_int_inverse(a_, b_, uniform01(g));
            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);
        }
    }
 
public:
    /**
     * 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; }
    /**
     * The 1-based mark of the interval last returned by `next`/`next_at`, or 0 when
     * this process carries none. At a Source the engine reads it to pick the class of
     * the arriving job, exactly as `Solver_ssj` consumes `arrivalPendingMark`.
     */
    int last_mark() const { return last_mark_; }
    /**
     * How many jobs the epoch last returned by `next`/`next_at` releases, or 1 when
     * this process carries no batch axis. At a Source the engine reads it to release
     * that many jobs at one instant, exactly as `Solver_ssj` consumes
     * `arrivalBatchSize`; at a bulk server it is how many completions the firing
     * clears.
     */
    int last_batch() const { return last_batch_ > 0 ? last_batch_ : 1; }
    bool marked() const { return !mark_.empty() || !seg_mark_.empty(); }
    /** True when the process carries its OWN batch sizes: a BMAP or a BMMAPt. */
    bool batched() const { return !batch_.empty() || !seg_batch_.empty(); }
    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");
        }
        // THE PER-MARK BLOCKS. `d.D1` is the aggregate sum_k D1k, which is what the
        // interval walk uses; `Dmark` tells the marks apart once a D1 transition has
        // fired. An MPH arrives here already lowered (D0 = S, D1k = s_k*alpha), so
        // the same blocks serve it and no separate renewal path is needed.
        mark_.clear();
        batch_.clear();
        if (lang::process_is_marked_stationary(type_)) {
            for (const Matrix<T>& Dk : d.Dmark) {
                Matrix<double> m(K, K, 0.0);
                for (std::size_t i = 0; i < K; ++i)
                    for (std::size_t j = 0; j < K; ++j)
                        m(i, j) = num_traits<T>::to_double(Dk(i, j));
                mark_.push_back(m);
            }
        }
        // A BMAP OVERLOADS THE SAME FIELD WITH BATCH SIZES: `Dmark[b]` is the block
        // that releases b+1 jobs. They used to be read into `mark_` and never
        // consulted, because `mark_` is filled only for a marked stationary process
        // and nothing else looked -- so this engine answered a BMAP model as its
        // AGGREGATE MAP, every batch collapsed to one job, in silence. The blocks now
        // drive `next_bmap_java`, which is the transcription of the reference's
        // `Map_sample.BmapSampler`.
        if (type_ == lang::ProcessType::BMAP) {
            for (const Matrix<T>& Db : d.Dmark) {
                Matrix<double> m(K, K, 0.0);
                for (std::size_t i = 0; i < K; ++i)
                    for (std::size_t j = 0; j < K; ++j)
                        m(i, j) = num_traits<T>::to_double(Db(i, j));
                batch_.push_back(m);
            }
        }
        // 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 ||
                      lang::process_is_marked_stationary(type_) ||
                      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;
    }
 
    /**
     * The 1-based mark of an arrival firing on transition (i, j), drawn with
     * probability D1k(i,j)/D1_agg(i,j). Transcribes `Mmap_sample.MmapSampler.sampleMark`,
     * including its fallbacks: an unmarked process reports mark 1, and so does a
     * transition on which every block is zero.
     */
    int sample_mark(std::size_t i, std::size_t j, Rng& g) {
        const std::size_t C = mark_.size();
        if (C == 0) return 1;
        double total = 0.0;
        for (std::size_t c = 0; c < C; ++c) total += std::max(0.0, mark_[c](i, j));
        if (!(total > 0.0)) return 1;
        const double r = g.aux.next_double() * total;
        double sum = 0.0;
        for (std::size_t c = 0; c < C; ++c) {
            sum += std::max(0.0, mark_[c](i, j));
            if (r < sum) return static_cast<int>(c) + 1;
        }
        return static_cast<int>(C);
    }
 
    /**
     * One marked interval, transcribed from `Mmap_sample.MmapSampler.next`.
     *
     * The interval walk is the unmarked one over the aggregate pair; the mark is a
     * SECOND draw taken on the transition that ended it, so an unmarked and a marked
     * process of the same pair do not share a sample path past the first arrival.
     * That is the reference's own ordering and must not be reordered to save a draw.
     *
     * Unlike `next_map_java` this does NOT carry the reference's fall-through repair:
     * `MmapSampler` leaves nextState at its initialized 2n-1 and reads it as a D1
     * transition into phase n-1, which is what is transcribed here.
     */
    double next_mmap_java(Rng& g) {
        const std::size_t n = map_.order();
        if (n == 1) {
            const double lambda = map_.D1(0, 0);
            const double time = -std::log(g.aux.next_double()) / lambda;
            last_mark_ = sample_mark(0, 0, g);
            return time;
        }
        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);
        while (true) {
            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) {
                    next_state = j;
                    break;
                }
            }
            if (next_state >= n) {
                const std::size_t dest = next_state - n;
                last_mark_ = sample_mark(phase_, dest, g);
                phase_ = dest;
                return sample;
            }
            phase_ = next_state;
        }
    }
 
    /**
     * One BATCH interval, transcribed from `Map_sample.BmapSampler.next`.
     *
     * Unlike the marked walk, the batch label is NOT a second draw: the block that
     * fires decides both the destination phase and how many jobs the epoch releases,
     * so the row scanned is D0 followed by the batch blocks BATCH-MAJOR then
     * destination, which is the reference's own `row[nphases + b*nphases + k]`
     * indexing. Reordering it would keep the same law and lose the sample path.
     *
     * The one-phase case is the reference's fast path and takes TWO draws: the
     * exponential interval, then the batch size against the per-batch rates.
     */
    double next_bmap_java(Rng& g) {
        const std::size_t n = map_.order();
        const std::size_t B = batch_.size();
        last_batch_ = 1;
        if (n == 1) {
            const double lambda = map_.D1(0, 0);
            const double time = -std::log(g.aux.next_double()) / lambda;
            double total = 0.0;
            for (std::size_t b = 0; b < B; ++b) total += batch_[b](0, 0);
            if (total > 0.0) {
                const double r = g.aux.next_double() * total;
                double cum = 0.0;
                for (std::size_t b = 0; b < B; ++b) {
                    cum += batch_[b](0, 0);
                    if (r < cum) {
                        last_batch_ = static_cast<int>(b) + 1;
                        break;
                    }
                }
            }
            return time;
        }
        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;
        }
        const std::size_t cols = n + B * n;
        std::vector<double> row(cols, 0.0);
        double sample = 0.0;
        while (true) {
            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[phase_] = 0.0;  // the diagonal is the exit rate, not a target
            for (std::size_t b = 0; b < B; ++b)
                for (std::size_t k = 0; k < n; ++k) row[n + b * n + k] = batch_[b](phase_, k);
            std::size_t next_state = phase_;
            bool fired = false;
            double sum = 0.0;
            const double r = g.aux.next_double();
            for (std::size_t j = 0; j < cols; ++j) {
                sum += row[j] / rate;
                if (r < sum) {
                    if (j < n) {
                        next_state = j;
                    } else {
                        const std::size_t idx = j - n;
                        last_batch_ = static_cast<int>(idx / n) + 1;
                        next_state = idx % n;
                        fired = true;
                    }
                    break;
                }
            }
            phase_ = next_state;
            if (fired) 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 (lang::process_is_marked_stationary(type_)) {
            return next_mmap_java(g);
        }
        if (type_ == lang::ProcessType::BMAP && !batch_.empty()) {
            return next_bmap_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);
        }
        // The per-mark blocks, segment-major here (`seg_mark_[k][c]`) because the walk
        // scans one segment at a time; the wire and the Distrib are mark-major.
        seg_mark_.clear();
        if (d.has_marked_schedule()) {
            const std::size_t C = d.sched_Dmark.size();
            for (std::size_t k = 0; k < d.sched_D0.size(); ++k) {
                const std::size_t H = d.sched_D0[k].rows();
                std::vector<Matrix<double>> per_mark;
                for (std::size_t c = 0; c < C; ++c) {
                    Matrix<double> m(H, H, 0.0);
                    for (std::size_t a = 0; a < H; ++a)
                        for (std::size_t b2 = 0; b2 < H; ++b2)
                            m(a, b2) = num_traits<T>::to_double(d.sched_Dmark[c][k](a, b2));
                    per_mark.push_back(m);
                }
                seg_mark_.push_back(per_mark);
            }
        }
        // The batch blocks, transposed the same way: the Distrib is mark-major then
        // batch then segment, the walk is segment-major, so this is
        // `seg_batch_[k][c][b]`. `seg_mark_` stays the batch-AGGREGATED per-mark
        // schedule and `segs_[k].D1` the aggregate over marks as well, so a consumer
        // that ignores batches walks exactly the MMAPt this hides down to.
        seg_batch_.clear();
        if (d.has_batch_schedule()) {
            const std::size_t C = d.sched_Dbatch.size();
            const std::size_t B = d.sched_Dbatch[0].size();
            for (std::size_t k = 0; k < d.sched_D0.size(); ++k) {
                const std::size_t H = d.sched_D0[k].rows();
                std::vector<std::vector<Matrix<double>>> per_mark;
                for (std::size_t c = 0; c < C; ++c) {
                    std::vector<Matrix<double>> per_batch;
                    for (std::size_t b = 0; b < B; ++b) {
                        Matrix<double> m(H, H, 0.0);
                        for (std::size_t a = 0; a < H; ++a)
                            for (std::size_t b2 = 0; b2 < H; ++b2)
                                m(a, b2) = num_traits<T>::to_double(d.sched_Dbatch[c][b][k](a, b2));
                        per_batch.push_back(m);
                    }
                    per_mark.push_back(per_batch);
                }
                seg_batch_.push_back(per_mark);
            }
        }
        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(); }
 
    /**
     * Position of t within one period, measured from bp_.front(), or -1 past a
     * non-cyclic horizon.
     *
     * The walk below carries this REDUCED offset rather than the wall clock: it
     * stays bounded by the period however far the simulation has run, which is
     * what keeps a boundary crossing representable.
     */
    double offset_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.0;
        }
        return offset;
    }
 
    /** Index of the segment holding an offset already reduced to [0, period). */
    int segment_of(double offset) const {
        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;
    }
 
    /** Index of the segment in force at t, or -1 past a non-cyclic horizon. */
    int segment_at(double t) const {
        const double offset = offset_at(t);
        return (offset < 0.0) ? -1 : segment_of(offset);
    }
 
    /**
     * 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 offset = offset_at(from);
        if (offset < 0.0) return 0.0;  // past a non-cyclic horizon: no more events
        std::size_t idx = static_cast<std::size_t>(segment_of(offset));
        double elapsed = 0.0;
        for (int guard = 0; guard < 1000000; ++guard) {
            const double to_boundary = (bp_[idx + 1] - bp_.front()) - offset;
            const mam::Map<double>& seg = segs_[idx];
            const double total = -seg.D0(phase_, phase_);
            // An absorbing phase never fires, so it always reaches the boundary.
            const double holding = (total > 0.0) ? -std::log(uniform01(g)) / total
                                                 : std::numeric_limits<double>::infinity();
            if (holding >= to_boundary) {
                ++idx;
                if (idx >= segs_.size()) {
                    if (!cyclic_) return 0.0;
                    idx = 0;
                }
                elapsed += to_boundary;
                // LAND EXACTLY ON THE BOUNDARY by taking the next segment's own
                // offset. Adding to_boundary to a running position instead would
                // stall the walk: far from the schedule origin the sliver left by
                // the modulo is below one ulp of that position, so the addition
                // does not move it and the loop spins on the same boundary until
                // the guard trips -- a hang, not a wrong number, and one that only
                // shows up once the simulation clock has run far enough.
                offset = bp_[idx] - bp_.front();
                continue;
            }
            elapsed += holding;
            offset += 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;
            // A MARKED schedule splits the arrival arm across the mark blocks at NO EXTRA
            // DRAW: the scan runs destination-major and mark-minor, so the running total
            // after all marks of a destination equals the aggregate total after that
            // destination. The winning destination is therefore the one the unmarked walk
            // would pick for the same uniform, which is what makes a K = 1 MMAPt reproduce
            // the MAPt with the same matrices sample path for sample path.
            if (!seg_mark_.empty()) {
                const std::vector<Matrix<double>>& marks = seg_mark_[idx];
                // The BATCH axis goes INSIDE the mark axis, so dropping it (every
                // batch size 1) leaves the scan order untouched and a B = 1 BMMAPt
                // reproduces the MMAPt sample path exactly as a K = 1 MMAPt
                // reproduces the MAPt. Neither label costs an extra draw.
                const bool batched_here = !seg_batch_.empty();
                const std::size_t B = batched_here ? seg_batch_[idx][0].size() : std::size_t(1);
                for (std::size_t j = 0; j < H; ++j) {
                    for (std::size_t c = 0; c < marks.size(); ++c) {
                        if (!batched_here) {
                            cum += marks[c](phase_, j);
                            if (u < cum) {
                                last_mark_ = static_cast<int>(c) + 1;
                                last_batch_ = 1;
                                phase_ = j;
                                return elapsed;
                            }
                            continue;
                        }
                        for (std::size_t b = 0; b < B; ++b) {
                            cum += seg_batch_[idx][c][b](phase_, j);
                            if (u < cum) {
                                last_mark_ = static_cast<int>(c) + 1;
                                last_batch_ = static_cast<int>(b) + 1;
                                phase_ = j;
                                return elapsed;
                            }
                        }
                    }
                }
                for (std::size_t j = 0; j < H; ++j) {
                    if (j == phase_) continue;
                    cum += seg.D0(phase_, j);
                    if (u < cum) {
                        chosen = static_cast<int>(H + j);
                        break;
                    }
                }
                if (chosen < 0) {
                    // Rounding shortfall: attribute the draw to the last positive arrival
                    // entry, as the unmarked walk below does.
                    for (std::size_t j = H; j-- > 0;)
                        for (std::size_t c = marks.size(); c-- > 0;) {
                            if (!batched_here) {
                                if (marks[c](phase_, j) > 0.0) {
                                    last_mark_ = static_cast<int>(c) + 1;
                                    last_batch_ = 1;
                                    phase_ = j;
                                    return elapsed;
                                }
                                continue;
                            }
                            for (std::size_t b = B; b-- > 0;)
                                if (seg_batch_[idx][c][b](phase_, j) > 0.0) {
                                    last_mark_ = static_cast<int>(c) + 1;
                                    last_batch_ = static_cast<int>(b) + 1;
                                    phase_ = j;
                                    return elapsed;
                                }
                        }
                    last_mark_ = 1;
                    last_batch_ = 1;
                    return elapsed;
                }
                phase_ = static_cast<std::size_t>(chosen) - H;
                continue;
            }
            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
    /// Per-SEGMENT, per-mark arrival blocks of a marked schedule; empty otherwise.
    std::vector<std::vector<Matrix<double>>> seg_mark_;
    /// Per-mark arrival blocks D1k; empty for an unmarked process.
    std::vector<Matrix<double>> mark_;
    /// Per-BATCH-SIZE blocks of a stationary BMAP, `batch_[b]` releasing b+1 jobs;
    /// empty for every other family. The schedule twin is `seg_batch_`.
    std::vector<Matrix<double>> batch_;
    /// Per-SEGMENT, per-mark, per-batch blocks of a BMMAPt; empty otherwise.
    std::vector<std::vector<std::vector<Matrix<double>>>> seg_batch_;
    /// The 1-based mark of the interval last returned; 0 before the first draw.
    int last_mark_ = 0;
    /// How many jobs the epoch last returned releases; 0 before the first draw.
    int last_batch_ = 0;
    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