lossn_asy.h

Go to the documentation of this file.

#ifndef LINE_API_LOSSN_LOSSN_ASY_H
#define LINE_API_LOSSN_LOSSN_ASY_H
 
/**
 * @file lossn_asy.h
 * @ingroup api_lossn
 * @brief Large-capacity (Kelly limit) approximation of loss networks.
 *
 * Templated port of matlab/src/api/lossn/lossn_asy.m.
 *
 * Scales offered loads nu and capacities C together, nu -> N nu and C -> N C,
 * and returns the limit of the class loss probabilities as N grows (Kelly
 * 1986; Ross 1995, Theorem 5.13). In this limit link j blocks independently
 * with probability a_j, and the a_j solve the complementary slackness
 * conditions
 *
 *   sum_k A(j,k) nu_k prod_i (1-a_i)^A(i,k) = C_j  if a_j > 0,
 *                                           <= C_j  if a_j = 0,
 *
 * the route loss being 1 - prod_j (1-a_j)^A(j,k). With y_j = -log(1-a_j)
 * these are the optimality conditions of the strictly convex problem
 *
 *   minimize  sum_k nu_k exp(-sum_j A(j,k) y_j) + sum_j C_j y_j,  y >= 0,
 *
 * solved here by cyclic coordinate descent, each coordinate by Newton's
 * method on a concave increasing function. The result needs no summation over
 * states and costs O(J R) per sweep, so it scales to any capacity; it reports
 * zero loss on links that are underloaded at the fluid scale, so it is
 * accurate in heavy (overloaded) traffic and for large C. A and C need not be
 * integral.
 *
 * ARITHMETIC. Tolerance-driven Newton steps on exponentials: gated on
 * num_traits<T>::has_transcendental.
 *
 * References: F. P. Kelly, "Blocking probabilities in large circuit-switched
 * networks", Adv. Appl. Probab. 18, 1986. F. P. Kelly, "Loss Networks", Ann.
 * Appl. Probab. 1(3), 1991. K. W. Ross, "Multiservice Loss Models for
 * Broadband Telecommunication Networks", Springer, 1995, Sec. 5.8.
 *
 * @see lossn_erlangfp, lossn_links
 */
 
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace lossn {
 
/** Result of lossn_asy. */
template <class T>
struct LossnAsyResult {
    std::vector<T> QLen;  ///< mean number of route-r calls in progress, nu_r (1 - Loss_r) (R)
    std::vector<T> Loss;  ///< limiting loss probability of route r (R)
    std::vector<T> a;     ///< limiting blocking probability of link j (J)
    int niter = 0;        ///< number of coordinate-descent sweeps
};
 
/**
 * @brief Fluid (Kelly) limit of a loss network by convex duality.
 *
 * @param nu offered load of route (class) r (R)
 * @param A  capacity requirement of link j for route r (J x R, nonnegative)
 * @param C  capacity of link j (J, nonnegative)
 */
template <class T>
LossnAsyResult<T> lossn_asy(const std::vector<T>& nu, const Matrix<T>& A, const std::vector<T>& C) {
    static_assert(num_traits<T>::has_transcendental,
                  "lossn_asy requires transcendental arithmetic (Newton steps on exponentials)");
    using std::abs;
    using std::exp;
    using std::expm1;
    const std::size_t R = nu.size();
    const std::size_t J = C.size();
    if (A.rows() != J || A.cols() != R)
        throw InputError("lossn_asy: A must be " + std::to_string(J) + "x" + std::to_string(R) +
                         " (J x R).");
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    for (std::size_t k = 0; k < R; ++k)
        if (nu[k] < zero) throw InputError("lossn_asy: nu, A and C must be nonnegative.");
    for (std::size_t j = 0; j < J; ++j) {
        if (C[j] < zero) throw InputError("lossn_asy: nu, A and C must be nonnegative.");
        for (std::size_t k = 0; k < R; ++k)
            if (A(j, k) < zero) throw InputError("lossn_asy: nu, A and C must be nonnegative.");
    }
 
    const T tol = num_traits<T>::from_double(1e-12);
    const T steptol = num_traits<T>::from_double(1e-14);
    const int maxsweeps = 10000;
 
    // A link of zero capacity blocks every route through it with certainty;
    // those routes carry nothing and drop out of the optimization
    LossnAsyResult<T> out;
    out.a.assign(J, zero);
    std::vector<char> dead(J, 0), lost(R, 0);
    for (std::size_t j = 0; j < J; ++j) {
        if (!(C[j] <= zero)) continue;
        for (std::size_t k = 0; k < R; ++k)
            if (A(j, k) > zero) dead[j] = 1;
        if (dead[j]) {
            out.a[j] = one;
            for (std::size_t k = 0; k < R; ++k)
                if (A(j, k) > zero) lost[k] = 1;
        }
    }
    std::vector<T> nul = nu;
    for (std::size_t k = 0; k < R; ++k)
        if (lost[k]) nul[k] = zero;
 
    std::vector<T> y(J, zero);
    std::vector<T> s(R, zero);  // s_k = sum_j A(j,k) y_j over live links
    std::vector<std::size_t> users;
    std::vector<T> Aj, rest, base, e;
    for (int sweep = 1; sweep <= maxsweeps; ++sweep) {
        out.niter = sweep;
        T change = zero;
        for (std::size_t j = 0; j < J; ++j) {
            if (dead[j]) continue;
            users.clear();
            for (std::size_t k = 0; k < R; ++k)
                if (A(j, k) > zero && nul[k] > zero) users.push_back(k);
            if (users.empty()) continue;
            const std::size_t U = users.size();
            Aj.assign(U, zero);
            rest.assign(U, zero);
            base.assign(U, zero);
            e.assign(U, zero);
            T sab = zero;
            for (std::size_t u = 0; u < U; ++u) {
                Aj[u] = A(j, users[u]);
                rest[u] = T(s[users[u]] - Aj[u] * y[j]);  // contribution of the other links
                base[u] = T(nul[users[u]] * exp(T(-rest[u])));
                sab += T(Aj[u] * base[u]);
            }
            // g(t) = C_j - sum_k A_jk base_k exp(-A_jk t) is concave increasing
            T t = zero;
            if (!(T(C[j] - sab) >= zero)) {
                t = y[j];
                T g0 = zero;
                for (std::size_t u = 0; u < U; ++u) g0 += T(Aj[u] * base[u] * exp(T(-Aj[u] * t)));
                if (T(C[j] - g0) > zero) t = zero;  // restart left of the root
                for (int it = 1; it <= 100; ++it) {
                    T g = C[j], dg = zero;
                    for (std::size_t u = 0; u < U; ++u) {
                        e[u] = T(base[u] * exp(T(-Aj[u] * t)));
                        g -= T(Aj[u] * e[u]);
                        dg += T(Aj[u] * Aj[u] * e[u]);
                    }
                    const T step = T(-g / dg);
                    t += step;
                    const T scale = t > one ? t : one;
                    if (abs(step) <= T(steptol * scale)) break;
                }
            }
            const T anew = T(-expm1(T(-t)));
            const T d = abs(T(anew + expm1(T(-y[j]))));
            if (d > change) change = d;
            for (std::size_t u = 0; u < U; ++u) s[users[u]] = T(rest[u] + Aj[u] * t);
            y[j] = t;
        }
        if (change <= tol) break;
    }
    for (std::size_t j = 0; j < J; ++j)
        if (!dead[j]) out.a[j] = T(-expm1(T(-y[j])));
 
    out.Loss.assign(R, zero);
    out.QLen.assign(R, zero);
    for (std::size_t k = 0; k < R; ++k) {
        if (lost[k]) {
            out.Loss[k] = one;
        } else {
            T ay = zero;
            for (std::size_t j = 0; j < J; ++j)
                if (!dead[j]) ay += T(A(j, k) * y[j]);
            out.Loss[k] = T(-expm1(T(-ay)));
        }
        out.QLen[k] = T(nu[k] * T(one - out.Loss[k]));
    }
    return out;
}
 
}  // namespace lossn
}  // namespace line
 
#endif  // LINE_API_LOSSN_LOSSN_ASY_H