pfqn_sens_momldmx.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_SENS_MOMLDMX_H
#define LINE_API_PFQN_SENS_MOMLDMX_H
 
/**
 * @file
 * @ingroup api_pfqn
 * Exact higher queue-length moments (up to order three) of mixed open/closed
 * product-form networks with limited load dependence.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_sens_momldmx.m. Theorem 1 of
 * Akyildiz and Strelen (IEEE TComm 39(6):828-832, 1991) states that one
 * further factor Q_(j,T) in a moment costs one derivative with respect to a
 * parameter y_j scaling the demands of the classes in T at station j. With
 * m = E[Q_(j,T)]:
 *
 *   Var[Q_jT]       = dm/dy_j
 *   E[Q_jT^3]       = d2m/dy_j^2 + (1 + 3m) dm/dy_j + m^3
 *   Cov[Q_iT,Q_jT'] = dm_iT/dy_(j,T')
 *
 * The derivatives come from second-order forward-mode differentiation of the
 * Bruell-Balbo-Afshari mixed load-dependent MVA (the recursion of
 * pfqn_mvaldmx and pfqn_sens_mvaldmx), with the effective-capacity second
 * derivatives supplied by pfqn_sens_ldmx_ec. The result struct is the one of
 * pfqn_sens_mom, so the two are interchangeable on closed load-independent
 * input.
 *
 * Arithmetic. Field-only apart from the skewness, as in pfqn_sens_mom; the
 * primal Pc(0) keeps the reference's eps floor.
 */
 
#include <cmath>
#include <cstddef>
#include <vector>
 
#include "line/api/pfqn/pfqn_sens_ldmx_ec.h"
#include "line/api/pfqn/pfqn_sens_mom.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
#include "line/util/population.h"
 
namespace line {
namespace pfqn {
 
/**
 * @brief Exact moments E[Q], Var[Q], E[Q^2], E[Q^3] and covariances of the
 *        grouped queue lengths of a mixed load-dependent product-form network.
 *
 * @param lambda (R) arrival rates, zero on closed classes
 * @param D      (M x R) service demands
 * @param N      (R) closed populations; a negative entry marks an open class
 * @param Z      (R) think times
 * @param mu     (M x max(1,sum N closed)) limited load-dependent rates
 * @param groups (R) 0-based group label of each class; empty means one group
 *               holding every class (per-station totals)
 */
template <class T>
SensMomResult<T> pfqn_sens_momldmx(const std::vector<T>& lambda, const Matrix<T>& D,
                                   const std::vector<int>& N, const std::vector<T>& Z,
                                   const Matrix<T>& mu, const std::vector<int>& groups) {
    const std::size_t M = D.rows();
    const std::size_t R = D.cols();
    if (lambda.size() != R || N.size() != R || Z.size() != R)
        throw InputError("pfqn_sens_momldmx: lambda, N, Z and D disagree on the class count");
    if (mu.rows() != M)
        throw InputError("pfqn_sens_momldmx: mu and D disagree on the station count");
 
    std::vector<int> grp = groups.empty() ? std::vector<int>(R, 0) : groups;
    if (grp.size() != R)
        throw InputError("pfqn_sens_momldmx: groups must have one label per class");
    std::size_t G = 0;
    for (int g : grp) {
        if (g < 0) throw InputError("pfqn_sens_momldmx: group labels start at zero");
        if (static_cast<std::size_t>(g) + 1 > G) G = static_cast<std::size_t>(g) + 1;
    }
    {
        std::vector<bool> seen(G, false);
        for (int g : grp) seen[static_cast<std::size_t>(g)] = true;
        for (std::size_t g = 0; g < G; ++g)
            if (!seen[g])
                throw InputError(
                    "pfqn_sens_momldmx: groups must label the classes consecutively with no empty group");
    }
 
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
    const T two = num_traits<T>::from_int(2);
 
    std::vector<std::size_t> openClasses, closedClasses;
    for (std::size_t r = 0; r < R; ++r) {
        if (N[r] < 0) {
            openClasses.push_back(r);
        } else {
            if (lambda[r] != zero && N[r] > 0)
                throw InputError("pfqn_sens_momldmx: an arrival rate cannot be set on a closed class");
            closedClasses.push_back(r);
        }
    }
    const std::size_t Cn = closedClasses.size();
 
    std::vector<int> Nc(Cn, 0);
    std::vector<T> Zc(Cn, zero);
    Matrix<T> Dc(M, Cn, zero);
    int NCtot = 0;
    for (std::size_t c = 0; c < Cn; ++c) {
        Nc[c] = N[closedClasses[c]];
        Zc[c] = Z[closedClasses[c]];
        NCtot += Nc[c];
        for (std::size_t i = 0; i < M; ++i) Dc(i, c) = D(i, closedClasses[c]);
    }
    if (static_cast<int>(mu.cols()) < (NCtot > 1 ? NCtot : 1))
        throw InputError(
            "pfqn_sens_momldmx: the load-dependent rates must be given up to the maximum closed "
            "population");
 
    Matrix<T> mux(M, mu.cols() + 1, zero);
    for (std::size_t i = 0; i < M; ++i) {
        for (std::size_t k = 0; k < mu.cols(); ++k) mux(i, k) = mu(i, k);
        mux(i, mu.cols()) = mu(i, mu.cols() - 1);
    }
    const SensLdmxEcResult<T> ec = pfqn_sens_ldmx_ec(lambda, D, mux);
 
    // ---- parameter list: y(h,g) multiplies D(h,r) for every r in group g ----
    const std::size_t P = M * G;
    const auto pidx = [&](std::size_t h, std::size_t g) { return h * G + g; };
    // eq. (21): Lo is linear in the open demands, so d2Lo = 0.
    Matrix<T> dLo(M, P, zero);
    for (std::size_t j = 0; j < M; ++j)
        for (std::size_t oi = 0; oi < openClasses.size(); ++oi) {
            const std::size_t o = openClasses[oi];
            dLo(j, pidx(j, static_cast<std::size_t>(grp[o]))) += lambda[o] * D(j, o);
        }
 
    // ---- population lattice over the closed classes, first class fastest ----
    const std::vector<std::size_t> prods = plane_sizes(Nc);
    const std::size_t NT = population_count(Nc);
    const std::size_t NL = static_cast<std::size_t>(NCtot) + 1;
 
    std::vector<Matrix<T>> Pc(NT, Matrix<T>(M, NL, zero));
    std::vector<std::vector<Matrix<T>>> dPc(NT, std::vector<Matrix<T>>(M, Matrix<T>(NL, P, zero)));
    std::vector<std::vector<Matrix<T>>> d2Pc(NT, std::vector<Matrix<T>>(M, Matrix<T>(NL, P, zero)));
    Matrix<T> x(NT, Cn, zero);
    std::vector<Matrix<T>> dx(NT, Matrix<T>(Cn, P, zero));
    std::vector<Matrix<T>> d2x(NT, Matrix<T>(Cn, P, zero));
    std::vector<Matrix<T>> w(NT, Matrix<T>(M, Cn, zero));
    std::vector<std::vector<Matrix<T>>> dw(NT, std::vector<Matrix<T>>(M, Matrix<T>(Cn, P, zero)));
    std::vector<std::vector<Matrix<T>>> d2w(NT, std::vector<Matrix<T>>(M, Matrix<T>(Cn, P, zero)));
 
    for (std::size_t i = 0; i < M; ++i) Pc[0](i, 0) = one;  // eq. (16)
 
    std::vector<T> dacc(P, zero), d2acc(P, zero);
    for (std::size_t k = 0; k < NT; ++k) {
        std::vector<int> nvec(Cn, 0);
        int nc = 0;
        for (std::size_t c = 0; c < Cn; ++c) {
            nvec[c] = static_cast<int>((k / prods[c]) % static_cast<std::size_t>(Nc[c] + 1));
            nc += nvec[c];
        }
        if (nc == 0) continue;  // the empty population keeps Pc(0) = 1
 
        // ---- residence times, eq. (12)-(13) ------------------------------------
        for (std::size_t i = 0; i < M; ++i) {
            for (std::size_t c = 0; c < Cn; ++c) {
                if (nvec[c] <= 0) continue;
                const std::size_t kc = k - prods[c];
                const std::size_t pc = pidx(i, static_cast<std::size_t>(grp[closedClasses[c]]));
                T acc = zero;
                for (std::size_t p = 0; p < P; ++p) {
                    dacc[p] = zero;
                    d2acc[p] = zero;
                }
                for (int n = 1; n <= nc; ++n) {
                    const std::size_t nu = static_cast<std::size_t>(n);
                    const T Pp = Pc[kc](i, nu - 1);
                    const T nT = num_traits<T>::from_int(n);
                    const T ECn = ec.EC(i, nu - 1);
                    const T dECn = ec.dEC(i, nu - 1);
                    const T d2ECn = ec.d2EC(i, nu - 1);
                    acc += nT * ECn * Pp;
                    for (std::size_t p = 0; p < P; ++p) {
                        const T dL = dLo(i, p);
                        const T dPp = dPc[kc][i](nu - 1, p);
                        dacc[p] += nT * (dECn * dL * Pp + ECn * dPp);
                        d2acc[p] += nT * (d2ECn * dL * dL * Pp + two * dECn * dL * dPp +
                                          ECn * d2Pc[kc][i](nu - 1, p));
                    }
                }
                w[k](i, c) = Dc(i, c) * acc;
                for (std::size_t p = 0; p < P; ++p) {
                    if (p == pc) {  // d(D y)/dy = D
                        dw[k][i](c, p) = Dc(i, c) * (dacc[p] + acc);
                        d2w[k][i](c, p) = Dc(i, c) * (d2acc[p] + two * dacc[p]);
                    } else {
                        dw[k][i](c, p) = Dc(i, c) * dacc[p];
                        d2w[k][i](c, p) = Dc(i, c) * d2acc[p];
                    }
                }
            }
        }
 
        // ---- throughputs, eq. (18) ---------------------------------------------
        for (std::size_t c = 0; c < Cn; ++c) {
            T sw = zero;
            for (std::size_t i = 0; i < M; ++i) sw += w[k](i, c);
            const T den = Zc[c] + sw;
            if (den == zero) {
                x(k, c) = zero;
                continue;
            }
            const T nT = num_traits<T>::from_int(nvec[c]);
            x(k, c) = nT / den;
            if (nvec[c] > 0) {
                const T den2 = den * den;
                const T den3 = den2 * den;
                for (std::size_t p = 0; p < P; ++p) {
                    T sdw = zero, sd2w = zero;
                    for (std::size_t i = 0; i < M; ++i) {
                        sdw += dw[k][i](c, p);
                        sd2w += d2w[k][i](c, p);
                    }
                    dx[k](c, p) = -nT * sdw / den2;
                    d2x[k](c, p) = -nT * sd2w / den2 + two * nT * sdw * sdw / den3;
                }
            }
        }
 
        // ---- conditional marginals, eq. (14)-(17) --------------------------------
        for (std::size_t i = 0; i < M; ++i) {
            for (int n = 1; n <= nc; ++n) {
                const std::size_t nu = static_cast<std::size_t>(n);
                for (std::size_t c = 0; c < Cn; ++c) {
                    if (nvec[c] <= 0) continue;
                    const std::size_t kc = k - prods[c];
                    const std::size_t pc = pidx(i, static_cast<std::size_t>(grp[closedClasses[c]]));
                    const T Pp = Pc[kc](i, nu - 1);
                    const T ECn = ec.EC(i, nu - 1);
                    const T dECn = ec.dEC(i, nu - 1);
                    const T d2ECn = ec.d2EC(i, nu - 1);
                    const T xc = x(k, c);
                    const T t = ECn * xc * Pp;
                    Pc[k](i, nu) += Dc(i, c) * t;
                    for (std::size_t p = 0; p < P; ++p) {
                        const T dL = dLo(i, p);
                        const T dE1 = dECn * dL;
                        const T d2E1 = d2ECn * dL * dL;
                        const T dxc = dx[k](c, p);
                        const T dPp = dPc[kc][i](nu - 1, p);
                        const T dt = dE1 * xc * Pp + ECn * dxc * Pp + ECn * xc * dPp;
                        const T d2t = d2E1 * xc * Pp + ECn * d2x[k](c, p) * Pp +
                                      ECn * xc * d2Pc[kc][i](nu - 1, p) +
                                      two * (dE1 * dxc * Pp + dE1 * xc * dPp + ECn * dxc * dPp);
                        if (p == pc) {
                            dPc[k][i](nu, p) += Dc(i, c) * (dt + t);
                            d2Pc[k][i](nu, p) += Dc(i, c) * (d2t + two * dt);
                        } else {
                            dPc[k][i](nu, p) += Dc(i, c) * dt;
                            d2Pc[k][i](nu, p) += Dc(i, c) * d2t;
                        }
                    }
                }
            }
            // eq. (17), with the reference's floor on the primal only
            T s1 = zero;
            for (int n = 1; n <= nc; ++n) s1 += Pc[k](i, static_cast<std::size_t>(n));
            const T epsT = num_traits<T>::from_double(2.220446049250313e-16);
            const T p0 = one - s1;
            Pc[k](i, 0) = p0 > epsT ? p0 : epsT;
            for (std::size_t p = 0; p < P; ++p) {
                T ds = zero, d2s = zero;
                for (int n = 1; n <= nc; ++n) {
                    ds += dPc[k][i](static_cast<std::size_t>(n), p);
                    d2s += d2Pc[k][i](static_cast<std::size_t>(n), p);
                }
                dPc[k][i](0, p) = -ds;
                d2Pc[k][i](0, p) = -d2s;
            }
        }
    }
 
    // ---- measures at the full population -------------------------------------
    const std::size_t kN = NT - 1;
    SensMomResult<T> res;
    res.XN.assign(R, zero);
    res.QN = Matrix<T>(M, R, zero);
    res.UN = Matrix<T>(M, R, zero);
    res.CN = Matrix<T>(M, R, zero);
    std::vector<Matrix<T>> dQN(P, Matrix<T>(M, R, zero));
    std::vector<Matrix<T>> d2QN(P, Matrix<T>(M, R, zero));
 
    // closed classes, eq. (32)
    for (std::size_t c = 0; c < Cn; ++c) {
        const std::size_t cls = closedClasses[c];
        res.XN[cls] = x(kN, c);
        const std::size_t kc = Nc[c] > 0 ? kN - prods[c] : kN;
        for (std::size_t i = 0; i < M; ++i) {
            res.CN(i, cls) = w[kN](i, c);
            res.QN(i, cls) = res.XN[cls] * res.CN(i, cls);
            for (std::size_t p = 0; p < P; ++p) {
                dQN[p](i, cls) = dx[kN](c, p) * w[kN](i, c) + x(kN, c) * dw[kN][i](c, p);
                d2QN[p](i, cls) = d2x[kN](c, p) * w[kN](i, c) + two * dx[kN](c, p) * dw[kN][i](c, p) +
                                  x(kN, c) * d2w[kN][i](c, p);
            }
            if (Nc[c] > 0) {
                T uacc = zero;
                for (int n = 1; n <= NCtot; ++n) {
                    const std::size_t nu = static_cast<std::size_t>(n);
                    uacc += Dc(i, c) * x(kN, c) * ec.Eprime(i, nu - 1) / ec.E(i, nu - 1) *
                            Pc[kc](i, nu - 1);
                }
                res.UN(i, cls) = uacc;
            }
        }
    }
 
    // open classes, eq. (33)
    for (std::size_t oi = 0; oi < openClasses.size(); ++oi) {
        const std::size_t r = openClasses[oi];
        res.XN[r] = lambda[r];
        for (std::size_t i = 0; i < M; ++i) {
            const std::size_t pr = pidx(i, static_cast<std::size_t>(grp[r]));
            T acc = zero;
            for (std::size_t p = 0; p < P; ++p) {
                dacc[p] = zero;
                d2acc[p] = zero;
            }
            for (int n = 0; n <= NCtot; ++n) {
                const std::size_t nu = static_cast<std::size_t>(n);
                const T Pn = Pc[kN](i, nu);
                const T nT = num_traits<T>::from_int(n + 1);
                const T ECn = ec.EC(i, nu);
                const T dECn = ec.dEC(i, nu);
                const T d2ECn = ec.d2EC(i, nu);
                acc += nT * ECn * Pn;
                for (std::size_t p = 0; p < P; ++p) {
                    const T dL = dLo(i, p);
                    const T dPn = dPc[kN][i](nu, p);
                    dacc[p] += nT * (dECn * dL * Pn + ECn * dPn);
                    d2acc[p] += nT * (d2ECn * dL * dL * Pn + two * dECn * dL * dPn +
                                      ECn * d2Pc[kN][i](nu, p));
                }
            }
            const T scale = lambda[r] * D(i, r);
            res.QN(i, r) = scale * acc;
            if (lambda[r] != zero) res.CN(i, r) = res.QN(i, r) / lambda[r];
            for (std::size_t p = 0; p < P; ++p) {
                if (p == pr) {
                    dQN[p](i, r) = scale * (dacc[p] + acc);
                    d2QN[p](i, r) = scale * (d2acc[p] + two * dacc[p]);
                } else {
                    dQN[p](i, r) = scale * dacc[p];
                    d2QN[p](i, r) = scale * d2acc[p];
                }
            }
            T uacc = zero;
            for (int n = 0; n <= NCtot; ++n) {
                const std::size_t nu = static_cast<std::size_t>(n);
                uacc += lambda[r] * ec.Eprime(i, nu + 1) / ec.E(i, nu + 1) * Pc[kN](i, nu);
            }
            res.UN(i, r) = uacc;
        }
    }
 
    // ---- group moments, packed exactly as pfqn_sens_mom ------------------------
    res.m = Matrix<T>(M, G, zero);
    res.d2m = Matrix<T>(M, G, zero);
    res.Cov = Matrix<T>(P, P, zero);
    res.dm = Matrix<T>(P, P, zero);
    res.Var = Matrix<T>(M, G, zero);
    res.M2 = Matrix<T>(M, G, zero);
    res.M3 = Matrix<T>(M, G, zero);
    res.Skew = Matrix<double>(M, G, 0.0);
    res.CovAsym = zero;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < R; ++r) {
            const std::size_t g = static_cast<std::size_t>(grp[r]);
            res.m(i, g) += res.QN(i, r);
            for (std::size_t j = 0; j < M; ++j)
                for (std::size_t g2 = 0; g2 < G; ++g2)
                    res.dm(pidx(i, g), pidx(j, g2)) += dQN[pidx(j, g2)](i, r);
            res.d2m(i, g) += d2QN[pidx(i, g)](i, r);
        }
    for (std::size_t a = 0; a < P; ++a)
        for (std::size_t b = 0; b < P; ++b) {
            const T d = num_abs(T(res.dm(a, b) - res.dm(b, a)));
            if (d > res.CovAsym) res.CovAsym = d;
        }
    for (std::size_t a = 0; a < P; ++a)
        for (std::size_t b = 0; b < P; ++b) res.Cov(a, b) = (res.dm(a, b) + res.dm(b, a)) / two;
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t g = 0; g < G; ++g) {
            const T d1 = res.dm(pidx(i, g), pidx(i, g));
            const T mg = res.m(i, g);
            res.Var(i, g) = d1;
            res.M2(i, g) = d1 + mg * mg;
            res.M3(i, g) = res.d2m(i, g) + (one + num_traits<T>::from_int(3) * mg) * d1 + mg * mg * mg;
            const T mu3 = res.M3(i, g) - num_traits<T>::from_int(3) * mg * res.M2(i, g) + two * mg * mg * mg;
            const double var = num_traits<T>::to_double(res.Var(i, g));
            res.Skew(i, g) = var > 0.0 ? num_traits<T>::to_double(mu3) / std::pow(var, 1.5)
                                       : std::nan("");
        }
    return res;
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_SENS_MOMLDMX_H