qsys_mmc_prio.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_QSYS_QSYS_MMC_PRIO_H
#define LINE_API_QSYS_QSYS_MMC_PRIO_H
 
/**
 * @file
 * @ingroup api_qsys
 * M/M/c with non-preemptive (head-of-line) priorities and a common exponential
 * service rate (Cobham, 1954).
 *
 * Port of matlab/src/api/qsys/qsys_mmc_prio.m. With A = sum(lambda)/mu, the
 * Erlang C probability of waiting C(c, A) and sigma_i = sum_{j<=i} lambda_j/(c mu),
 *
 *   Wq_i = (C / (c mu)) / ((1 - sigma_{i-1})(1 - sigma_i)),   W_i = Wq_i + 1/mu.
 *
 * ARITHMETIC: field. Erlang B is a rational recursion, so the computation is
 * exact for T = Rational.
 */
 
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/num/number.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** MATLAB's [W, rho, res] outputs. */
template <class T>
struct MmcPrioResult {
    std::vector<T> W; ///< per-class mean response time, class 1 first
    T rho;            ///< Q/(1+Q)
    QsysResT<T> res;  ///< meanSojournTime, meanWaitingTime, meanQueueLength, probWait, utilization, analyzer
};
 
/**
 * @brief M/M/c with head-of-line priorities.
 *
 * @param lambda per-class arrival rates, highest priority first
 * @param mu     common service rate
 * @param c      number of servers
 */
template <class T>
MmcPrioResult<T> qsys_mmc_prio(const std::vector<T>& lambda, const T& mu, unsigned c) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    for (const T& l : lambda)
        if (l <= zero)
            throw InputError(ErrorKind("Argument"), "qsys_mmc_prio: lambda must be positive and mu a positive scalar");
    if (mu <= zero)
        throw InputError(ErrorKind("Argument"), "qsys_mmc_prio: lambda must be positive and mu a positive scalar");
    if (c < 1) throw InputError(ErrorKind("Argument"), "qsys_mmc_prio: c must be a positive integer");
    const std::size_t R = lambda.size();
    const T cT = num_traits<T>::from_int(static_cast<long>(c));
    T lsum = zero;
    for (const T& l : lambda) lsum += l;
    const T A = T(lsum / mu);
    if (A >= cT) resdetail::throw_unstable("qsys_mmc_prio", "lambda/(c mu)", num_traits<T>::to_double(T(A / cT)));
    T B = one;
    for (unsigned k = 1; k <= c; ++k) B = T(A * B / (num_traits<T>::from_int(static_cast<long>(k)) + A * B));
    const T C = T(cT * B / (cT - A * (one - B)));
    std::vector<T> sig(R + 1, zero);
    for (std::size_t i = 0; i < R; ++i) sig[i + 1] = T(sig[i] + lambda[i]);
    for (T& s : sig) s = T(s / (cT * mu));
    MmcPrioResult<T> r;
    std::vector<T> Wq(R), Qi(R);
    r.W.assign(R, zero);
    T Q = zero;
    for (std::size_t i = 0; i < R; ++i) {
        Wq[i] = T((C / (cT * mu)) / ((one - sig[i]) * (one - sig[i + 1])));
        r.W[i] = T(Wq[i] + one / mu);
        Qi[i] = T(lambda[i] * r.W[i]);
        Q += Qi[i];
    }
    r.rho = T(Q / (one + Q));
    r.res.set("meanSojournTime", r.W);
    r.res.set("meanWaitingTime", Wq);
    r.res.set("meanQueueLength", Qi);
    r.res.set("probWait", C);
    r.res.set("utilization", T(A / cT));
    r.res.analyzer = "exact:M/M/" + std::to_string(c) + "-prio";
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMC_PRIO_H