ldes_stats.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_STATS_H
#define LINE_SOLVERS_LDES_LDES_STATS_H
 
/**
 * @file
 * @ingroup line_solvers
 * The estimators of the native LDES engine: running integrals, the MSER-5
 * warmup filter, and the batch-means confidence intervals.
 *
 * THE INTEGRALS ARE LAZY. `tot_qlen(i,r)` is the time integral of the class-r
 * queue length at station i, advanced only when that pair CHANGES. Advancing
 * every pair on every event would be O(M*K) per event on a hot loop where the
 * event touches one pair; the reference does the same thing with
 * `lastQueueUpdateTime`, and the invariant is that `update_qlen` is called
 * BEFORE any write to `qlen`, never after.
 *
 * THE WARMUP FILTER IS MSER-5 over EVENT-SPACED observations, not time-spaced
 * ones. That matters for what the estimators may then do with the series: an
 * unweighted mean of per-interval averages overweights congested epochs,
 * because a congested epoch contains more events and therefore more intervals
 * per unit time. So the means below are integral DIFFERENCES over elapsed
 * time, and only the CI -- which needs a sequence of comparable observations,
 * not a mean -- reads the per-interval series directly.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/sim/sim_dist.h"
#include "line/solvers/wrappers/ldes/ldes_options.h"
#include "line/util/matrix.h"
 
namespace line {
namespace ldes {
namespace engine {
 
/** The running per-(station, class) integrals and tallies. */
struct Accum {
    Accum(std::size_t M, std::size_t K)
        : qlen(M, std::vector<double>(K, 0.0)),
          held(M, std::vector<double>(K, 0.0)),
          busy(M, std::vector<double>(K, 0.0)),
          tot_qlen(M, std::vector<double>(K, 0.0)),
          tot_busy(M, std::vector<double>(K, 0.0)),
          last_qlen(M, std::vector<double>(K, 0.0)),
          last_busy(M, std::vector<double>(K, 0.0)),
          busy_scale(M, 1.0),
          util_peak(M, 1.0),
          resp_sum(M, std::vector<double>(K, 0.0)),
          resp_cnt(M, std::vector<double>(K, 0.0)),
          resp_sq(M, std::vector<double>(K, 0.0)),
          completed(M, std::vector<double>(K, 0.0)),
          arrived(M, std::vector<double>(K, 0.0)),
          join_dropped(M, std::vector<double>(K, 0.0)) {}
 
    /** Advance the queue-length integral of (i,r) to `now`. Call BEFORE writing qlen. */
    void update_qlen(std::size_t i, std::size_t r, double now) {
        const double dt = now - last_qlen[i][r];
        if (dt > 0.0) tot_qlen[i][r] += (qlen[i][r] + held[i][r]) * dt;
        last_qlen[i][r] = now;
    }
    /**
     * Advance the busy-server integral of (i,r) to `now`. Call BEFORE writing busy.
     *
     * The integrand is the WORK being delivered, `busy * busy_scale`, not the
     * head count: a load-dependent server running alpha(n) times faster does the
     * same work in less time, so plain busy TIME reads it as no busier than one
     * at its nominal rate. `busy_scale` is 1 unless the station is load
     * dependent, so this is the identity everywhere else.
     */
    void update_busy(std::size_t i, std::size_t r, double now) {
        const double dt = now - last_busy[i][r];
        if (dt > 0.0) tot_busy[i][r] += busy[i][r] * busy_scale[i] * dt;
        last_busy[i][r] = now;
    }
    /**
     * Install the load-dependent speed station `i` runs at from `now` on.
     *
     * Flushes every class FIRST, so the interval that just ended is credited at
     * the scale that was actually in force over it. The engine calls this after
     * the population is written, from the same chokepoint that re-times the
     * departures.
     */
    void set_busy_scale(std::size_t i, double v, double now) {
        if (v == busy_scale[i]) return;
        for (std::size_t r = 0; r < busy[i].size(); ++r) update_busy(i, r, now);
        busy_scale[i] = v;
    }
 
    std::vector<std::vector<double>> qlen;
    /**
     * Callers parked on their server awaiting a synchronous REPLY (sn.syncreply), per (station, class).
     *
     * They are counted in the QLen INTEGRAL, as the reference counts currentBlockedServers in
     * Solver_ssj.updateQueueStats, but kept OUT of `qlen`, which every decision reads: capacity,
     * load-dependent rates, JSQ and setup all use currentQueueLength there, which excludes them.
     * Written only after update_qlen, like qlen.
     */
    std::vector<std::vector<double>> held;
    std::vector<std::vector<double>> busy;
    std::vector<std::vector<double>> tot_qlen, tot_busy;
    std::vector<std::vector<double>> last_qlen, last_busy;
    /**
     * The speed each station is running at right now (load dependence only),
     * and the peak capacity that normalizes its utilization,
     * max(c, max(alpha)).
     *
     * `util_peak` is CTMC's own `ceff`, which is what makes LDES report the same
     * utilization as the analytic solvers on a load-dependent station. Both are
     * 1 and c respectively when the station is not load dependent, leaving every
     * other model's numbers untouched.
     */
    std::vector<double> busy_scale, util_peak;
    std::vector<std::vector<double>> resp_sum, resp_cnt, completed;
    /** Sum of squared response times, for the iid interval the reference takes off its response-time Tally. */
    std::vector<std::vector<double>> resp_sq;
    /** Tally one response time `x` at (i,r): the sum, its square and the count move together. */
    void tally_resp(std::size_t i, std::size_t r, double x) {
        resp_sum[i][r] += x;
        resp_sq[i][r] += x * x;
        resp_cnt[i][r] += 1.0;
    }
    /**
     * Jobs that ARRIVED at each (station, class), which is not what `completed`
     * counts: the reference reports AN as the arrivals after the truncation point
     * over the truncated window (cumulative, differenced like `completed`), and the
     * two differ by whatever is still in the station at the end and by anything
     * dropped, balked or reneged. Reporting completions as the arrival rate
     * would hide exactly the loss those features exist to measure.
     */
    std::vector<std::vector<double>> arrived;
    /**
     * Siblings a QUORUM Join discarded, per (station, class): a sibling that
     * reaches the Join after its parent already fired lost the race and is
     * thrown away. It is a completion of nothing, so it belongs neither in
     * `completed` nor in `arrived`, and it is what a Join row's loss rate is.
     * Cumulative, and differenced against the truncation point like `completed`.
     */
    std::vector<std::vector<double>> join_dropped;
};
 
/**
 * MSER-5 truncation point over a series of observations, in BATCHES.
 *
 * Returns the batch index d minimising var(Z_{d..N-1}) / (N-d)^2 over the
 * batch means Z, and 0 when there are fewer than four batches -- below that
 * the criterion is estimated from too few terms to be a truncation rule.
 * Transcribes `computeMSER5TruncationPoint`.
 */
inline int mser5_truncation(const std::vector<double>& obs, int batch_size) {
    const std::size_t n = obs.size();
    if (batch_size <= 0 || n < static_cast<std::size_t>(batch_size) * 4) return 0;
    const int num_batches = static_cast<int>(n / static_cast<std::size_t>(batch_size));
    std::vector<double> means(static_cast<std::size_t>(num_batches), 0.0);
    for (int j = 0; j < num_batches; ++j) {
        double sum = 0.0;
        for (int i = 0; i < batch_size; ++i)
            sum += obs[static_cast<std::size_t>(j * batch_size + i)];
        means[static_cast<std::size_t>(j)] = sum / batch_size;
    }
    double min_mser = std::numeric_limits<double>::max();
    int optimal_d = 0;
    const int max_d = num_batches / 2;
    for (int d = 0; d < max_d; ++d) {
        const int remaining = num_batches - d;
        if (remaining < 2) break;
        double sum = 0.0;
        for (int j = d; j < num_batches; ++j) sum += means[static_cast<std::size_t>(j)];
        const double mean = sum / remaining;
        double variance = 0.0;
        for (int j = d; j < num_batches; ++j) {
            const double diff = means[static_cast<std::size_t>(j)] - mean;
            variance += diff * diff;
        }
        variance /= (remaining - 1);
        const double mser = variance / (static_cast<double>(remaining) * remaining);
        if (mser < min_mser) {
            min_mser = mser;
            optimal_d = d;
        }
    }
    return optimal_d;
}
 
/** Where the warmup ended, and whether a truncation was applied at all. */
struct Truncation {
    std::size_t index = 0;    ///< observation index of the truncation point
    double warmup_end = 0.0;  ///< the instant of that observation
    bool applied = false;
};
 
/** The event-spaced observation series MSER-5 and the CI both read. */
struct Observations {
    Observations(std::size_t M, std::size_t K, bool mser_on, int batch)
        : enabled(mser_on),
          batch_size(batch),
          qlen(M, std::vector<std::vector<double>>(K)),
          qt(M, std::vector<std::vector<double>>(K)),
          bt(M, std::vector<std::vector<double>>(K)),
          cmp(M, std::vector<std::vector<double>>(K)),
          drp(M, std::vector<std::vector<double>>(K)),
          arv(M, std::vector<std::vector<double>>(K)),
          last_qt(M, std::vector<double>(K, 0.0)),
          in_mser(M, static_cast<char>(1)),
          nstations(M),
          nclasses(K) {}
 
    /**
     * Record one observation.
     *
     * The queue-length entry is the INTERVAL TIME-AVERAGE, not the instantaneous
     * count: it is the integral delivered since the previous observation over
     * the interval length, which is what makes a sequence of them comparable
     * even though the intervals differ in duration. The other three are
     * CUMULATIVE, because the estimators difference them against the truncation
     * point rather than averaging them.
     */
    void collect(const Accum& acc, double now) {
        const double dt = now - last_time;
        time.push_back(now);
        for (std::size_t i = 0; i < nstations; ++i)
            for (std::size_t r = 0; r < nclasses; ++r) {
                if (dt > 0.0) {
                    qlen[i][r].push_back((acc.tot_qlen[i][r] - last_qt[i][r]) / dt);
                    last_qt[i][r] = acc.tot_qlen[i][r];
                } else {
                    qlen[i][r].push_back(acc.qlen[i][r]);
                }
                qt[i][r].push_back(acc.tot_qlen[i][r]);
                bt[i][r].push_back(acc.tot_busy[i][r]);
                cmp[i][r].push_back(acc.completed[i][r]);
                drp[i][r].push_back(acc.join_dropped[i][r]);
                arv[i][r].push_back(acc.arrived[i][r]);
            }
        last_time = now;
    }
 
    /**
     * The truncation point, on the AGGREGATE queue length first.
     *
     * The per-series fallback is not a refinement, it is the closed-network
     * case: a closed model holds a constant total population, so the aggregate
     * series is flat and carries no transient signal at all, however long the
     * warmup actually was. The most conservative per-series point is taken then.
     *
     * BOTH RUN OVER THE SERVICE STATIONS ONLY, which is the range
     * `applyMSER5Truncation` walks. A Join is measured like a station but is
     * not one of them, and folding its series in would move the truncation
     * point of every fork-join model away from the reference's.
     */
    Truncation truncate() const {
        Truncation t;
        if (!enabled || time.empty()) return t;
        std::vector<double> aggregate(time.size(), 0.0);
        for (std::size_t k = 0; k < time.size(); ++k) {
            double tot = 0.0;
            for (std::size_t i = 0; i < nstations; ++i) {
                if (!in_mser[i]) continue;
                for (std::size_t r = 0; r < nclasses; ++r)
                    if (k < qlen[i][r].size()) tot += qlen[i][r][k];
            }
            aggregate[k] = tot;
        }
        int batch = mser5_truncation(aggregate, batch_size);
        if (batch == 0) {
            for (std::size_t i = 0; i < nstations; ++i) {
                if (!in_mser[i]) continue;
                for (std::size_t r = 0; r < nclasses; ++r) {
                    const int b = mser5_truncation(qlen[i][r], batch_size);
                    if (b > batch) batch = b;
                }
            }
        }
        const std::size_t idx =
            static_cast<std::size_t>(batch) * static_cast<std::size_t>(batch_size);
        if (idx < time.size()) {
            t.index = idx;
            t.warmup_end = time[idx];
            t.applied = true;
        }
        return t;
    }
 
    bool enabled = true;
    int batch_size = 5;
    std::vector<double> time;
    std::vector<std::vector<std::vector<double>>> qlen, qt, bt, cmp, drp, arv;
    std::vector<std::vector<double>> last_qt;
    /// Whether station i feeds the truncation criterion; every one is still recorded.
    std::vector<char> in_mser;
    double last_time = 0.0;
    std::size_t nstations = 0, nclasses = 0;
};
 
/** Grand mean, standard error and degrees of freedom of a batch-means estimate. */
struct StatTriple {
    double mean = 0.0;
    double stderr_ = 0.0;
    int df = 0;
    bool ok = false;
};
 
/**
 * The variance inflation of OVERLAPPING batch means.
 *
 * Overlapping batches are positively correlated, so the naive sample variance
 * of their means underestimates the true one; 4/3 is the standard asymptotic
 * correction at 50% overlap, interpolated linearly below it and reducing to 1
 * at zero overlap, where the batches are disjoint and no correction is due.
 * Transcribes `computeOverlapAdjustmentFactor`.
 */
inline double overlap_adjustment(double overlap) {
    if (overlap <= 0.0) return 1.0;
    if (overlap >= 0.5) return 4.0 / 3.0;
    return 1.0 + (overlap / 0.5) * (4.0 / 3.0 - 1.0);
}
 
/** Non-overlapping batch means. Transcribes `computeBMStatisticsInternal`. */
inline StatTriple bm_statistics(const std::vector<double>& obs, int batch_size) {
    StatTriple s;
    const std::size_t n = obs.size();
    const int nb = (batch_size > 0) ? static_cast<int>(n / static_cast<std::size_t>(batch_size)) : 0;
    if (nb < 2) return s;
    std::vector<double> means(static_cast<std::size_t>(nb), 0.0);
    for (int i = 0; i < nb; ++i) {
        double sum = 0.0;
        for (int j = 0; j < batch_size; ++j)
            sum += obs[static_cast<std::size_t>(i * batch_size + j)];
        means[static_cast<std::size_t>(i)] = sum / batch_size;
    }
    double grand = 0.0;
    for (int i = 0; i < nb; ++i) grand += means[static_cast<std::size_t>(i)];
    grand /= nb;
    double ss = 0.0;
    for (int i = 0; i < nb; ++i) {
        const double d = means[static_cast<std::size_t>(i)] - grand;
        ss += d * d;
    }
    s.mean = grand;
    s.stderr_ = std::sqrt((ss / (nb - 1)) / nb);
    s.df = nb - 1;
    s.ok = true;
    return s;
}
 
/** Overlapping batch means. Transcribes `computeOBMStatisticsInternal`. */
inline StatTriple obm_statistics(const std::vector<double>& obs, int batch_size, double overlap) {
    StatTriple s;
    const std::size_t n = obs.size();
    if (batch_size <= 0 || n < static_cast<std::size_t>(batch_size) * 2) return s;
    int step = static_cast<int>(batch_size * (1.0 - overlap));
    if (step < 1) step = 1;
    const int nb = static_cast<int>((n - static_cast<std::size_t>(batch_size)) /
                                    static_cast<std::size_t>(step)) + 1;
    if (nb < 2) return s;
    std::vector<double> means(static_cast<std::size_t>(nb), 0.0);
    for (int i = 0; i < nb; ++i) {
        const std::size_t start = static_cast<std::size_t>(i) * static_cast<std::size_t>(step);
        double sum = 0.0;
        for (int j = 0; j < batch_size; ++j)
            if (start + static_cast<std::size_t>(j) < n) sum += obs[start + static_cast<std::size_t>(j)];
        means[static_cast<std::size_t>(i)] = sum / batch_size;
    }
    double grand = 0.0;
    for (int i = 0; i < nb; ++i) grand += means[static_cast<std::size_t>(i)];
    grand /= nb;
    double ss = 0.0;
    for (int i = 0; i < nb; ++i) {
        const double d = means[static_cast<std::size_t>(i)] - grand;
        ss += d * d;
    }
    const double adj = overlap_adjustment(overlap);
    s.mean = grand;
    s.stderr_ = std::sqrt((adj * ss / (nb - 1)) / nb);
    s.df = static_cast<int>((nb - 1) / adj);
    if (s.df < 1) s.df = 1;
    s.ok = true;
    return s;
}
 
/**
 * The two-sided t critical value at confidence `level` with `df` degrees of freedom: the EXACT Student-t
 * quantile t_{(1+level)/2, df}, as the Java engine's `Solver_ssj.getTCriticalValue` computes it.
 *
 * Both engines used a three-decimal table, bucketed to the 90/95/99% levels and saturating at df = 30, until
 * 2026-10-28; that moved every half-width by parts in a thousand and kept df > 30 from tending to the normal
 * quantile. Below one degree of freedom there is no t law, and the normal quantile at the same level stands in.
 */
inline double t_critical(double level, int df) {
    const double p = 0.5 + 0.5 * level;
    if (df < 1) return sim::sim_norminv(p);
    return sim::sim_tinv(p, static_cast<double>(df));
}
 
 
/**
 * Heidelberger-Welch spectral estimate of the variance of the sample mean.
 *
 * The variance of a mean over a correlated series is 2*pi*S(0)/m, with S(0)
 * the spectral density at frequency zero. S(0) cannot be read off the
 * periodogram directly -- the ordinate at f=0 is exactly the quantity being
 * estimated and carries all the bias -- so the LOG-periodogram over the lowest
 * frequencies is fitted with a quadratic and the intercept extrapolated back
 * to zero. That is what makes this estimator different from batch means rather
 * than a re-spelling of it: it models the correlation instead of trying to
 * batch it away.
 *
 * FALLS BACK TO PLAIN BATCH MEANS, not to an error, whenever the fit cannot be
 * trusted: too few batches, too few regression points, a zero periodogram
 * ordinate (log undefined) or a non-positive intercept. A spectral estimate
 * from a singular fit is not conservative, it is arbitrary. Transcribes
 * `computeSpectralStatisticsInternal`.
 */
inline StatTriple spectral_statistics(const std::vector<double>& obs, int batch_size,
                                      double low_freq_frac) {
    const std::size_t n = obs.size();
    const int m = (batch_size > 0) ? static_cast<int>(n / static_cast<std::size_t>(batch_size)) : 0;
    if (m < 4) return bm_statistics(obs, batch_size);
 
    std::vector<double> means(static_cast<std::size_t>(m), 0.0);
    for (int i = 0; i < m; ++i) {
        double sum = 0.0;
        for (int j = 0; j < batch_size; ++j)
            sum += obs[static_cast<std::size_t>(i * batch_size + j)];
        means[static_cast<std::size_t>(i)] = sum / batch_size;
    }
    double grand = 0.0;
    for (int i = 0; i < m; ++i) grand += means[static_cast<std::size_t>(i)];
    grand /= m;
 
    const int nf = m / 2;
    if (nf < 2) return bm_statistics(obs, batch_size);
    std::vector<double> per(static_cast<std::size_t>(nf), 0.0);
    for (int k = 1; k <= nf; ++k) {
        double re = 0.0, im = 0.0;
        const double f = 2.0 * 3.14159265358979323846 * k / m;
        for (int t = 0; t < m; ++t) {
            const double c = means[static_cast<std::size_t>(t)] - grand;
            re += c * std::cos(f * t);
            im += c * std::sin(f * t);
        }
        per[static_cast<std::size_t>(k - 1)] = (re * re + im * im) / m;
    }
 
    int npts = static_cast<int>(nf * low_freq_frac);
    if (npts < 3) npts = 3;
    if (npts > nf) return bm_statistics(obs, batch_size);
    for (int k = 0; k < npts; ++k)
        if (!(per[static_cast<std::size_t>(k)] > 0.0)) return bm_statistics(obs, batch_size);
 
    // Ordinary least squares for log I(f) = b0 + b1 f + b2 f^2, solved through
    // the 3x3 normal equations by Gaussian elimination: the design is tiny and
    // pulling in a general solver here would only add a failure mode.
    double A[3][4];
    for (int r = 0; r < 3; ++r)
        for (int c = 0; c < 4; ++c) A[r][c] = 0.0;
    for (int k = 0; k < npts; ++k) {
        const double f = 2.0 * 3.14159265358979323846 * (k + 1) / m;
        const double x[3] = {1.0, f, f * f};
        const double y = std::log(per[static_cast<std::size_t>(k)]);
        for (int r = 0; r < 3; ++r) {
            for (int c = 0; c < 3; ++c) A[r][c] += x[r] * x[c];
            A[r][3] += x[r] * y;
        }
    }
    for (int i = 0; i < 3; ++i) {
        int piv = i;
        for (int r = i + 1; r < 3; ++r)
            if (std::fabs(A[r][i]) > std::fabs(A[piv][i])) piv = r;
        if (!(std::fabs(A[piv][i]) > 1e-300)) return bm_statistics(obs, batch_size);
        if (piv != i)
            for (int c = 0; c < 4; ++c) std::swap(A[i][c], A[piv][c]);
        for (int r = 0; r < 3; ++r) {
            if (r == i) continue;
            const double fct = A[r][i] / A[i][i];
            for (int c = i; c < 4; ++c) A[r][c] -= fct * A[i][c];
        }
    }
    const double b0 = A[0][3] / A[0][0];
    const double s0 = std::exp(b0);
    if (!(s0 > 0.0) || !std::isfinite(s0)) return bm_statistics(obs, batch_size);
 
    StatTriple st;
    st.mean = grand;
    st.stderr_ = std::sqrt(2.0 * 3.14159265358979323846 * s0 / m);
    st.df = npts - 3;
    if (st.df < 1) st.df = 1;
    st.ok = true;
    return st;
}
 
/** Dispatch on `cimethod`. */
inline StatTriple ci_statistics(const std::vector<double>& obs, int batch_size,
                                const LdesOptions& o) {
    if (o.cimethod == "bm") return bm_statistics(obs, batch_size);
    if (o.cimethod == "spectral")
        return spectral_statistics(obs, batch_size, o.spectral_low_freq_frac);
    return obm_statistics(obs, batch_size, o.obmoverlap);
}
 
 
/**
 * Convergence-based stopping: batch the four metrics and stop once EVERY
 * active (station, class) pair has reached the requested relative precision.
 *
 * THE CONJUNCTION IS THE POINT. A run stops when the WORST metric at the WORST
 * pair is precise enough, not when the aggregate is: stopping on an average
 * relative precision leaves the lightly loaded stations, whose estimates
 * converge slowest, arbitrarily wrong while the report looks converged.
 *
 * A pair whose mean is effectively zero counts as converged: a relative
 * precision against zero is not a bound, it is a division. Transcribes
 * `checkConvergence` and `isMetricConverged`.
 */
class Convergence {
public:
    void init(std::size_t M, std::size_t K, const LdesOptions& o, std::uint64_t max_events) {
        enabled_ = o.cnvgon;
        if (!enabled_) return;
        tol_ = o.cnvgtol;
        min_batches_ = o.cnvgbatch;
        interval_ = (o.cnvgchk > 0) ? static_cast<std::uint64_t>(o.cnvgchk)
                                    : std::max<std::uint64_t>(1, max_events / 50);
        level_ = (o.confint > 0.0) ? o.confint : 0.95;
        q_.assign(M, std::vector<std::vector<double>>(K));
        u_.assign(M, std::vector<std::vector<double>>(K));
        r_.assign(M, std::vector<std::vector<double>>(K));
        t_.assign(M, std::vector<std::vector<double>>(K));
        start_qt_.assign(M, std::vector<double>(K, 0.0));
        start_bt_.assign(M, std::vector<double>(K, 0.0));
        start_cmp_.assign(M, std::vector<double>(K, 0.0));
        start_rsum_.assign(M, std::vector<double>(K, 0.0));
        start_rcnt_.assign(M, std::vector<double>(K, 0.0));
        nstations_ = M;
        nclasses_ = K;
    }
 
    bool enabled() const { return enabled_; }
    std::uint64_t interval() const { return interval_; }
 
    /** Close the current batch at `now` and record one batch mean per metric. */
    void finalize_batch(const Accum& acc, const std::vector<std::size_t>& nservers, double now) {
        if (!enabled_) return;
        const double dur = now - batch_start_;
        if (!(dur > 0.0)) {
            reset_batch(acc, now);
            return;
        }
        for (std::size_t i = 0; i < nstations_; ++i)
            for (std::size_t k = 0; k < nclasses_; ++k) {
                const double ql = (acc.tot_qlen[i][k] - start_qt_[i][k]) / dur;
                q_[i][k].push_back(ql);
                const double c = acc.util_peak[i] > 0.0
                                     ? acc.util_peak[i]
                                     : static_cast<double>(nservers[i] > 0 ? nservers[i] : 1);
                u_[i][k].push_back((acc.tot_busy[i][k] - start_bt_[i][k]) / (dur * c));
                t_[i][k].push_back((acc.completed[i][k] - start_cmp_[i][k]) / dur);
                const double dn = acc.resp_cnt[i][k] - start_rcnt_[i][k];
                r_[i][k].push_back(dn > 0.0 ? (acc.resp_sum[i][k] - start_rsum_[i][k]) / dn : 0.0);
            }
        reset_batch(acc, now);
    }
 
    /** True once every active pair has reached the tolerance on all four metrics. */
    bool converged(const std::vector<std::vector<bool>>& off) const {
        if (!enabled_) return false;
        if (nstations_ == 0 || nclasses_ == 0) return false;
        if (static_cast<int>(q_[0][0].size()) < min_batches_) return false;
        for (std::size_t i = 0; i < nstations_; ++i)
            for (std::size_t k = 0; k < nclasses_; ++k) {
                if (off[i][k]) continue;
                if (!metric_ok(q_[i][k])) return false;
                if (!metric_ok(u_[i][k])) return false;
                if (!metric_ok(t_[i][k])) return false;
                bool any_positive = false;
                for (double v : r_[i][k])
                    if (v > 0.0) {
                        any_positive = true;
                        break;
                    }
                if (any_positive && !metric_ok(r_[i][k])) return false;
            }
        return true;
    }
 
    int batches() const { return (nstations_ && nclasses_) ? static_cast<int>(q_[0][0].size()) : 0; }
 
private:
    void reset_batch(const Accum& acc, double now) {
        batch_start_ = now;
        for (std::size_t i = 0; i < nstations_; ++i)
            for (std::size_t k = 0; k < nclasses_; ++k) {
                start_qt_[i][k] = acc.tot_qlen[i][k];
                start_bt_[i][k] = acc.tot_busy[i][k];
                start_cmp_[i][k] = acc.completed[i][k];
                start_rsum_[i][k] = acc.resp_sum[i][k];
                start_rcnt_[i][k] = acc.resp_cnt[i][k];
            }
    }
 
    bool metric_ok(const std::vector<double>& b) const {
        const std::size_t n = b.size();
        if (n < 2) return false;
        double sum = 0.0;
        for (double x : b) sum += x;
        const double mean = sum / n;
        // A relative precision against a zero mean is a division, not a bound.
        if (std::fabs(mean) < 1e-12) return true;
        double var = 0.0;
        for (double x : b) var += (x - mean) * (x - mean);
        var /= (n - 1);
        const double se = std::sqrt(var / n);
        const double half = t_critical(level_, static_cast<int>(n) - 1) * se;
        return half / std::fabs(mean) <= tol_;
    }
 
    bool enabled_ = false;
    double tol_ = 0.05, level_ = 0.95, batch_start_ = 0.0;
    int min_batches_ = 20;
    std::uint64_t interval_ = 1;
    std::size_t nstations_ = 0, nclasses_ = 0;
    std::vector<std::vector<std::vector<double>>> q_, u_, r_, t_;
    std::vector<std::vector<double>> start_qt_, start_bt_, start_cmp_, start_rsum_, start_rcnt_;
};
 
/**
 * Fill the half-width matrices of `res` from the post-warmup observation
 * series. Transcribes `computeOBMConfidenceIntervals`.
 *
 * The three series are not treated alike, and that is the reference's design:
 *
 *  - QUEUE LENGTH has one observation per interval already, so it batches
 *    directly;
 *  - THROUGHPUT is stored CUMULATIVE, so it is differenced into per-interval
 *    rates first -- batching the cumulative counts would estimate the variance
 *    of a random walk;
 *  - UTILIZATION is not measured separately at all. Its interval is the
 *    throughput's, carried through the utilization law U = T*S/c (`u_law`, the
 *    S/c of each row, with c the peak capacity and 1 at a delay), as MATLAB's
 *    law and the JAR's TNCI/(mu*c) read it, not through the measured UN/TN;
 *  - RESPONSE TIME is not batched: as the reference's `responseTimeTally`, it is
 *    the iid interval t(n-1) * sd / sqrt(n) over the tallied response times,
 *    reported only above 30 of them, over the same tally `RN` is the mean of.
 *
 * A series shorter than `ciminobs` yields NO interval rather than a wide one:
 * a half-width computed from a handful of batches is not conservative, it is
 * arbitrary.
 */
inline void batch_means_ci(const Observations& obs, const Truncation& tr, const Accum& acc,
                           const LdesOptions& o, LdesResult& res, const Matrix<double>& u_law,
                           const std::vector<std::vector<bool>>* lost_unserved = nullptr) {
    if (o.cimethod == "none" || !(o.confint > 0.0) || obs.time.empty()) return;
    const std::size_t M = obs.nstations, K = obs.nclasses;
    res.QNCI = Matrix<double>(M, K, 0.0);
    res.UNCI = Matrix<double>(M, K, 0.0);
    res.RNCI = Matrix<double>(M, K, 0.0);
    res.TNCI = Matrix<double>(M, K, 0.0);
    const std::size_t t0 = tr.applied ? tr.index : 0;
 
    for (std::size_t i = 0; i < M; ++i) {
        for (std::size_t r = 0; r < K; ++r) {
            if (t0 >= obs.qlen[i][r].size()) continue;
            const std::vector<double> q(obs.qlen[i][r].begin() + static_cast<std::ptrdiff_t>(t0),
                                        obs.qlen[i][r].end());
            if (q.size() >= static_cast<std::size_t>(o.ciminobs)) {
                const int b = std::max(o.ciminbatch, static_cast<int>(std::sqrt(
                                                         static_cast<double>(q.size()))));
                const StatTriple s = ci_statistics(q, b, o);
                if (s.ok) res.QNCI(i, r) = t_critical(o.confint, s.df) * s.stderr_;
            }
 
            const bool little = lost_unserved != nullptr && (*lost_unserved)[i][r];
            const double n = acc.resp_cnt[i][r];
            if (!little && n > 30.0) {
                const double var = (acc.resp_sq[i][r] - acc.resp_sum[i][r] * acc.resp_sum[i][r] / n) / (n - 1.0);
                if (var > 0.0)
                    res.RNCI(i, r) = t_critical(o.confint, static_cast<int>(n) - 1) * std::sqrt(var / n);
            }
 
            std::vector<double> rates;
            for (std::size_t k = t0 + 1; k < obs.cmp[i][r].size(); ++k) {
                const double dt = obs.time[k] - obs.time[k - 1];
                if (dt > 0.0) rates.push_back((obs.cmp[i][r][k] - obs.cmp[i][r][k - 1]) / dt);
            }
            if (rates.size() >= static_cast<std::size_t>(o.ciminobs)) {
                const int b = std::max(o.ciminbatch, static_cast<int>(std::sqrt(
                                                         static_cast<double>(rates.size()))));
                const StatTriple s = ci_statistics(rates, b, o);
                if (s.ok) {
                    const double half = t_critical(o.confint, s.df) * s.stderr_;
                    res.TNCI(i, r) = half;
                    if (u_law.rows() > i && u_law(i, r) > 0.0) res.UNCI(i, r) = half * u_law(i, r);
                }
            }
            // RN = QN/TN where jobs left unserved, so its half-width is the delta-method one from QNCI and TNCI, as
            // the JAR propagates it, not the served-job tally's.
            if (little && res.QN(i, r) > 0.0 && res.TN(i, r) > 0.0) {
                const double rq = res.QNCI(i, r) / res.QN(i, r), rt = res.TNCI(i, r) / res.TN(i, r);
                res.RNCI(i, r) = (res.QN(i, r) / res.TN(i, r)) * std::sqrt(rq * rq + rt * rt);
            }
        }
    }
 
    // HOW LONG THE RUN SHOULD HAVE BEEN, when the caller asked for it. The
    // half-widths just computed pin the ASYMPTOTIC variance of each estimator,
    // which is the quantity a run length is planned from -- not the stationary
    // variance, which on M/M/1 differs from it by a factor blowing up like
    // (1-rho)^-2.
    if (o.run_length_plan_precision > 0.0) {
        // The ACTUAL number of events simulated where the result carries it,
        // not the budget: LDES stops early on convergence, and planning from a
        // budget it never spent would overstate N and so overstate sigma^2.
        const double used = res.total_simulated_events > 0
                                ? static_cast<double>(res.total_simulated_events)
                                : static_cast<double>(o.events > 0 ? o.events : o.samples);
        if (used > 0.0) {
            res.run_length_plan = sim::sim_runlength_plan<double>(
                res.QN, res.QNCI, used, o.run_length_plan_precision, o.confint);
            res.has_run_length_plan = true;
        }
    }
}
 
}  // namespace engine
}  // namespace ldes
}  // namespace line
 
#endif  // LINE_SOLVERS_LDES_LDES_STATS_H