pfqn_mmsample2.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_MMSAMPLE2_H
#define LINE_API_PFQN_MMSAMPLE2_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Sampled McKenna-Mitra integral form of the normalizing constant of a
 * repairman (single-queue plus delay) model.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_mmsample2.m, cross-checked
 * against jar/src/main/java/jline/api/pfqn/nc/Pfqn_mmsample2.java.
 *
 * For a single queue the McKenna-Mitra form collapses to the one-dimensional
 * integral
 *
 *   G(N) = 1/prod_r N_r! int_0^inf e^{-u} prod_r (Z_r + L_r u)^{N_r} du,
 *
 * whose log-integrand at the grid point u is
 *
 *   f(u) = -u + sum_r N_r log(Z_r + L_r u).
 *
 * The reference lays a grid of ceil(samples/2) uniform points on [0,1) plus
 * ceil(samples/2) points log-spaced over [1, 1e5], SORTS it, rescales the
 * demands so that every coefficient is at least one, takes the panel widths
 * du = [v(1), diff(v)] (the first panel covering [0, v(1)]) and reports the
 * quadrature sum in the log domain,
 *
 *   lG = log sum_i du_i exp(f(v_i)) - sum_r log N_r! + sum(N) log(scaleFactor).
 *
 * (An earlier reference took the maximum of the weighted log-integrand over an
 * unsorted grid with the integrand (Z_r + L_r) u; MATLAB has since been
 * corrected and this port follows the corrected estimator.)
 *
 * Only the FIRST row of L is read, matching the reference: this is a
 * single-queue model and additional rows are ignored.
 *
 * Arithmetic: INEXACT BY CONSTRUCTION. The grid is random, the integrand is
 * evaluated in the log domain, and the log-spaced block is generated by
 * exponentiation.
 *
 * RNG contract: see pfqn_mc_common.h. Comparable to MATLAB only in
 * distribution, never stream for stream; reproducible within this port only
 * when the generator is passed in the same state.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/api/pfqn/pfqn_ca.h"
#include "line/api/pfqn/pfqn_mc_common.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/**
 * @brief Sampled McKenna-Mitra integral form of the normalizing constant of a
 *        repairman (single-queue plus delay) model.
 *
 * @param L       (M x R) demands; only row 0 is read, as in the reference
 * @param N       (R) population per class
 * @param Z       (R) think times
 * @param samples grid size; half uniform on [0,1), half log-spaced on [1,1e5]
 * @param rng     explicit generator, advanced by the call
 */
template <class T>
NcResult<T> pfqn_mmsample2(const Matrix<T>& L, const std::vector<int>& N, const std::vector<T>& Z,
                           std::size_t samples, McRng& rng) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_mmsample2 requires transcendental arithmetic: it samples a random "
                  "quadrature grid and evaluates a log-domain integrand on it");
 
    const std::size_t R = N.size();
    if (L.empty()) throw InputError("pfqn_mmsample2: empty demand matrix");
    if (L.cols() != R) throw InputError("pfqn_mmsample2: L and N disagree on the class count");
    if (Z.size() != R) throw InputError("pfqn_mmsample2: Z has the wrong length");
    for (int n : N)
        if (n < 0) throw InputError("pfqn_mmsample2: negative population");
    if (samples == 0) throw InputError("pfqn_mmsample2: at least one grid point is required");
 
    // ---- rescale so that every coefficient is at least one -----------------
    T mn = L(0, 0);
    for (std::size_t i = 0; i < L.rows(); ++i)
        for (std::size_t r = 0; r < R; ++r)
            if (L(i, r) < mn) mn = L(i, r);
    for (std::size_t r = 0; r < R; ++r)
        if (Z[r] < mn) mn = Z[r];
    const T scale = num_traits<T>::from_double(1e-7) + mn;
    if (!(scale > num_traits<T>::from_int(0)))
        throw NumericError(ErrorKind{"Degenerate"}, "pfqn_mmsample2: non-positive scale factor, demands must be positive");
 
    // scaled coefficients Z_r/scale and L_r/scale (first row of L)
    std::vector<double> Zs(R), Ls(R);
    for (std::size_t r = 0; r < R; ++r) {
        Zs[r] = num_traits<T>::to_double(T(Z[r] / scale));
        Ls[r] = num_traits<T>::to_double(T(L(0, r) / scale));
    }
    long Ntot = 0;
    for (int n : N) Ntot += n;
 
    // ---- the grid: uniform block plus the log-spaced block, SORTED ---------
    const std::size_t nu = static_cast<std::size_t>(
        std::ceil(0.5 * static_cast<double>(samples)));
    const std::size_t nl = static_cast<std::size_t>(
        std::ceil(0.5 * static_cast<double>(samples)));
    std::vector<double> v;
    v.reserve(nu + nl);
    for (std::size_t i = 0; i < nu; ++i) v.push_back(mc_uniform01(rng));
    for (std::size_t i = 0; i < nl; ++i) {
        const double t = nl == 1 ? 0.0 : 5.0 * static_cast<double>(i) / static_cast<double>(nl - 1);
        v.push_back(std::pow(10.0, t));
    }
    std::sort(v.begin(), v.end());
 
    // lterms_i = log(du_i) - v_i + sum_r N_r log(Z_r + L_r v_i), du = [v(1), diff(v)]
    std::vector<double> lterms;
    lterms.reserve(v.size());
    double lmax = -std::numeric_limits<double>::infinity();
    for (std::size_t i = 0; i < v.size(); ++i) {
        const double du = i == 0 ? v[0] : v[i] - v[i - 1];
        double f = std::log(du) - v[i];
        for (std::size_t r = 0; r < R; ++r) {
            if (N[r] == 0) continue;
            f += static_cast<double>(N[r]) * std::log(Zs[r] + Ls[r] * v[i]);
        }
        lterms.push_back(f);
        if (f > lmax) lmax = f;
    }
    double ssum = 0.0;
    for (double f : lterms)
        if (std::isfinite(f)) ssum += std::exp(f - lmax);
 
    double lG = lmax + std::log(ssum);
    for (std::size_t r = 0; r < R; ++r) lG -= mc_log_factorial<T>(N[r]);
    lG += static_cast<double>(Ntot) * num_traits<T>::log_as_double(scale);
    return {mc_exp<T>(lG), lG};
}
 
/** Reference call shape with an explicit grid size. */
template <class T>
NcResult<T> pfqn_mmsample2(const Matrix<T>& L, const std::vector<int>& N, const std::vector<T>& Z,
                           McRng& rng) {
    return pfqn_mmsample2(L, N, Z, static_cast<std::size_t>(100000), rng);
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_MMSAMPLE2_H