pfqn_ble.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_PFQN_BLE_H
#define LINE_API_PFQN_PFQN_BLE_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Logistic expansion with the eps->0 bias correction (BLE).
 *
 * Templated port of matlab/src/api/pfqn/pfqn_ble.m. Cas17 Theorem 4.1 holds for
 * eps >= eps_N > 0; the K(1+eps*N) self-looping populations are what make the
 * integrand concentrate. Evaluated at eps->0, as pfqn_le does, the curvature at
 * the saddle tends to 1 rather than growing with N, so Laplace's method has no
 * asymptotic regime there and carries an O(1) relative bias per Laplaced
 * direction. With kappa = r(1) = 1 - log(2 pi)/2 and r(a) the Stirling remainder
 *   r(a) = gammaln(a) - ((a - 1/2) log a - a + log(2 pi)/2),
 * the two branch constants are (as in MATLAB pfqn_ble.m):
 *   Z = 0 : M*kappa - r(M), the M-1 simplex directions plus the Jacobian, the
 *           radial integral being exact as Gamma(N+M);
 *   Z > 0 : M*kappa, the radial direction is Laplaced in t = log(v) and
 *           contributes exactly +r(M), cancelling the -r(M).
 * The published expansion is NOT in error; see _kb/03-api-layer.md.
 *
 * ARITHMETIC. Inherited from pfqn_le: the correction is a logarithm, so the
 * routine is gated on num_traits<T>::has_transcendental for the same reason.
 *
 * DEGENERATE BRANCH. When there is no queueing station to expand around, pfqn_le
 * returns the exact delay term and there is no Laplace step to correct, so
 * pfqn_ble returns it unchanged; this is the MATLAB branching.
 */
 
#include <cmath>
#include <cstddef>
#include <vector>
 
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/api/pfqn/pfqn_le.h"
#include "line/lang/lang_types.h"
#include "line/num/number.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** Return value of pfqn_ble, mirroring [Gn, lGn]. */
template <class T>
using BleResult = LeResult<T>;
 
/**
 * Logistic expansion estimate of the normalizing constant, bias-corrected.
 *
 * @param L (M x R) demands, @param N (R) population, @param Z (R) think times
 *          (pass an empty vector or all zeros for the Z = 0 branch)
 */
template <class T>
BleResult<T> pfqn_ble(const Matrix<T>& L, const std::vector<T>& N, const std::vector<T>& Z) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_ble requires transcendental arithmetic (Laplace approximation of an integral)");
    using std::exp;
    using std::log;
    const std::size_t M = L.rows(), R = L.cols();
    const T zero = num_traits<T>::from_int(0);
 
    BleResult<T> res = pfqn_le(L, N, Z);
 
    T Ntot = zero, Lsum = zero, Zsum = zero;
    for (const T& x : N) Ntot += x;
    for (const T& x : Z) Zsum += x;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < R; ++r) Lsum += L(i, r);
    if (M == 0 || N.empty() || Ntot == zero || num_traits<T>::to_double(Lsum) < 1e-4) {
        return res;
    }
 
    // Same predicate pfqn_le branches on, so the constant always matches the branch.
    const bool no_delay =
        Z.empty() || num_traits<T>::to_double(Zsum) < lang::GlobalConstants::Zero;
    const T one = num_traits<T>::from_int(1);
    const T twopi = num_traits<T>::from_double(6.283185307179586476925286766559);
    const T kappa = T(one - log(twopi) / num_traits<T>::from_int(2));
    const T Mt = num_traits<T>::from_int(static_cast<long>(M));
    // MATLAB: Z=0 adds M*kappa - stirlingrem(M); Z>0 adds M*kappa.
    res.lG += Mt * kappa;
    if (no_delay) {
        // stirlingrem(a) = gammaln(a) - ((a-0.5)*log(a) - a + 0.5*log(2*pi))
        const T half = num_traits<T>::from_rational(1, 2);
        const T rem = T(detail::num_lgamma<T>(Mt) - T(T(Mt - half) * log(Mt) - Mt + half * log(twopi)));
        res.lG -= rem;
    }
    res.G = exp(res.lG);
    return res;
}
 
template <class T>
BleResult<T> pfqn_ble(const Matrix<T>& L, const std::vector<T>& N) {
    return pfqn_ble(L, N, std::vector<T>());
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_BLE_H