pfqn_qdamva.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_QDAMVA_H
#define LINE_API_PFQN_QDAMVA_H
 
/**
 * @file
 * @ingroup api_pfqn
 * QD-AMVA: queue-dependent approximate mean value analysis.
 *
 * Port of matlab/src/api/pfqn/pfqn_qdamva.m, the queue-dependent AMVA of
 * Casale, Perez and Wang (IFIP PERFORMANCE 2015), on a closed multiclass
 * product-form network.
 *
 * A Schweitzer/Bard core in which the class-r demand at station k is scaled by
 * the queue-dependence term g_k evaluated at the ARRIVAL-INSTANT total queue
 * length, `g = pfqn_lldfun(1 + delta * rowsum(Q), mu)`.
 *
 * SETTING `mu` TO A CONSTANT ROW RECOVERS PLAIN SCHWEITZER AMVA ONLY FOR A
 * SINGLE CLASS. `pfqn_lldfun` does skip a row of ONES, so g == 1 there, but the
 * residence time that remains is `1 + delta * rowsum(Q)` with ONE aggregate
 * delta = (sum(N)-1)/sum(N) applied to the whole arrival-instant queue, where
 * Bard-Schweitzer shrinks the TAGGED class alone:
 *
 *   1 + sum_{s != r} Q(k,s) + (N(r)-1)/N(r) * Q(k,r)
 *
 * The two coincide iff K == 1. Measured over 40 random three-class instances,
 * pfqn_qdamva(L,N,Z,ones) departs from pfqn_bs by up to 0.217 in absolute queue
 * length, and is the LESS accurate of the two on single-server multiclass models
 * (mean relative error on Q 0.069 against 0.056 at R = 3), the aggregate delta
 * buying nothing once g == 1. This is the QD-AMVA closure, not a defect of the
 * port, but do not use the function as a Schweitzer oracle for K > 1.
 *
 * MU IS A DIMENSIONLESS RATE MULTIPLIER, NOT A RATE. mu(k,n) is the factor by
 * which station k serves faster when it holds n jobs. Two traps follow from
 * `pfqn_lldfun` and are the reference's, not this port's:
 *
 *   - it SKIPS a station whose mu row is identically ONE, so a single-server
 *     station must be a row of ones and a c-server station `min(1..smax, c)`.
 *     Until 2026-09-13 the gate was `range(...) > 0` and skipped EVERY constant
 *     row, so a uniform multiplier c != 1 silently returned g = 1 -- the station
 *     ran unscaled, against what solver_mvald and SolverCTMC give.
 *   - smax = mu.cols() must be at least ceil(sum(N)) or the interpolation
 *     clamps the population and the top of the rate curve is never reached.
 *
 * Delay stations are carried in Z, not as rows of L. Closed classes only: an
 * infinite N(r) is not supported.
 *
 * Arithmetic: TRANSCENDENTAL-GATED, inherited whole from `pfqn_lldfun` -- the
 * softmin it evaluates is not an element of the field generated by the inputs.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <vector>
 
#include "line/api/pfqn/pfqn_lldfun.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
#include "line/num/number.h"
 
namespace line {
namespace pfqn {
 
/** What `pfqn_qdamva` returns: the fixed point and how it was reached. */
template <class T>
struct QdAmvaResult {
    Matrix<T> Q;  ///< (M x R) mean queue lengths
    Matrix<T> X;  ///< (1 x R) per-class throughputs
    Matrix<T> U;  ///< (M x R) per-class utilizations, carrying the g scaling
    Matrix<T> R;  ///< (M x R) per-class residence times, Q = X .* R
    std::size_t iter = 0;
};
 
/**
 * @brief QD-AMVA: queue-dependent approximate mean value analysis.
 *
 * @param L (M x R) service demand matrix.
 * @param N (R) population vector, finite.
 * @param Z (R) think time vector; empty means no think time.
 * @param mu (M x smax) queue-dependent rate multipliers; empty means none.
 * @param Q0 (M x R) initial guess; empty means the reference's demand split.
 * @param tol convergence tolerance on the queue lengths (default 1e-6).
 * @param maxiter maximum number of iterations (default 10000).
 */
template <class T>
QdAmvaResult<T> pfqn_qdamva(const Matrix<T>& L, const std::vector<T>& N, const std::vector<T>& Z,
                            const Matrix<T>& mu, const Matrix<T>& Q0, double tol = 1e-6,
                            std::size_t maxiter = 10000) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_qdamva requires transcendental arithmetic");
    const std::size_t M = static_cast<std::size_t>(L.rows());
    const std::size_t K = static_cast<std::size_t>(L.cols());
    if (N.size() != K)
        throw InputError("pfqn_qdamva: the population vector must have one entry per class");
    if (!Z.empty() && Z.size() != K)
        throw InputError("pfqn_qdamva: the think-time vector must have one entry per class");
 
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
 
    QdAmvaResult<T> out;
    out.Q = Matrix<T>(M, K, zero);
    out.X = Matrix<T>(1, K, zero);
    out.U = Matrix<T>(M, K, zero);
    out.R = Matrix<T>(M, K, zero);
 
    T Ntot = zero;
    for (std::size_t r = 0; r < K; ++r) Ntot += N[r];
    // delta is undefined on an empty population, and the reference returns the
    // zero queue rather than dividing by it.
    if (!(num_traits<T>::to_double(Ntot) > 0.0)) return out;
 
    if (Q0.rows() == static_cast<int>(M) && Q0.cols() == static_cast<int>(K)) {
        out.Q = Q0;
    } else {
        // Ltot = 0 for a class with no demand anywhere: L/Ltot is a NaN the
        // iteration never recovers from. Such a column arises routinely in a
        // layered fixed point, where a caller can start with no work at the
        // layer station, so the column is left at zero instead.
        for (std::size_t r = 0; r < K; ++r) {
            T tot = zero;
            for (std::size_t k = 0; k < M; ++k) tot += L(k, r);
            if (!(num_traits<T>::to_double(tot) > 0.0)) continue;
            for (std::size_t k = 0; k < M; ++k) out.Q(k, r) = T(L(k, r) / tot * N[r]);
        }
    }
 
    const T delta = T((Ntot - one) / Ntot);
    // Q*10 as the sentinel, as the reference notes, stalls on an all-zero seed:
    // the loop would exit before its first pass. Offset instead.
    Matrix<T> Qprev = out.Q;
    for (std::size_t k = 0; k < M; ++k)
        for (std::size_t r = 0; r < K; ++r)
            Qprev(k, r) = T(out.Q(k, r) + num_traits<T>::from_double(10.0 * (1.0 + tol)));
 
    std::vector<T> Ak(M, zero);
    while (out.iter < maxiter) {
        double gap = 0.0;
        for (std::size_t k = 0; k < M; ++k)
            for (std::size_t r = 0; r < K; ++r)
                gap = std::max(gap, std::fabs(num_traits<T>::to_double(out.Q(k, r)) -
                                              num_traits<T>::to_double(Qprev(k, r))));
        if (out.iter > 0 && gap <= tol) break;
        ++out.iter;
        Qprev = out.Q;
 
        // The arrival-instant total queue length, class independent.
        for (std::size_t k = 0; k < M; ++k) {
            T rowsum = zero;
            for (std::size_t r = 0; r < K; ++r) rowsum += out.Q(k, r);
            Ak[k] = T(one + delta * rowsum);
        }
        const std::vector<T> g = pfqn_lldfun(Ak, mu, std::vector<double>());
 
        for (std::size_t r = 0; r < K; ++r) {
            for (std::size_t k = 0; k < M; ++k) {
                T rowsum = zero;
                for (std::size_t s = 0; s < K; ++s) rowsum += out.Q(k, s);
                out.R(k, r) = T(L(k, r) * g[k] * (one + delta * rowsum));
            }
            T denom = Z.empty() ? zero : Z[r];
            for (std::size_t k = 0; k < M; ++k) denom += out.R(k, r);
            out.X(0, r) = (num_traits<T>::to_double(denom) > 0.0) ? T(N[r] / denom) : zero;
            for (std::size_t k = 0; k < M; ++k) {
                out.Q(k, r) = T(out.X(0, r) * out.R(k, r));
                out.U(k, r) = T(L(k, r) * g[k] * out.X(0, r));
            }
        }
    }
    return out;
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_QDAMVA_H