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^{-v} prod_r ((Z_r + L_r) v)^{N_r} dv,
 *
 * whose log-integrand at the grid point v is
 *
 *   f(v) = -v + sum(N) log v + sum_r N_r log(Z_r + L_r).
 *
 * The reference lays a grid of ceil(samples/2) uniform points on [0,1) plus
 * ceil(samples/2) points log-spaced over [1, 1e5], rescaling the demands so
 * that every coefficient is at least one, and reports
 *
 *   lG = max_i (du_i + f(v_i)) - sum_r log N_r! + sum(N) log(scaleFactor).
 *
 * DEFECT IN THE REFERENCE, reproduced here rather than repaired. The estimator
 * takes the MAXIMUM of the weighted log-integrand instead of summing the
 * quadrature contributions log sum_i du_i exp(f(v_i)). Two consequences:
 *
 *  - the width of the integrand is discarded, so the estimate is short by the
 *    Laplace factor sqrt(2 pi sum(N)): with the exact value being
 *    log((sum N)!) + sum_r N_r log(Z_r+L_r) - sum_r log N_r!, the returned
 *    value is Stirling's leading term without the sqrt(2 pi n) correction, a
 *    systematic relative shortfall of about 1/sqrt(2 pi sum(N)) in the log and
 *    a factor sqrt(2 pi sum(N)) in G itself;
 *  - the grid v is NOT sorted (uniform points first, then an increasing
 *    log-spaced block), so du = [0, diff(v)] alternates in sign and is not a
 *    quadrature weight at all; it merely perturbs the argmax.
 *
 * Because of this the routine does not converge to the exact constant as the
 * sample count grows, and its test in this tree asserts the reference's own
 * behaviour and the size of the shortfall, not agreement with pfqn_ca.
 *
 * 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 <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("pfqn_mmsample2: non-positive scale factor, demands must be positive");
 
    // sum_r N_r log((Z_r + L_r)/scale), the v-independent part of f
    double lcoef = 0.0;
    for (std::size_t r = 0; r < R; ++r) {
        if (N[r] == 0) continue;
        const T coef = (Z[r] + L(0, r)) / scale;
        lcoef += N[r] * num_traits<T>::log_as_double(coef);
    }
    long Ntot = 0;
    for (int n : N) Ntot += n;
 
    // ---- the grid: uniform block, then the log-spaced block ----------------
    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));
    }
 
    double best = -std::numeric_limits<double>::infinity();
    for (std::size_t i = 0; i < v.size(); ++i) {
        // du = [0, diff(v)]: the reference's (sign-indefinite) weight term.
        const double du = i == 0 ? 0.0 : v[i] - v[i - 1];
        if (v[i] <= 0.0) continue;  // log(0) would be -inf and never wins the max
        const double f = du - v[i] + static_cast<double>(Ntot) * std::log(v[i]) + lcoef;
        if (f > best) best = f;
    }
 
    double lG = best;
    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