lossn_erlangfp.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_ERLANGFP_H
#define LINE_API_LOSSN_LOSSN_ERLANGFP_H
 
/**
 * @file
 * @ingroup api_lossn
 * Erlang fixed-point (reduced-load) approximation for a loss network.
 *
 * Templated port of matlab/src/api/lossn/lossn_erlangfp.m. Link j sees class
 * r offered at nu_r thinned by the blocking of r on every OTHER link,
 *   rho_(j,r) = nu_r prod_(i != j, A(i,r) > 0) (1 - B(i,r)),
 * and its occupancy distribution is the Kaufman-Roberts one (Kaufman 1981,
 * Roberts 1981), n q(n) = sum_r A(j,r) rho_(j,r) q(n - A(j,r)), so the
 * blocking of class r on link j is B(j,r) = P(more than C_j - A(j,r) units
 * busy). The (J x R) matrix B is the fixed point of that map, reached by the
 * shared iteration in line::da::da_fpi. With unit requirements this is Kelly's
 * Erlang fixed point; the earlier form fed the multi-unit load to Erlang B,
 * which counted a call of size a as one unit.
 *
 * A fractional row is the same constraint scaled to integers: row j is
 * multiplied by the least common denominator of A(j,:) and C_j (each entry
 * read as the rational of MATLAB's rat(x, 1e-12)), which leaves the admissible
 * set unchanged and gives the recursion whole units to count.
 *
 * The recursion itself is rational, but the fixed point around it stops on a
 * tolerance, so the function stays transcendental-gated: an exact
 * instantiation would promise more than the algorithm delivers.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <vector>
 
#include "line/api/da/da_fpi.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace lossn {
 
template <class T>
struct ErlangFpResult {
    std::vector<T> QLen;  ///< carried traffic per class
    std::vector<T> Loss;  ///< blocking probability per class
    Matrix<T> E;          ///< (J x R) blocking of class r on link j (the fixed point)
    std::size_t iterations = 0;
    bool converged = false;
};
 
/**
 * Erlang's loss formula B(nu, C), evaluated through logs as MATLAB does so the
 * factorials stay in range for large C.
 */
template <class T>
T erlang_b(const T& nu, int C) {
    static_assert(num_traits<T>::has_transcendental,
                  "erlang_b as ported evaluates through log/exp; use the rational recursion "
                  "B_k = nu B_{k-1} / (k + nu B_{k-1}) for an exact variant");
    if (C < 0) throw InputError("erlang_b: negative capacity");
    using std::exp;
    using std::log;
    const T lnu = log(nu);
    T den = num_traits<T>::from_int(0);
    for (int i = 0; i <= C; ++i)
        den += exp(num_traits<T>::from_int(i) * lnu - log(num_factorial<T>(static_cast<unsigned>(i))));
    const T lb = num_traits<T>::from_int(C) * lnu -
                 log(num_factorial<T>(static_cast<unsigned>(C))) - log(den);
    return exp(lb);
}
 
/**
 * Erlang B for a REAL capacity, the reference's ErlangB(nu, C) evaluated as it
 * stands: the sum runs over i = 0..floor(C) and factln(C) is lgamma(C+1).
 */
template <class T>
T erlang_b_real(const T& nu, const T& C) {
    using std::exp;
    using std::lgamma;
    using std::log;
    const double cd = num_traits<T>::to_double(C);
    if (!(cd >= 0.0)) return num_traits<T>::from_double(std::numeric_limits<double>::quiet_NaN());
    const T lnu = log(nu);
    T den = num_traits<T>::from_int(0);
    const long cfloor = static_cast<long>(std::floor(cd));
    for (long i = 0; i <= cfloor; ++i)
        den += exp(num_traits<T>::from_int(i) * lnu - log(num_factorial<T>(static_cast<unsigned>(i))));
    const T lb = C * lnu - num_traits<T>::from_double(lgamma(cd + 1.0)) - log(den);
    return exp(lb);
}
 
/**
 * The reference's fixed-point settings: iter_max 10000, iter_tol 1e-10, the
 * max-norm of the increment (`da_norm = @(d) max(abs(d(:)))`, the driver's
 * default) and `da_nanstop = true`.
 */
inline da::FpiOptions lossn_erlangfp_options() {
    da::FpiOptions o;
    o.iter_max = 10000;
    o.iter_tol = 1e-10;
    o.nanstop = true;
    return o;
}
 
namespace detail {
 
/**
 * Denominator of MATLAB's rat(x, tol): the continued-fraction convergent n/d
 * of x is extended until |x - n/d| <= tol.
 */
inline long long erlangfp_rat_den(double x, double tol = 1e-12) {
    if (!std::isfinite(x)) throw InputError("lossn_erlangfp: non-finite circuit requirement or capacity");
    double y = x;
    long long n0 = 1, d0 = 0, n1 = static_cast<long long>(std::floor(y)), d1 = 1;
    double frac = y - std::floor(y);
    for (int k = 0; k < 64 && std::fabs(x - static_cast<double>(n1) / static_cast<double>(d1)) > tol; ++k) {
        if (frac == 0.0) break;
        y = 1.0 / frac;
        const long long a = static_cast<long long>(std::floor(y));
        frac = y - std::floor(y);
        const long long n2 = a * n1 + n0, d2 = a * d1 + d0;
        n0 = n1; d0 = d1; n1 = n2; d1 = d2;
    }
    return d1 < 0 ? -d1 : d1;
}
 
inline long long erlangfp_gcd(long long a, long long b) {
    while (b != 0) { const long long t = a % b; a = b; b = t; }
    return a < 0 ? -a : a;
}
 
/**
 * Per-class blocking on one link of capacity c, class r asking a[r] units at
 * offered load rho[r]: the Kaufman-Roberts recursion n q(n) = sum_r a_r rho_r
 * q(n - a_r), normalized, and B_r = sum_{n > c - a_r} q(n).
 */
template <class T>
std::vector<T> kaufman_roberts_blocking(const std::vector<T>& rho, const std::vector<long long>& a, long long c) {
    const std::size_t R = rho.size();
    const T zero = num_traits<T>::from_int(0);
    std::vector<T> q(static_cast<std::size_t>(c + 1), zero);
    q[0] = num_traits<T>::from_int(1);
    for (long long n = 1; n <= c; ++n) {
        T s = zero;
        for (std::size_t r = 0; r < R; ++r)
            if (a[r] > 0 && n - a[r] >= 0)
                s += num_traits<T>::from_int(static_cast<int>(a[r])) * rho[r] *
                     q[static_cast<std::size_t>(n - a[r])];
        q[static_cast<std::size_t>(n)] = T(s / num_traits<T>::from_int(static_cast<int>(n)));
        // rescale; only the normalized distribution is used
        if (num_traits<T>::to_double(q[static_cast<std::size_t>(n)]) > 1e250) {
            const T scale = q[static_cast<std::size_t>(n)];
            for (long 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;
    std::vector<T> b(R, zero);
    for (std::size_t r = 0; r < R; ++r) {
        if (a[r] <= 0) continue;
        T acc = zero;
        for (long long n = std::max<long long>(0, c - a[r] + 1); n <= c; ++n) acc += q[static_cast<std::size_t>(n)];
        b[r] = T(acc / tot);
    }
    return b;
}
 
}  // namespace detail
 
/**
 * @brief Erlang fixed-point (reduced-load) approximation for a loss network.
 *
 * @param nu offered load per class (R)
 * @param A  (J x R) capacity requirement of link j for class r; fractional
 *           entries are scaled to integers row by row (see the file comment)
 * @param C  (J) link capacities
 * @param options fixed-point options; the default is the reference's
 *           (max-norm, iter_tol 1e-10, nanstop), see lossn_erlangfp_options()
 *
 * Returns the carried traffic QLen_r = nu_r prod_j (1 - B(j,r)), the loss
 * 1 - prod_j (1 - B(j,r)) and the (J x R) blocking matrix E = B.
 */
template <class T>
ErlangFpResult<T> lossn_erlangfp(const std::vector<T>& nu, const Matrix<T>& A,
                                 const std::vector<T>& C,
                                 const da::FpiOptions& options = lossn_erlangfp_options()) {
    static_assert(num_traits<T>::has_transcendental,
                  "lossn_erlangfp requires transcendental arithmetic (a tolerance-driven fixed "
                  "point around the Kaufman-Roberts recursion)");
    const std::size_t R = nu.size();
    const std::size_t J = C.size();
    if (A.rows() != J || A.cols() != R)
        throw InputError("lossn_erlangfp: the route matrix does not match nu and C");
 
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
 
    // Integer rows: multiply row j by the lcm of the denominators of A(j,:), C_j.
    std::vector<std::vector<long long>> Aint(J, std::vector<long long>(R, 0));
    std::vector<long long> Cint(J, 0);
    for (std::size_t j = 0; j < J; ++j) {
        long long m = 1;
        for (std::size_t r = 0; r <= R; ++r) {
            const double x = (r < R) ? num_traits<T>::to_double(A(j, r)) : num_traits<T>::to_double(C[j]);
            const long long d = detail::erlangfp_rat_den(x);
            if (d > 0) m = m / detail::erlangfp_gcd(m, d) * d;
        }
        for (std::size_t r = 0; r < R; ++r)
            Aint[j][r] = std::llround(static_cast<double>(m) * num_traits<T>::to_double(A(j, r)));
        Cint[j] = static_cast<long long>(std::floor(static_cast<double>(m) * num_traits<T>::to_double(C[j]) + 1e-9));
        if (Cint[j] < 0) throw InputError("lossn_erlangfp: negative link capacity");
    }
 
    // One sweep of the reduced-load map over the row-major (J x R) iterate.
    const std::function<std::pair<std::vector<T>, std::vector<T>>(const std::vector<T>&, std::size_t)>
        sweep = [&](const std::vector<T>& B, std::size_t) {
            std::vector<T> Bnew = B;
            for (std::size_t j = 0; j < J; ++j) {
                std::vector<T> rho = nu;
                for (std::size_t r = 0; r < R; ++r)
                    for (std::size_t i = 0; i < J; ++i)
                        if (i != j && Aint[i][r] > 0) rho[r] = T(rho[r] * T(one - B[i * R + r]));
                const std::vector<T> bj = detail::kaufman_roberts_blocking<T>(rho, Aint[j], Cint[j]);
                for (std::size_t r = 0; r < R; ++r) Bnew[j * R + r] = bj[r];
            }
            return std::make_pair(Bnew, B);
        };
 
    const std::vector<T> B0(J * R, zero);
    da::FpiResult<T> fp = da::da_fpi<T>(sweep, B0, options);
 
    ErlangFpResult<T> r;
    r.E = Matrix<T>(J, R, zero);
    for (std::size_t j = 0; j < J; ++j)
        for (std::size_t c = 0; c < R; ++c) r.E(j, c) = fp.x[j * R + c];
    r.iterations = fp.iterations;
    r.converged = fp.converged;
    r.QLen.assign(R, zero);
    r.Loss.assign(R, zero);
    for (std::size_t c = 0; c < R; ++c) {
        T pass = one;
        for (std::size_t j = 0; j < J; ++j) pass = T(pass * T(one - r.E(j, c)));
        r.QLen[c] = T(nu[c] * pass);
        r.Loss[c] = T(one - pass);  // blocking probability = 1 - carried/offered
    }
    return r;
}
 
/** lossn_erlangfp with integer capacities. */
template <class T>
ErlangFpResult<T> lossn_erlangfp(const std::vector<T>& nu, const Matrix<T>& A,
                                 const std::vector<int>& C,
                                 const da::FpiOptions& options = lossn_erlangfp_options()) {
    std::vector<T> Ct(C.size());
    for (std::size_t j = 0; j < C.size(); ++j) Ct[j] = num_traits<T>::from_int(C[j]);
    return lossn_erlangfp(nu, A, Ct, options);
}
 
}  // namespace lossn
}  // namespace line
 
#endif  // LINE_API_LOSSN_LOSSN_ERLANGFP_H