pfqn_divdiff_ld.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_DIVDIFF_LD_H
#define LINE_API_PFQN_PFQN_DIVDIFF_LD_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Exact load-dependent normalizing constant by divided differences of a
 * single-class constant.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_divdiff_ld.m. Evaluates the
 * multiclass load-dependent normalizing constant as the N-th order finite
 * difference, in the class populations, of a SINGLE-CLASS load-dependent
 * constant:
 *
 *   G(N) = sum_{0<=n<=N} (-1)^(|N|-|n|) / (N_1! ... N_R!)
 *                        * prod_r nchoosek(N_r,n_r) * Gld_n(|N|)
 *
 * where |N| = sum_r N_r, |n| = sum_r n_r and Gld_n(|N|) is the normalizing
 * constant of the single-class model over the same M stations, with the same
 * load-dependent rates alpha_i(.) = mu(i,.), total population |N| and
 * aggregated demands rho_i(n) = sum_r n_r L(i,r). Since
 * nchoosek(N_r,n_r)/N_r! = 1/(n_r! (N_r-n_r)!), the coefficient is evaluated
 * in the equivalent factorial form, which needs no binomial.
 *
 * The identity holds for arbitrary rate functions alpha_i(.), hence it covers
 * load-independent, infinite-server, multiserver and limited load-dependent
 * stations alike. Cost is prod_r (N_r+1) evaluations of pfqn_gldsingle.
 *
 * NUMERICS. The sum alternates in sign and its terms are much larger than the
 * result, so it is evaluated as a signed log-sum-exp, which removes the range
 * problem but not the cancellation. The reference raises warnings past 15
 * decimal digits lost and when the sum comes out negative; warnings become
 * outputs on this side: lossDigits is reported and a negative sum returns NaN
 * exactly as MATLAB does, so the routine is gated on
 * num_traits<T>::has_transcendental like pfqn_explicit_ld.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/api/pfqn/pfqn_asympt_common.h"
#include "line/api/pfqn/pfqn_explicit.h"
#include "line/api/pfqn/pfqn_gldsingle.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
#include "line/util/population.h"
 
namespace line {
namespace pfqn {
 
/** Return value of pfqn_divdiff_ld, mirroring [lG, G, lossDigits]. */
template <class T>
struct DivdiffLdResult {
    T lG;                     ///< logarithm of the normalizing constant
    T G;                      ///< the normalizing constant
    double lossDigits = 0.0;  ///< decimal digits lost to cancellation in the alternating sum
};
 
/**
 * @brief Exact load-dependent normalizing constant by divided differences of a
 *        single-class constant.
 *
 * @param L  (M x R) service demands
 * @param N  (R) population vector
 * @param Z  (D x R) think times, folded in as an infinite server; empty for none
 * @param mu (M x >= sum(N)) load-dependent rates, alpha_i(j) = mu(i,j); empty
 *           for load-independent stations (all ones)
 */
template <class T>
DivdiffLdResult<T> pfqn_divdiff_ld(const Matrix<T>& L, const std::vector<int>& N,
                                   const Matrix<T>& Z, const Matrix<T>& mu) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_divdiff_ld requires transcendental arithmetic: the alternating sum is "
                  "carried as a signed log-sum-exp. Use pfqn_gld for the same constant in exact "
                  "arithmetic");
    using std::exp;
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
    const double dinf = std::numeric_limits<double>::infinity();
    const double dnan = std::numeric_limits<double>::quiet_NaN();
    const std::size_t R = N.size();
 
    DivdiffLdResult<T> res;
    res.lossDigits = 0.0;
    long Nt = 0;
    for (std::size_t r = 0; r < R; ++r) Nt += N[r];
    if (Nt < 0) {
        res.lG = num_traits<T>::from_double(-dinf);
        res.G = zero;
        return res;
    }
    if (Nt == 0) {
        res.lG = zero;
        res.G = one;
        return res;
    }
    const std::size_t M0 = L.empty() ? 0 : L.rows();
    const std::size_t Ntu = static_cast<std::size_t>(Nt);
 
    // mu defaults to ones(M, Nt); otherwise checked and trimmed to Nt columns so
    // that every rate used downstream is positive.
    Matrix<T> mu2(M0, Ntu, one);
    if (!mu.empty()) {
        if (mu.rows() != M0)
            throw InputError(
                "pfqn_divdiff_ld: the load-dependent rate matrix must have one row per station "
                "of L.");
        if (mu.cols() < Ntu)
            throw InputError(
                "pfqn_divdiff_ld: the load-dependent rate matrix must have at least sum(N) "
                "columns.");
        for (std::size_t i = 0; i < M0; ++i)
            for (std::size_t j = 0; j < Ntu; ++j) mu2(i, j) = mu(i, j);
    }
 
    // Z folded in as an infinite-server station with rates 1..Nt
    T Zsum = zero;
    std::vector<T> Zrow(R, zero);
    if (!Z.empty()) {
        for (std::size_t d = 0; d < Z.rows(); ++d)
            for (std::size_t r = 0; r < Z.cols() && r < R; ++r) {
                Zrow[r] += Z(d, r);
                Zsum += Z(d, r);
            }
    }
    const bool withZ = !Z.empty() && Zsum > zero;
    const std::size_t M = M0 + (withZ ? 1 : 0);
    Matrix<T> L2(M, R, zero);
    for (std::size_t i = 0; i < M0; ++i)
        for (std::size_t r = 0; r < R; ++r) L2(i, r) = L(i, r);
    Matrix<T> mu3(M, Ntu, one);
    for (std::size_t i = 0; i < M0; ++i)
        for (std::size_t j = 0; j < Ntu; ++j) mu3(i, j) = mu2(i, j);
    if (withZ) {
        for (std::size_t r = 0; r < R; ++r) L2(M0, r) = Zrow[r];
        for (std::size_t j = 0; j < Ntu; ++j)
            mu3(M0, j) = num_traits<T>::from_int(static_cast<long>(j) + 1);
    }
 
    if (M == 0) {
        res.lG = num_traits<T>::from_double(-dinf);
        res.G = zero;
        return res;
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < R; ++r)
            if (L2(i, r) < zero)
                throw InputError("pfqn_divdiff_ld: the demand matrix must be nonnegative.");
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < Ntu; ++j)
            if (!(mu3(i, j) > zero))
                throw InputError(
                    "pfqn_divdiff_ld: the load-dependent rates must be strictly positive.");
 
    if (R == 1) {
        // in a single class the difference operator is the identity, so the
        // recursion is used directly and no cancellation is incurred
        const NcResult<T> g = pfqn_gldsingle(L2, static_cast<int>(Nt), mu3);
        res.lG = num_traits<T>::from_double(g.lG);
        res.G = exp(res.lG);
        return res;
    }
 
    // ---- alternating sum over the sub-populations 0 <= n <= N ----
    std::vector<T> lterm;
    std::vector<double> sterm;
    std::vector<int> n(R, 0);
    const T minf = num_traits<T>::from_double(-dinf);
    while (true) {
        long sn = 0;
        for (std::size_t r = 0; r < R; ++r) sn += n[r];
        if (sn > 0) {
            // n = 0 leaves every induced demand at zero, so its single-class
            // constant vanishes at sum(N)>0 and the term is left out of the sum
            Matrix<T> Ln(M, 1, zero);
            for (std::size_t i = 0; i < M; ++i)
                for (std::size_t r = 0; r < R; ++r)
                    Ln(i, 0) += L2(i, r) * num_traits<T>::from_int(n[r]);
            const NcResult<T> g = pfqn_gldsingle(Ln, static_cast<int>(Nt), mu3);
            T lt = num_traits<T>::from_double(g.lG);
            for (std::size_t r = 0; r < R; ++r) {
                lt -= detail::num_factln<T>(num_traits<T>::from_int(n[r]));
                lt -= detail::num_factln<T>(num_traits<T>::from_int(N[r] - n[r]));
            }
            lterm.push_back(lt);
            sterm.push_back(((Nt - sn) % 2 == 0) ? 1.0 : -1.0);
        } else {
            lterm.push_back(minf);
            sterm.push_back(0.0);
        }
        if (!next_pop(n, N)) break;
    }
 
    const detail::SignedLse<T> s = detail::explicit_signed_logsumexp(lterm, sterm);
    res.lossDigits = s.lossDigits;
    if (s.sgn == 0) {
        res.lG = minf;
        res.G = zero;
        return res;
    } else if (s.sgn < 0) {
        // The alternating sum returned a negative value: double precision is
        // exhausted by cancellation. MATLAB warns and returns NaN; lossDigits
        // carries the diagnostic here.
        res.lG = num_traits<T>::from_double(dnan);
        res.G = num_traits<T>::from_double(dnan);
        return res;
    }
    res.lG = s.lS;
    res.G = exp(res.lG);
    // lossDigits > 15: the reference warns that the result is unreliable in double
    // precision; the count is returned for the caller to act on.
    return res;
}
 
template <class T>
DivdiffLdResult<T> pfqn_divdiff_ld(const Matrix<T>& L, const std::vector<int>& N,
                                   const Matrix<T>& Z) {
    return pfqn_divdiff_ld(L, N, Z, Matrix<T>());
}
 
template <class T>
DivdiffLdResult<T> pfqn_divdiff_ld(const Matrix<T>& L, const std::vector<int>& N) {
    return pfqn_divdiff_ld(L, N, Matrix<T>(), Matrix<T>());
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_DIVDIFF_LD_H