qsys_mg1_prio_pr.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_MG1_PRIO_PR_H
#define LINE_API_QSYS_QSYS_MG1_PRIO_PR_H
 
/**
 * @file
 * @ingroup api_qsys
 * M/G/1 with preemptive-resume priorities: per-class mean (and optionally
 * variance of the) response time.
 *
 * Port of matlab/src/api/qsys/qsys_mg1_prio_pr.m. Class 1 is the highest
 * priority. With sigma_i = sum_{j<=i} rho_j and M2_i = sum_{j<=i} lambda_j E[S_j^2],
 *
 *   W_i = E[S_i]/(1 - sigma_{i-1}) + M2_i / (2 (1 - sigma_{i-1})(1 - sigma_i)),
 *
 * and, given third moments, the variance of Takagi (1991) Sec. 3.5. The
 * returned rho is the qsys family's Q/(1+Q) with Q = sum_i lambda_i W_i.
 *
 * ARITHMETIC: field; exact for T = Rational.
 */
 
#include <cstddef>
#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 Mg1PrioPrResult {
    std::vector<T> W; ///< per-class mean response time, class 1 first
    T rho;            ///< Q/(1+Q)
    QsysResT<T> res;  ///< meanSojournTime, meanWaitingTime, meanQueueLength, utilization[, varSojournTime], analyzer
};
 
/**
 * @brief M/G/1 with preemptive-resume priorities.
 *
 * @param lambda per-class arrival rates, highest priority first
 * @param mu     per-class service rates
 * @param cs     per-class coefficients of variation of the service time
 * @param ES3    per-class third moments of the service time (MATLAB's 'ES3'
 *               option); empty when not given
 */
template <class T>
Mg1PrioPrResult<T> qsys_mg1_prio_pr(const std::vector<T>& lambda, const std::vector<T>& mu,
                                    const std::vector<T>& cs, const std::vector<T>& ES3 = std::vector<T>()) {
    const std::size_t R = lambda.size();
    if (mu.size() != R || cs.size() != R)
        throw InputError(ErrorKind("Argument"), "qsys_mg1_prio_pr: lambda, mu and cs must have the same length");
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const T two = num_traits<T>::from_int(2), three = num_traits<T>::from_int(3), four = num_traits<T>::from_int(4);
    for (std::size_t i = 0; i < R; ++i)
        if (lambda[i] <= zero || mu[i] <= zero || cs[i] < zero)
            throw InputError(ErrorKind("Argument"),
                             "qsys_mg1_prio_pr: lambda and mu must be positive and cs non-negative");
    std::vector<T> ES(R), ES2(R), rho_i(R);
    T rhoSum = zero;
    for (std::size_t i = 0; i < R; ++i) {
        ES[i] = T(one / mu[i]);
        ES2[i] = T((one + cs[i] * cs[i]) / (mu[i] * mu[i]));
        rho_i[i] = T(lambda[i] * ES[i]);
        rhoSum += rho_i[i];
    }
    if (rhoSum >= one) resdetail::throw_unstable("qsys_mg1_prio_pr", "rho", num_traits<T>::to_double(rhoSum));
    std::vector<T> sig(R + 1, zero), M2(R + 1, zero);  // sig[i] = sigma_{i-1} in 1-based terms
    for (std::size_t i = 0; i < R; ++i) {
        sig[i + 1] = T(sig[i] + rho_i[i]);
        M2[i + 1] = T(M2[i] + lambda[i] * ES2[i]);
    }
    Mg1PrioPrResult<T> r;
    r.W.assign(R, zero);
    T Q = zero;
    for (std::size_t i = 0; i < R; ++i) {
        r.W[i] = T(ES[i] / (one - sig[i]) + M2[i + 1] / (two * (one - sig[i]) * (one - sig[i + 1])));
        Q += T(lambda[i] * r.W[i]);
    }
    r.rho = T(Q / (one + Q));
    std::vector<T> Wq(R), Qi(R);
    for (std::size_t i = 0; i < R; ++i) {
        Wq[i] = T(r.W[i] - ES[i]);
        Qi[i] = T(lambda[i] * r.W[i]);
    }
    r.res.set("meanSojournTime", r.W);
    r.res.set("meanWaitingTime", Wq);
    r.res.set("meanQueueLength", Qi);
    r.res.set("utilization", rhoSum);
    if (!ES3.empty()) {
        if (ES3.size() != R)
            throw InputError(ErrorKind("Argument"), "qsys_mg1_prio_pr: ES3 must have one entry per class");
        std::vector<T> M3(R + 1, zero);
        for (std::size_t i = 0; i < R; ++i) M3[i + 1] = T(M3[i] + lambda[i] * ES3[i]);
        std::vector<T> V(R);
        for (std::size_t i = 0; i < R; ++i) {
            const T a = T(one - sig[i]), b = T(one - sig[i + 1]);
            V[i] = T((ES2[i] - ES[i] * ES[i]) / (a * a) + ES[i] * M2[i] / (a * a * a) +
                     M3[i + 1] / (three * a * a * b) + M2[i + 1] * M2[i + 1] / (four * a * a * b * b) +
                     M2[i + 1] * M2[i] / (two * a * a * a * b));
        }
        r.res.set("varSojournTime", V);
    }
    r.res.analyzer = "exact:M/G/1-prio-pr";
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_PRIO_PR_H