pfqn_lld.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_LLD_H
#define LINE_API_PFQN_PFQN_LLD_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Normalizing constant of a multiclass LIMITED load-dependent closed model.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_lld.m. Same recursion, same
 * arithmetic and the same result as the reference's pfqn_gld, but with the
 * rate shift saturated at the LIMITED LOAD-DEPENDENCE threshold, which turns
 * the recursion's state space finite and lets it be memoised:
 *
 *   g(m,n,j) = g(m-1,n,0) + sum_r L(m,r)/alpha_m(j+1) * g(m,n-e_r,min(j+1,s_m-1))
 *
 * with j the number of shifts row m has taken. Once j >= s_m-1, where s_m is
 * the population past which alpha_m stays constant, a further shift leaves the
 * row unchanged over the columns the recursion can still read, so the state
 * (m,n,j) repeats and one memo answers what the unmemoised recursion recomputes
 * down an exponential tree. The state space is M * prod_r(N_r+1) * max_k s_k.
 *
 * Every shortcut of the reference's pfqn_gld cascade is kept and evaluated on
 * the same materialised arguments: the single-station closed form, the
 * single-class kernel (pfqn_lldsingle), and the exact fixed-rate constant
 * (pfqn_ca, MATLAB's pfqn_nc 'exact') when no station of the block is load
 * dependent. The third output s holds the detected per-station thresholds.
 *
 * ARITHMETIC. Field operations only: the closed form at m = 1 is formed as a
 * multinomial times integer powers rather than through exp/log, so the routine
 * is EXACT in rational arithmetic, like pfqn_gld and pfqn_lldsingle here. The
 * threshold scan tests exact equality first and then an eps-relative tie on
 * FINITE pairs, as the reference does.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <map>
#include <vector>
 
#include "line/api/pfqn/pfqn_amva_common.h"
#include "line/api/pfqn/pfqn_ca.h"
#include "line/api/pfqn/pfqn_lldsingle.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** Return value of pfqn_lld, mirroring [G, lG, s]. */
template <class T>
struct LldResult {
    T G;                 ///< normalizing constant in the requested arithmetic
    double lG;           ///< log of the constant
    std::vector<int> s;  ///< (M) per-station thresholds, alpha_i(n) = alpha_i(s_i) for n >= s_i
};
 
/**
 * @brief Normalizing constant of a multiclass LIMITED load-dependent closed model.
 *
 * @param L  (M x R) service demands
 * @param N  (R) population vector
 * @param mu (M x >= sum(N)) load-dependent rates, alpha_i(j) = mu(i,j); empty for
 *           all ones (every station a single server)
 */
template <class T>
LldResult<T> pfqn_lld(const Matrix<T>& L, const std::vector<int>& N, const Matrix<T>& mu) {
    const std::size_t M = L.rows();
    const std::size_t R = N.size();
    if (!L.empty() && L.cols() != R)
        throw InputError("pfqn_lld: demand matrix and population vector disagree on the class count");
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
 
    long Ntot = 0;
    for (int v : N) {
        if (v < 0) throw InputError("pfqn_lld: negative population");
        Ntot += v;
    }
    const std::size_t Ntu = static_cast<std::size_t>(Ntot);
 
    // As in pfqn_gld, the mu default has to precede every branch that reads mu.
    Matrix<T> mu2 = mu.empty() ? Matrix<T>(M, Ntu, one) : mu;
    if (!mu2.empty() && mu2.rows() != M) throw InputError("pfqn_lld: mu has the wrong station count");
    const std::size_t ncols = mu2.empty() ? 0 : mu2.cols();
 
    // ---- s_k: the smallest offset past which the rate row is constant ----
    // Equality is tested first so that an infinite rate ties with itself instead
    // of producing Inf-Inf; the tolerance is confined to FINITE pairs.
    std::vector<int> s(M, 1);
    if (ncols >= 1) {
        const double eps = std::numeric_limits<double>::epsilon();
        for (std::size_t m = 0; m < M; ++m) {
            s[m] = static_cast<int>(ncols);
            const T tail = mu2(m, ncols - 1);
            const double tl = num_traits<T>::to_double(tail);
            for (std::size_t n = ncols; n >= 2; --n) {
                const T prev = mu2(m, n - 2);
                const double pv = num_traits<T>::to_double(prev);
                const bool tie = prev == tail ||
                                 (std::isfinite(tl) && std::isfinite(pv) &&
                                  std::fabs(pv - tl) <= eps * std::max(std::fabs(tl), 1.0));
                if (tie)
                    s[m] = static_cast<int>(n - 1);
                else
                    break;
            }
        }
    }
    for (std::size_t m = 0; m < M; ++m)
        if (s[m] < 1) s[m] = 1;
 
    // memo over (m, j, n)
    std::map<std::vector<int>, T> memo;
 
    // materialize(m, n, j): rows 0..m-2 are never shifted, row m-1 is shifted j
    // times, every row truncated by one column per job already placed.
    const auto materialize = [&](std::size_t m, const std::vector<int>& n, int j, Matrix<T>& Lm,
                                 Matrix<T>& mum) {
        long sn = 0;
        for (int v : n) sn += v;
        const long cols = static_cast<long>(ncols) - (Ntot - sn);
        Lm = Matrix<T>(m, R, zero);
        for (std::size_t i = 0; i < m; ++i)
            for (std::size_t r = 0; r < R; ++r) Lm(i, r) = L(i, r);
        const std::size_t cu = cols > 0 ? static_cast<std::size_t>(cols) : 0;
        mum = Matrix<T>(m, cu, zero);
        if (cols > 0 && m > 0) {
            for (std::size_t i = 0; i + 1 < m; ++i)
                for (std::size_t c = 0; c < cu; ++c) mum(i, c) = mu2(i, c);
            for (std::size_t c = 0; c < cu; ++c)
                mum(m - 1, c) = mu2(m - 1, static_cast<std::size_t>(j) + c);
        }
    };
 
    std::function<T(std::size_t, const std::vector<int>&, int)> node;
    std::function<T(const Matrix<T>&, const std::vector<int>&, const Matrix<T>&, std::size_t, int)>
        gld_node;
 
    node = [&](std::size_t m, const std::vector<int>& n, int j) -> T {
        std::vector<int> key;
        key.reserve(2 + R);
        key.push_back(static_cast<int>(m));
        key.push_back(j);
        for (int v : n) key.push_back(v);
        const auto it = memo.find(key);
        if (it != memo.end()) return it->second;
        Matrix<T> Lm, mum;
        materialize(m, n, j, Lm, mum);
        const T g = gld_node(Lm, n, mum, m, j);
        memo[key] = g;
        return g;
    };
 
    gld_node = [&](const Matrix<T>& Lm, const std::vector<int>& n, const Matrix<T>& mum,
                   std::size_t m, int j) -> T {
        long sn = 0;
        for (int v : n) sn += v;
        if (m == 1) {
            for (std::size_t r = 0; r < R; ++r)
                if (n[r] > 0 && Lm(0, r) == zero) return zero;
            // exp(factln(sum n) - sum factln(n) + sum_r n_r log L_r - sum_k log mu(1,k))
            T g = num_multinomial<T>(n);
            for (std::size_t r = 0; r < R; ++r)
                if (n[r] > 0) g *= num_pow_int(Lm(0, r), static_cast<unsigned>(n[r]));
            for (long k = 0; k < sn; ++k) g = T(g / mum(0, static_cast<std::size_t>(k)));
            return g;
        }
        if (R == 1) {
            return pfqn_lldsingle(Lm, static_cast<int>(sn), mum).G;
        }
        if (Lm.empty()) return zero;
        bool isLoadDep = false;
        std::vector<bool> isInfServer(m, false);
        for (std::size_t ist = 0; ist < m; ++ist) {
            bool allOne = true, allRamp = true;
            for (long k = 0; k < sn; ++k) {
                const T& v = mum(ist, static_cast<std::size_t>(k));
                if (v != one) allOne = false;
                if (v != num_traits<T>::from_int(k + 1)) allRamp = false;
            }
            // MATLAB: min(row)==1 & max(row)==1 is [] on an empty row, so the
            // branch falls through to all([]) = true, the infinite-server arm.
            if (sn > 0 && allOne) {
                isInfServer[ist] = false;
            } else if (allRamp) {
                isInfServer[ist] = true;
            } else {
                isLoadDep = true;
            }
        }
        if (!isLoadDep) {
            std::size_t nli = 0;
            for (std::size_t ist = 0; ist < m; ++ist)
                if (!isInfServer[ist]) ++nli;
            Matrix<T> Lli(nli > 0 ? nli : 1, R, zero);
            Matrix<T> Zli(1, R, zero);
            std::size_t q = 0;
            for (std::size_t ist = 0; ist < m; ++ist) {
                if (isInfServer[ist]) {
                    for (std::size_t r = 0; r < R; ++r) Zli(0, r) += Lm(ist, r);
                } else {
                    for (std::size_t r = 0; r < R; ++r) Lli(q, r) = Lm(ist, r);
                    ++q;
                }
            }
            return pfqn_ca(Lli, n, Zli).G;
        }
        if (m == 0) return zero;
        if (sn == 0) return one;
        // The recursion itself, memoised through node(). The shift saturates at
        // s(m)-1, past which the row the child would see is the one it sees now.
        T g = node(m - 1, n, 0);
        for (std::size_t r = 0; r < R; ++r) {
            if (n[r] > 0) {
                const int jn = std::min(j + 1, s[m - 1] - 1);
                g += T(T(Lm(m - 1, r) / mum(m - 1, 0)) * node(m, oner(n, r + 1), jn));
            }
        }
        return g;
    };
 
    LldResult<T> res;
    res.s = s;
    res.G = node(M, N, 0);
    res.lG = num_traits<T>::log_as_double(res.G);
    return res;
}
 
/** Overload with all rates equal to one, i.e. every station a single server. */
template <class T>
LldResult<T> pfqn_lld(const Matrix<T>& L, const std::vector<int>& N) {
    return pfqn_lld(L, N, Matrix<T>());
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_LLD_H