lossn_revsens.h

Go to the documentation of this file.

#ifndef LINE_API_LOSSN_LOSSN_REVSENS_H
#define LINE_API_LOSSN_LOSSN_REVSENS_H
 
/**
 * @file lossn_revsens.h
 * @ingroup api_lossn
 * @brief Revenue sensitivities and implied costs of a loss network.
 *
 * Templated port of matlab/src/api/lossn/lossn_revsens.m.
 *
 * A route-r call in progress earns revenue at rate w_r, so the long-run
 * revenue rate is W = sum_r w_r E[n_r]. Under the product form the
 * sensitivity of W to the offered load of route k is (Kelly 1988; Ross 1995,
 * Sec. 6.3)
 *
 *      dW/dnu_k = (1 - B_k) (w_k - c_k),   c_k = W(C) - W(C - A e_k),
 *
 * where the implied cost c_k is the revenue lost by permanently removing the
 * units one route-k call holds. A positive dW_k means that admitting more
 * route-k traffic raises the revenue; a negative value means route k displaces
 * more valuable traffic than it earns. The link implied cost
 * clink(j,r) = W(C) - W(C - A(j,r) e_j) is the revenue lost by removing the
 * units of one route-r call from link j alone.
 *
 * Methods:
 *  - 'exact' (also 'manjunath', 'ms'): implied costs by one exact solution per
 *    route and per distinct (link, requirement) pair, needs integral A and C.
 *  - 'rec': as exact with MDD-rec, also for fractional A and C.
 *  - 'erlangfp': Kelly's approximation from the Erlang fixed point. The link
 *    implied costs solve the linear system
 *        c(j,k) = sum_l eta(j,l) [B_jl(C_j - A_jk) - B_jl(C_j)]
 *                        (w_l - sum_{i ~= j} c(i,l))
 *    with eta(j,l) the load of route l reduced by the other links and
 *    B_jl(.) the Kaufman-Roberts blocking of link j, and
 *    c_k = sum_j c(j,k), dW_k = prod_j (1 - E(j,k)) (w_k - c_k). It reproduces
 *    Ross Table 5.1.
 *  - 'mci': the score-function estimator dW/dnu_l = Cov(n_l, r(n))/nu_l of
 *    lossn_mci, with confidence intervals in ci.
 *
 * ARITHMETIC. The fixed point and the Monte Carlo estimator are
 * transcendental, so the header is gated on num_traits<T>::has_transcendental.
 *
 * References: F. P. Kelly, "Routing in circuit-switched networks: optimization,
 * shadow prices and decentralization", Adv. Appl. Probab. 20, 1988. K. W. Ross,
 * "Multiservice Loss Models for Broadband Telecommunication Networks",
 * Springer, 1995, Sec. 6.3 and Table 5.1.
 *
 * @see lossn_links, lossn_mci, lossn_erlangfp, lossn_kr
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/lossn/lossn_erlangfp.h"
#include "line/api/lossn/lossn_kr.h"
#include "line/api/lossn/lossn_links.h"
#include "line/api/lossn/lossn_manjunath.h"
#include "line/api/lossn/lossn_mci.h"
#include "line/api/lossn/lossn_rec.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace lossn {
 
/** Result of lossn_revsens. */
template <class T>
struct LossnRevsensResult {
    std::vector<T> dW;     ///< (R) dW/dnu_k
    std::vector<T> c;      ///< (R) route implied cost, NaN for a route that never fits
    T W = num_traits<T>::from_int(0);  ///< revenue rate sum_r w_r E[n_r]
    Matrix<T> clink;       ///< (J x R) link implied cost, NaN where A(j,r) > C(j)
    // Filled only by 'mci' (hasCI):
    bool hasCI = false;
    std::vector<T> revsens;  ///< (R) the point estimates of dW
    Matrix<T> revsensCI;     ///< (R x 2) interval on dW
    T level = num_traits<T>::from_int(0);  ///< confidence level of the interval
};
 
namespace detail {
 
/** The erlangfp branch: Kelly's implied-cost linear system at the Erlang fixed point. */
template <class T>
LossnRevsensResult<T> revsens_erlangfp(const std::vector<T>& nu, const Matrix<T>& A,
                                       const std::vector<T>& C, const std::vector<T>& w) {
    const std::size_t R = nu.size(), J = C.size();
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    // Same integral rescaling of each link as lossn_erlangfp, so that the
    // Kaufman-Roberts recursion counts whole units
    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 = erlangfp_rat_den(x);
            if (d > 0) m = m / 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));
    }
    const ErlangFpResult<T> fp = lossn_erlangfp(nu, A, C);
    const Matrix<T>& E = fp.E;
    LossnRevsensResult<T> out;
    for (std::size_t r = 0; r < R; ++r) out.W += T(w[r] * fp.QLen[r]);
 
    // Reduced load of route l on link j: thinned by every other link it uses
    Matrix<T> eta(J, R, zero);
    for (std::size_t j = 0; j < J; ++j)
        for (std::size_t l = 0; l < R; ++l)
            if (Aint[j][l] > 0) {
                T t = nu[l];
                for (std::size_t i = 0; i < J; ++i)
                    if (i != j && Aint[i][l] > 0) t = T(t * T(one - E(i, l)));
                eta(j, l) = t;
            }
 
    // Unknowns c(j,k) for every pair with Aint(j,k) > 0, in column-major order
    std::vector<std::size_t> pj, pk;
    std::vector<std::vector<long>> idx(J, std::vector<long>(R, -1));
    for (std::size_t k = 0; k < R; ++k)
        for (std::size_t j = 0; j < J; ++j)
            if (Aint[j][k] > 0) {
                idx[j][k] = static_cast<long>(pj.size());
                pj.push_back(j);
                pk.push_back(k);
            }
    const std::size_t P = pj.size();
    Matrix<T> M = eye<T>(P);
    std::vector<T> rhs(P, zero);
    for (std::size_t j = 0; j < J; ++j) {
        std::vector<T> etaj(R), bj(R);
        std::vector<long long> levels;
        for (std::size_t l = 0; l < R; ++l) {
            etaj[l] = eta(j, l);
            bj[l] = num_traits<T>::from_int(static_cast<long>(Aint[j][l]));
            if (Aint[j][l] > 0) levels.push_back(Aint[j][l]);
        }
        if (levels.empty()) continue;
        std::sort(levels.begin(), levels.end());
        levels.erase(std::unique(levels.begin(), levels.end()), levels.end());
        const std::vector<T> Q0 = lossn_kr(etaj, bj, static_cast<long>(Cint[j])).B;
        for (const long long a : levels) {
            const std::vector<T> Qa =
                Cint[j] - a < 0 ? std::vector<T>(R, one)
                                : lossn_kr(etaj, bj, static_cast<long>(Cint[j] - a)).B;
            for (std::size_t k = 0; k < R; ++k) {
                if (Aint[j][k] != a) continue;
                const std::size_t p = static_cast<std::size_t>(idx[j][k]);
                for (std::size_t l = 0; l < R; ++l) {
                    if (Aint[j][l] <= 0) continue;
                    const T coef = T(T(Qa[l] - Q0[l]) * eta(j, l));  // extra blocking of route l
                    rhs[p] += T(coef * w[l]);
                    for (std::size_t i = 0; i < J; ++i)
                        if (i != j && Aint[i][l] > 0) {
                            const std::size_t q = static_cast<std::size_t>(idx[i][l]);
                            M(p, q) = T(M(p, q) + coef);
                        }
                }
            }
        }
    }
    const std::vector<T> cvec = P > 0 ? solve(M, rhs) : std::vector<T>();
    out.clink = Matrix<T>(J, R, zero);
    for (std::size_t p = 0; p < P; ++p) out.clink(pj[p], pk[p]) = cvec[p];
    out.c.assign(R, zero);
    out.dW.assign(R, zero);
    for (std::size_t k = 0; k < R; ++k) {
        T pass = one;
        for (std::size_t j = 0; j < J; ++j) {
            out.c[k] += out.clink(j, k);
            pass = T(pass * T(one - E(j, k)));
        }
        out.dW[k] = T(pass * T(w[k] - out.c[k]));
    }
    return out;
}
 
}  // namespace detail
 
/**
 * @brief Revenue sensitivities and implied costs of a loss network.
 *
 * @param nu     offered load of route r (R)
 * @param A      capacity requirement of link j for route r (J x R)
 * @param C      capacity of link j (J)
 * @param w      revenue rate of a route-r call in progress (R)
 * @param method 'exact' (default), 'rec', 'erlangfp' or 'mci'
 * @param mciopt options for 'mci' (samples, alpha, gamma); its revenue field is
 *               overwritten with w
 * @param seed   seed of the 'mci' generator
 */
template <class T>
LossnRevsensResult<T> lossn_revsens(const std::vector<T>& nu, const Matrix<T>& A,
                                    const std::vector<T>& C, const std::vector<T>& w,
                                    const std::string& method = "exact",
                                    const LossnMciOptions<T>& mciopt = LossnMciOptions<T>(),
                                    std::uint64_t seed = 0) {
    static_assert(num_traits<T>::has_transcendental,
                  "lossn_revsens requires transcendental arithmetic");
    const std::size_t R = nu.size(), J = C.size();
    if (A.rows() != J || A.cols() != R)
        throw InputError("lossn_revsens: A must be " + std::to_string(J) + "x" + std::to_string(R) +
                         " (J x R).");
    if (w.size() != R)
        throw InputError("lossn_revsens: w must have " + std::to_string(R) + " entries.");
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const T nan = num_traits<T>::from_double(std::numeric_limits<double>::quiet_NaN());
    const std::string m = detail::lossn_lower(method.empty() ? std::string("exact") : method);
 
    if (m == "exact" || m == "manjunath" || m == "ms" || m == "rec") {
        const bool rec = (m == "rec");
        auto qlen = [&](const std::vector<T>& Cv) {
            return rec ? lossn_rec(nu, A, Cv).QLen : lossn_manjunath(nu, A, Cv).QLen;
        };
        auto revenue = [&](const std::vector<T>& Q) {
            T s = zero;
            for (std::size_t r = 0; r < R; ++r) s += T(w[r] * Q[r]);
            return s;
        };
        std::vector<T> QLen, Loss;
        if (rec) {
            const LossnRecResult<T> r0 = lossn_rec(nu, A, C);
            QLen = r0.QLen;
            Loss = r0.Loss;
        } else {
            const LossnManjunathResult<T> r0 = lossn_manjunath(nu, A, C);
            QLen = r0.QLen;
            Loss = r0.Loss;
        }
        LossnRevsensResult<T> out;
        out.W = revenue(QLen);
        out.c.assign(R, nan);
        out.dW.assign(R, zero);
        for (std::size_t k = 0; k < R; ++k) {
            std::vector<T> Ck(J);
            bool fits = true, allzero = true;
            for (std::size_t j = 0; j < J; ++j) {
                Ck[j] = T(C[j] - A(j, k));
                if (Ck[j] < zero) fits = false;
                if (A(j, k) != zero) allzero = false;
            }
            if (!fits) continue;  // never admitted: B_k = 1 and dW_k = 0
            out.c[k] = allzero ? zero : T(out.W - revenue(qlen(Ck)));
            out.dW[k] = T(T(one - Loss[k]) * T(w[k] - out.c[k]));
        }
        out.clink = Matrix<T>(J, R, zero);
        for (std::size_t j = 0; j < J; ++j)
            for (const T& a : detail::lossn_row_levels(A, j)) {
                T val = nan;
                if (!(T(C[j] - a) < zero)) {
                    std::vector<T> Cs = C;
                    Cs[j] = T(C[j] - a);
                    val = T(out.W - revenue(qlen(Cs)));
                }
                for (std::size_t r = 0; r < R; ++r)
                    if (A(j, r) == a) out.clink(j, r) = val;
            }
        return out;
    }
    if (m == "erlangfp") return detail::revsens_erlangfp(nu, A, C, w);
    if (m == "mci") {
        LossnMciOptions<T> o = mciopt;
        o.revenue = w;
        const LossnMciResult<T> cim = lossn_mci(nu, A, C, o, seed);
        LossnRevsensResult<T> out;
        out.W = cim.revenue;
        out.dW = cim.revsens;
        out.c = cim.impliedCost;
        out.clink = cim.linkImpliedCost;
        out.hasCI = true;
        out.revsens = cim.revsens;
        out.revsensCI = cim.revsensCI;
        out.level = cim.level;
        return out;
    }
    throw InputError("lossn_revsens: Unknown method '" + method +
                     "'; use exact, rec, erlangfp or mci.");
}
 
}  // namespace lossn
}  // namespace line
 
#endif  // LINE_API_LOSSN_LOSSN_REVSENS_H