lossn_kr.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_LOSSN_LOSSN_KR_H
#define LINE_API_LOSSN_LOSSN_KR_H
 
/**
 * @file
 * @ingroup api_lossn
 * The multirate single-link loss system (stochastic knapsack): exact blocking
 * by the Kaufman-Roberts recursion, dimensioning against per-class blocking
 * targets, and the critical-loading normal approximation.
 *
 * Templated ports of matlab/src/api/lossn/lossn_kr.m,
 * lossn_kr_dim.m and lossn_kr_asy.m.
 *
 * A single link of C bandwidth units is offered K classes of calls. Class k
 * arrives as a Poisson process with offered load rho_k (arrival rate times
 * mean holding time) and holds b_k units for the whole call; a call that finds
 * fewer than b_k free units is lost. The occupancy distribution
 * q(n) = P(n units busy), n = 0..C, satisfies the Kaufman-Roberts recursion
 *
 *      n q(n) = sum_k b_k rho_k q(n - b_k),   q(n) = 0 for n < 0,
 *
 * and a class-k call is blocked when more than C - b_k units are busy,
 * B_k = sum_{n = C-b_k+1}^{C} q(n).
 *
 * ARITHMETIC. The recursion and the dimensioning search are rational (the
 * 1e250 rescale only compares a double image of the running term), so
 * lossn_kr and lossn_kr_dim instantiate at every
 * arithmetic. The critical-loading approximation evaluates the normal density
 * and distribution function and is gated on num_traits<T>::has_transcendental.
 *
 * Reference: J. S. Kaufman, "Blocking in a Shared Resource Environment", IEEE
 * Trans. Commun. 29(10), 1981. J. W. Roberts, "A Service System with
 * Heterogeneous User Requirements", 1981. K. W. Ross, "Multiservice Loss Models
 * for Broadband Telecommunication Networks", Springer, 1995, Sec. 2.2-2.7.
 * M. I. Reiman, "A critically loaded multirate link with trunk reservation",
 * Queueing Systems 9, 1991.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace lossn {
 
/** Result of lossn_kr. */
template <class T>
struct KrResult {
    std::vector<T> q;     ///< occupancy distribution, q[n] = P(n units busy), n = 0..C
    std::vector<T> B;     ///< blocking probability of class k
    std::vector<T> QLen;  ///< mean number of class-k calls in progress, rho_k (1 - B_k)
};
 
/** Result of lossn_kr_dim. */
template <class T>
struct KrDimResult {
    long C = 0;           ///< smallest capacity with B_k(C) <= Btarget_k for all k
    std::vector<T> B;     ///< blocking probabilities at that capacity
};
 
/** Result of lossn_kr_asy. */
template <class T>
struct KrAsyResult {
    std::vector<T> B;  ///< approximate blocking probability of class k
    T alpha;           ///< normalized spare capacity (C - sum b rho)/sqrt(C)
    T sigma;           ///< normalized occupancy standard deviation
    T delta;           ///< blocking per unit of call size, times sqrt(C)
};
 
namespace detail {
 
/** The shared argument checks of the three functions: sizes and signs of rho and b. */
template <class T>
void kr_check_rho_b(const char* fn, const std::vector<T>& rho, const std::vector<T>& b) {
    if (b.size() != rho.size())
        throw InputError(std::string(fn) + ": b has " + std::to_string(b.size()) +
                         " entries but rho has " + std::to_string(rho.size()) + ".");
    const T zero = num_traits<T>::from_int(0);
    for (const T& r : rho)
        if (r < zero) throw InputError(std::string(fn) + ": rho must be nonnegative.");
}
 
/** b checked to be nonnegative integers within 1e-9, then rounded. */
template <class T>
std::vector<long> kr_integral_b(const char* fn, const std::vector<T>& b) {
    std::vector<long> out(b.size(), 0);
    for (std::size_t k = 0; k < b.size(); ++k) {
        const double x = num_traits<T>::to_double(b[k]);
        if (!(x >= 0.0) || std::fabs(x - std::round(x)) > 1e-9)
            throw InputError(std::string(fn) + ": b must contain nonnegative integers.");
        out[k] = static_cast<long>(std::llround(x));
    }
    return out;
}
 
/** min(1, x), as MATLAB caps the blocking sums. */
template <class T>
T kr_cap1(const T& x) {
    const T one = num_traits<T>::from_int(1);
    return x > one ? one : x;
}
 
/**
 * phi(x)/Phi(x). For x far in the lower tail both factors underflow, so the
 * ratio is taken from the asymptotic Mills-ratio series instead.
 */
template <class T>
T kr_normal_hazard(const T& x) {
    using std::erfc;
    using std::exp;
    using std::sqrt;
    if (x >= num_traits<T>::from_int(-35)) {
        const T two = num_traits<T>::from_int(2);
        const T pi = num_traits<T>::from_double(3.14159265358979323846);
        return T(exp(T(-(x * x) / two)) / sqrt(T(two * pi)) /
                 T(num_traits<T>::from_double(0.5) * erfc(T(-x / sqrt(two)))));
    }
    const T one = num_traits<T>::from_int(1);
    const T x2 = T(x * x), x4 = T(x2 * x2), x6 = T(x4 * x2);
    return T(-x / T(one - one / x2 + num_traits<T>::from_int(3) / x4 -
                     num_traits<T>::from_int(15) / x6));
}
 
}  // namespace detail
 
/**
 * @brief Exact occupancy distribution and per-class blocking of a stochastic knapsack.
 *
 * The result is insensitive to the holding-time distributions beyond their
 * means. A class with b_k = 0 is never blocked; a class with b_k > C always is.
 *
 * @param rho offered load of class k (K, nonnegative)
 * @param b   bandwidth units held by a class-k call (K, nonnegative integers)
 * @param C   link capacity in bandwidth units (nonnegative integer)
 * @return q (C+1), B (K) and QLen = rho (1 - B) (K)
 */
template <class T>
KrResult<T> lossn_kr(const std::vector<T>& rho, const std::vector<T>& b,
                                             const T& C) {
    const char* fn = "lossn_kr";
    detail::kr_check_rho_b(fn, rho, b);
    const std::vector<long> bi = detail::kr_integral_b(fn, b);
    const double cd = num_traits<T>::to_double(C);
    if (!(cd >= 0.0) || std::fabs(cd - std::round(cd)) > 1e-9)
        throw InputError(std::string(fn) + ": C must be a nonnegative integer.");
    const long Ci = static_cast<long>(std::llround(cd));
    const std::size_t K = rho.size();
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
 
    KrResult<T> out;
    std::vector<T>& q = out.q;
    q.assign(static_cast<std::size_t>(Ci + 1), zero);
    q[0] = one;
    for (long n = 1; n <= Ci; ++n) {
        T s = zero;
        for (std::size_t k = 0; k < K; ++k)
            if (bi[k] > 0 && n - bi[k] >= 0)
                s += T(num_traits<T>::from_int(bi[k]) * rho[k] *
                       q[static_cast<std::size_t>(n - bi[k])]);
        q[static_cast<std::size_t>(n)] = T(s / num_traits<T>::from_int(n));
        if (num_traits<T>::to_double(q[static_cast<std::size_t>(n)]) > 1e250) {
            // rescale; only the normalized distribution is used
            const T scale = q[static_cast<std::size_t>(n)];
            for (long m = 0; m <= n; ++m)
                q[static_cast<std::size_t>(m)] = T(q[static_cast<std::size_t>(m)] / scale);
        }
    }
    T tot = zero;
    for (const T& v : q) tot += v;
    for (T& v : q) v = T(v / tot);
 
    out.B.assign(K, zero);
    out.QLen.assign(K, zero);
    for (std::size_t k = 0; k < K; ++k) {
        if (bi[k] > 0) {
            T acc = zero;
            for (long n = std::max<long>(0, Ci - bi[k] + 1); n <= Ci; ++n)
                acc += q[static_cast<std::size_t>(n)];
            out.B[k] = detail::kr_cap1(acc);
        }
        out.QLen[k] = T(rho[k] * T(one - out.B[k]));
    }
    return out;
}
 
/** lossn_kr with an integer capacity. */
template <class T>
KrResult<T> lossn_kr(const std::vector<T>& rho, const std::vector<T>& b,
                                             long C) {
    return lossn_kr(rho, b, num_traits<T>::from_int(C));
}
 
/**
 * @brief Smallest link capacity meeting per-class blocking targets.
 *
 * Returns the smallest integer capacity C such that every class of the
 * stochastic knapsack has blocking B_k(C) no larger than its target. The
 * unnormalized Kaufman-Roberts terms n q(n) = sum_k b_k rho_k q(n - b_k),
 * q(0) = 1, do not depend on C, so the search extends a single recursion one
 * unit at a time and evaluates B_k(C) = sum_{n>C-b_k} q(n) / sum_{n<=C} q(n) at
 * each step, at a total cost of O(C K). The search scans every C upwards rather
 * than bisecting because the blocking of a wide class need not decrease
 * monotonically in C (Ross 1995, Sec. 2.6).
 *
 * @param rho     offered load of class k (K, nonnegative)
 * @param b       bandwidth units held by a class-k call (K, nonnegative integers)
 * @param Btarget blocking target per class (K), each in (0,1]
 */
template <class T>
KrDimResult<T> lossn_kr_dim(const std::vector<T>& rho,
                                                    const std::vector<T>& b,
                                                    const std::vector<T>& Btarget) {
    const char* fn = "lossn_kr_dim";
    detail::kr_check_rho_b(fn, rho, b);
    const std::vector<long> bi = detail::kr_integral_b(fn, b);
    const std::size_t K = rho.size();
    if (Btarget.size() != K)
        throw InputError(std::string(fn) + ": Btarget has " + std::to_string(Btarget.size()) +
                         " entries but rho has " + std::to_string(K) + ".");
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    for (const T& t : Btarget)
        if (!(t > zero) || t > one) throw InputError(std::string(fn) + ": Btarget must lie in (0,1].");
 
    const long Cmax = 100000000L;  // far beyond any knapsack the recursion can be asked to size
    std::vector<T> q(1024, zero);
    q[0] = one;
    T total = one;
    long C = 0;
    KrDimResult<T> out;
    while (true) {
        std::vector<T> B(K, zero);
        bool ok = true;
        for (std::size_t k = 0; k < K; ++k) {
            if (bi[k] > 0) {
                T acc = zero;
                for (long n = std::max<long>(0, C - bi[k] + 1); n <= C; ++n)
                    acc += q[static_cast<std::size_t>(n)];
                B[k] = detail::kr_cap1(T(acc / total));
            }
            if (B[k] > Btarget[k]) ok = false;
        }
        if (ok) {
            out.C = C;
            out.B = B;
            return out;
        }
        if (C >= Cmax)
            throw NumericError(ErrorKind{"Capacity"},
                std::string(fn) + ": No capacity up to " + std::to_string(Cmax) +
                " meets the blocking targets.");
        ++C;
        if (static_cast<std::size_t>(C + 1) > q.size()) q.resize(2 * q.size(), zero);
        T s = zero;
        for (std::size_t k = 0; k < K; ++k)
            if (bi[k] > 0 && C - bi[k] >= 0)
                s += T(num_traits<T>::from_int(bi[k]) * rho[k] *
                       q[static_cast<std::size_t>(C - bi[k])]);
        q[static_cast<std::size_t>(C)] = T(s / num_traits<T>::from_int(C));
        total += q[static_cast<std::size_t>(C)];
        if (num_traits<T>::to_double(q[static_cast<std::size_t>(C)]) > 1e250) {
            // only ratios of the terms are used
            const T scale = q[static_cast<std::size_t>(C)];
            for (long m = 0; m <= C; ++m)
                q[static_cast<std::size_t>(m)] = T(q[static_cast<std::size_t>(m)] / scale);
            total = T(total / scale);
        }
    }
}
 
/** lossn_kr_dim with one blocking target shared by every class. */
template <class T>
KrDimResult<T> lossn_kr_dim(const std::vector<T>& rho,
                                                    const std::vector<T>& b, const T& Btarget) {
    return lossn_kr_dim(rho, b, std::vector<T>(rho.size(), Btarget));
}
 
/**
 * @brief Normal approximation of multirate single-link blocking for large C.
 *
 * In the critically loaded regime, where C and the offered loads grow together
 * with sum_k b_k rho_k / C -> 1, the normalized occupancy of a stochastic
 * knapsack converges to a truncated multivariate normal law (Kelly 1991; Ross
 * 1995, Theorem 2.8) and blocking is asymptotically proportional to the call
 * size (Ross 1995, Corollary 2.7). For a knapsack of fixed capacity C
 *
 *      alpha   = (C - sum_k b_k rho_k) / sqrt(C),
 *      sigma^2 = sum_k b_k^2 rho_k / C,
 *      delta   = phi(alpha/sigma) / (sigma Phi(alpha/sigma)),
 *      B_k     = b_k delta / sqrt(C),
 *
 * with phi and Phi the standard normal density and distribution function. B_k
 * is capped at 1. The approximation is accurate when |alpha| is small (roughly
 * below 1); far from critical loading use lossn_kr for the exact
 * values, or lossn_asy for the fluid (Kelly) limit.
 *
 * @param rho offered load of class k (K, nonnegative)
 * @param b   bandwidth units held by a class-k call (K, nonnegative)
 * @param C   link capacity (positive)
 */
template <class T>
KrAsyResult<T> lossn_kr_asy(const std::vector<T>& rho,
                                                    const std::vector<T>& b, const T& C) {
    static_assert(num_traits<T>::has_transcendental,
                  "lossn_kr_asy evaluates the normal density and distribution function");
    using std::sqrt;
    const char* fn = "lossn_kr_asy";
    const std::size_t K = rho.size();
    if (b.size() != K)
        throw InputError(std::string(fn) + ": b has " + std::to_string(b.size()) +
                         " entries but rho has " + std::to_string(K) + ".");
    const T zero = num_traits<T>::from_int(0);
    for (std::size_t k = 0; k < K; ++k)
        if (rho[k] < zero || b[k] < zero)
            throw InputError(std::string(fn) + ": rho and b must be nonnegative.");
    if (!(C > zero)) throw InputError(std::string(fn) + ": C must be a positive scalar.");
 
    T sbr = zero, sb2r = zero;
    for (std::size_t k = 0; k < K; ++k) {
        sbr += T(b[k] * rho[k]);
        sb2r += T(b[k] * b[k] * rho[k]);
    }
    const T sqC = sqrt(C);
    KrAsyResult<T> out;
    out.alpha = T(T(C - sbr) / sqC);
    out.sigma = sqrt(T(sb2r / C));
    out.B.assign(K, zero);
    if (out.sigma == zero) {
        // no class asks for any bandwidth: nothing is ever blocked
        out.delta = zero;
        return out;
    }
    out.delta = T(detail::kr_normal_hazard(T(out.alpha / out.sigma)) / out.sigma);
    for (std::size_t k = 0; k < K; ++k) out.B[k] = detail::kr_cap1(T(b[k] * out.delta / sqC));
    return out;
}
 
}  // namespace lossn
}  // namespace line
 
#endif  // LINE_API_LOSSN_LOSSN_KR_H