pfqn_debemva.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_PFQN_DEBEMVA_H
#define LINE_API_PFQN_PFQN_DEBEMVA_H
 
/**
 * @file
 * @ingroup api_pfqn
 * debeMVA: MVA-based approximation of closed networks with finite (BAS) buffers.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_debemva.m, i.e. debeMVA (Algorithms
 * 1 and 2) of V. de Nitto Persone and A. Di Lonardo, "Approximating finite
 * resources: An approach based on MVA", Performance Evaluation 131:1-21, 2019.
 *
 * Closed single-class network of M stations. A queue is a single exponential
 * server with total capacity K[i]; a delay is an infinite server of any
 * distribution, never full. Blocking is Blocking-After-Service with
 * first-blocked-first-unblocked release.
 *
 * MVA on k = 1..N with residence time S*(z + n(k-1)) + BT at a queue and
 * S + BT/k at a delay. The arrival job function z (start 1) and the blocking
 * time BT (start 0) persist across populations. Algorithm 2 first adjusts them
 * around the slowest finite stations once k > K[j]+1; each pass then moves the
 * share 1-K[j]/k of z[j] upstream for every station with n[j] > K[j]+tol,
 * charges the blocking time P(i,j)*S[j]^2/(S[i]+S[j]), and restores any queue
 * that fell below its value at k-1 by taking the difference from its upstream
 * stations. tol = 0 is the paper's strict test. A capacity K[i] >= N is
 * treated as unbounded.
 *
 * Results follow the paper, so a blocked job is counted at the SOURCE, whereas
 * LINE's CTMC/JMT/LDES count it at the full destination.
 *
 * ARITHMETIC. Only field operations and comparisons, so every T is admitted,
 * Rational included. Capacities are doubles because an unbounded one is +inf.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <string>
#include <utility>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
/** Return value of pfqn_debemva, mirroring [X,Q,U,R,BT,z,theta,iter,converged]; vectors are (M). */
template <class T>
struct DebemvaResult {
    std::vector<T> X;   ///< throughput of each station
    std::vector<T> Q;   ///< mean number of jobs, blocked jobs at the source
    std::vector<T> U;   ///< utilization X.*S, blocked time excluded
    std::vector<T> R;   ///< mean residence time per visit, blocking time included
    std::vector<T> BT;  ///< blocking time accumulated at each station
    std::vector<T> z;   ///< arrival job function at each station
    long theta = 0;     ///< number of capacity violations corrected (the paper's theta)
    long iter = 0;      ///< total number of passes over all populations
    bool converged = true;  ///< false if a population exhausted iter_max passes
};
 
/**
 * Runs debeMVA.
 *
 * @param N        population (positive integer)
 * @param S        (M) mean service time per visit
 * @param P        (M x M) station routing matrix, row-stochastic and irreducible
 * @param K        (M) capacity, +inf for unbounded; ignored at delays
 * @param isdelay  (M) true at infinite-server stations, empty for none
 * @param tol      slack on the capacity test n[j] > K[j]+tol (the paper: 0)
 * @param iter_max maximum passes at one population
 */
template <class T>
DebemvaResult<T> pfqn_debemva(long N, const std::vector<T>& S, const Matrix<T>& P,
                              const std::vector<double>& K, std::vector<bool> isdelay,
                              double tol = 0.0, long iter_max = 100000) {
    using nt = num_traits<T>;
    const std::size_t M = S.size();
    const T zero = nt::from_int(0);
    const T one = nt::from_int(1);
    const double inf = std::numeric_limits<double>::infinity();
    if (isdelay.empty()) isdelay.assign(M, false);
    const std::string caller = "pfqn_debemva: ";
    if (N <= 0)
        throw InputError(ErrorKind{"Argument"}, caller + "debeMVA requires a positive integer population N.");
    if (P.rows() != M || P.cols() != M || K.size() != M || isdelay.size() != M)
        throw InputError(ErrorKind{"Argument"},
            caller + "debeMVA requires S, K and isdelay of length M and an M x M routing matrix P.");
    std::vector<double> B(M, inf);
    for (std::size_t i = 0; i < M; ++i) {
        const double s = nt::to_double(S[i]);
        if (!(s >= 0.0) || !std::isfinite(s))
            throw InputError(ErrorKind{"Argument"}, caller + "debeMVA requires finite nonnegative service times.");
        if (!isdelay[i] && !(K[i] >= 1.0))
            throw InputError(ErrorKind{"Argument"},
                caller + "debeMVA requires every queue capacity K(i) to be at least 1.");
        B[i] = (isdelay[i] || K[i] >= static_cast<double>(N)) ? inf : K[i];
        double row = 0.0;
        for (std::size_t j = 0; j < M; ++j) {
            const double p = nt::to_double(P(i, j));
            if (p < 0.0)
                throw InputError(ErrorKind{"Argument"}, caller + "debeMVA requires a row-stochastic routing matrix P.");
            row += p;
        }
        if (std::fabs(row - 1.0) > 1e-8)
            throw InputError(ErrorKind{"Argument"}, caller + "debeMVA requires a row-stochastic routing matrix P.");
    }
    std::vector<std::vector<bool>> reach(M, std::vector<bool>(M, false));
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j) reach[i][j] = i == j || P(i, j) > zero;
    const int squarings = static_cast<int>(std::ceil(std::log2(static_cast<double>(M < 2 ? 2 : M)))) + 1;
    for (int h = 0; h < squarings; ++h) {
        std::vector<std::vector<bool>> next(M, std::vector<bool>(M, false));
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t j = 0; j < M; ++j) {
                bool r = reach[i][j];
                for (std::size_t l = 0; l < M && !r; ++l) r = reach[i][l] && reach[l][j];
                next[i][j] = r;
            }
        reach.swap(next);
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j)
            if (!reach[i][j])
                throw InputError(ErrorKind{"Argument"},
                    caller + "debeMVA requires a routing matrix that is irreducible on the stations.");
    // minimum total capacity of a directed cycle of finite stations (Floyd-Warshall)
    std::vector<std::vector<double>> W(M, std::vector<double>(M, inf));
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j)
            if (i != j && P(i, j) > zero && std::isfinite(B[j])) W[i][j] = B[j];
    for (std::size_t h = 0; h < M; ++h)
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t j = 0; j < M; ++j)
                if (W[i][h] + W[h][j] < W[i][j]) W[i][j] = W[i][h] + W[h][j];
    double cyc = inf;
    for (std::size_t i = 0; i < M; ++i)
        if (W[i][i] < cyc) cyc = W[i][i];
    if (static_cast<double>(N) >= cyc)
        throw UnsupportedError(ErrorKind{"Structure"},
            caller + "debeMVA requires N below the total capacity " + std::to_string(cyc) +
            " of every cycle of finite stations, to rule out a deadlock of the blocking network; N = " +
            std::to_string(N) + ".");
 
    // visit ratios e = e*P, normalized at station 1; every output is invariant to the scale.
    // A = P' - I with row 1 replaced by e_1, solved by Gaussian elimination with partial pivoting.
    Matrix<T> A(M, M, zero);
    std::vector<T> e(M, zero);
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t j = 0; j < M; ++j) A(i, j) = T(P(j, i) - (i == j ? one : zero));
    for (std::size_t j = 0; j < M; ++j) A(0, j) = j == 0 ? one : zero;
    e[0] = one;
    for (std::size_t c = 0; c < M; ++c) {
        std::size_t piv = c;
        for (std::size_t r = c + 1; r < M; ++r)
            if (std::fabs(nt::to_double(A(r, c))) > std::fabs(nt::to_double(A(piv, c)))) piv = r;
        if (piv != c) {
            for (std::size_t j = 0; j < M; ++j) std::swap(A(c, j), A(piv, j));
            std::swap(e[c], e[piv]);
        }
        for (std::size_t r = 0; r < M; ++r) {
            if (r == c || A(r, c) == zero) continue;
            const T f = T(A(r, c) / A(c, c));
            for (std::size_t j = c; j < M; ++j) A(r, j) = T(A(r, j) - f * A(c, j));
            e[r] = T(e[r] - f * e[c]);
        }
    }
    for (std::size_t i = 0; i < M; ++i) e[i] = T(e[i] / A(i, i));
 
    // capacities as T, valid where B is finite
    std::vector<T> Bt(M, zero);
    std::vector<bool> fin(M, false);
    std::vector<int> nd(M, 0);
    bool anyFin = false;
    T smax = zero;
    for (std::size_t i = 0; i < M; ++i) {
        fin[i] = std::isfinite(B[i]);
        if (fin[i]) {
            Bt[i] = nt::from_double(B[i]);
            if (!anyFin || S[i] > smax) smax = S[i];
            anyFin = true;
        }
        for (std::size_t l = 0; l < M; ++l)
            if (P(i, l) > zero) ++nd[i];
    }
    std::vector<bool> slow(M, false);
    for (std::size_t i = 0; i < M; ++i) slow[i] = fin[i] && S[i] == smax;
    const T tolT = nt::from_double(tol);
 
    DebemvaResult<T> res;
    std::vector<T> z(M, one), bt(M, zero), n(M, zero), t(S), X(M, zero), beta(M, zero);
    std::vector<bool> first(M, false);
    for (long k = 1; k <= N; ++k) {
        const std::vector<T> nprev = n;
        const T kT = nt::from_int(k);
        for (std::size_t j = 0; j < M; ++j)
            beta[j] = (fin[j] && B[j] < static_cast<double>(k)) ? T(one - Bt[j] / kT) : zero;
        // Algorithm 2: initialization around the slowest finite stations
        for (std::size_t j = 0; j < M; ++j) {
            if (slow[j] && static_cast<double>(k) > B[j] + 1.0 && !first[j]) {
                first[j] = true;
                for (std::size_t i = 0; i < M; ++i)
                    if (i != j && P(i, j) > zero) {
                        z[i] = T(z[i] + P(i, j) * z[j] * beta[j]);
                        bt[i] = T(bt[i] + P(i, j) * S[j] * T(S[j] / T(S[j] + S[i])));
                    }
                z[j] = T(z[j] - z[j] * beta[j]);
            }
        }
        for (std::size_t j = 0; j < M; ++j) {
            if (!(slow[j] && static_cast<double>(k) > B[j] + 1.0)) continue;
            for (std::size_t i = 0; i < M; ++i) {
                if (i == j || !(P(i, j) > zero) || nd[i] <= 1) continue;
                for (std::size_t l = 0; l < M; ++l)
                    if (l != j && P(i, l) > zero) z[l] = T(z[l] + z[i] * beta[j] * P(i, l));
                z[i] = T(z[i] - z[i] * beta[j]);
            }
        }
        // Algorithm 1: passes until no capacity is exceeded
        long passes = 0;
        while (true) {
            ++passes;
            ++res.iter;
            T denom = zero;
            for (std::size_t i = 0; i < M; ++i) {
                t[i] = isdelay[i] ? T(S[i] + bt[i] / kT) : T(S[i] * (z[i] + nprev[i]) + bt[i]);
                denom = T(denom + e[i] * t[i]);
            }
            for (std::size_t i = 0; i < M; ++i) {
                X[i] = T(kT * e[i] / denom);
                n[i] = T(t[i] * X[i]);
            }
            bool rep = false;
            for (std::size_t j = 0; j < M; ++j) {
                if (fin[j] && n[j] > T(Bt[j] + tolT)) {
                    rep = true;
                    ++res.theta;
                    for (std::size_t i = 0; i < M; ++i)
                        if (i != j && P(i, j) > zero && (isdelay[i] || !fin[i] || n[i] < Bt[i])) {
                            z[i] = T(z[i] + P(i, j) * z[j] * beta[j]);
                            bt[i] = T(bt[i] + P(i, j) * S[j] * T(S[j] / T(S[j] + S[i])));
                        }
                    z[j] = T(z[j] - z[j] * beta[j]);
                } else if (n[j] < nprev[j]) {
                    const T d = T(nprev[j] - n[j]);
                    n[j] = nprev[j];
                    T den = zero;
                    for (std::size_t i = 0; i < M; ++i)
                        if (i != j && P(i, j) > zero) den = T(den + P(i, j));
                    if (den > zero)
                        for (std::size_t i = 0; i < M; ++i)
                            if (i != j && P(i, j) > zero) n[i] = T(n[i] - P(i, j) * d / den);
                }
            }
            if (!rep) break;
            if (passes >= iter_max) {
                res.converged = false;
                break;
            }
        }
    }
    res.X = X;
    res.Q = n;
    res.R = t;
    res.BT = bt;
    res.z = z;
    res.U.assign(M, zero);
    for (std::size_t i = 0; i < M; ++i) res.U[i] = T(X[i] * S[i]);
    return res;
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_PFQN_DEBEMVA_H