qsys_mer1.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_MER1_H
#define LINE_API_QSYS_QSYS_MER1_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the M/Er/1 queue (Erlang-r service of total rate mu/r per
 * phase-sum, i.e. each phase has rate mu).
 *
 * Port of matlab/src/api/qsys/qsys_mer1.m. The number of phases of work in the
 * system has the distribution p_n = sum_k c_k x_k^n, where x_1..x_r are the
 * roots inside the unit disk of mu x^{r+1} - (lambda + mu) x^r + lambda and
 * c_k = (1 - rho) / prod_{j != k} (1 - x_j / x_k). The number of customers
 * groups the phases in blocks of r. The waiting time is a mixture of
 * exponentials and the sojourn time a mixture of Erlangs.
 *
 * The roots and coefficients are complex: they are returned beside `res`
 * (MATLAB's `res.roots` and `res.coeffs`), which holds the real measures.
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdio>
#include <functional>
#include <vector>
 
#include "line/api/qsys/qsys_mg1.h"
#include "line/api/qsys/qsys_mg1_ps.h"
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** The M/Er/1 `res`: the real measures plus the complex `roots` and `coeffs`. */
struct Mer1Result {
    QsysRes res;                             ///< every real field of MATLAB's res, in order
    std::vector<std::complex<double>> roots;  ///< res.roots, the x_k inside the unit disk
    std::vector<std::complex<double>> coeffs; ///< res.coeffs, the c_k
};
 
/**
 * @brief Exact M/Er/1 queue.
 *
 * @param lambda arrival rate
 * @param mu     rate of each Erlang phase (mean service r/mu)
 * @param r      number of phases
 * @param o      tPoints, percentiles and nmax (number of customer levels)
 */
inline Mer1Result qsys_mer1(double lambda, double mu, unsigned r, const QsysResOptions& o = QsysResOptions()) {
    using namespace resdetail;
    using Cplx = std::complex<double>;
    check_percentiles(o.percentiles, "qsys_mer1");
    if (!(lambda > 0.0) || !(mu > 0.0)) throw Error("qsys_mer1: lambda and mu must be positive scalars");
    if (r < 1) throw Error("qsys_mer1: r must be a positive integer");
    const double rho = lambda * r / mu;
    if (rho >= 1.0) throw_unstable("qsys_mer1", "rho", rho);
 
    // roots([mu, -(lambda + mu), zeros(1, r-1), lambda])
    std::vector<double> poly(r + 2, 0.0);
    poly[0] = mu;
    poly[1] = -(lambda + mu);
    poly[r + 1] = lambda;
    std::vector<Cplx> all = mg1psdetail::poly_roots(poly);
    // Newton polish against the original coefficients.
    for (Cplx& z : all) {
        for (int it = 0; it < 4; ++it) {
            Cplx f(0.0, 0.0), df(0.0, 0.0);
            for (double c : poly) {
                df = df * z + f;
                f = f * z + c;
            }
            if (std::abs(df) < 1e-300) break;
            z -= f / df;
        }
    }
    std::vector<Cplx> x;
    for (const Cplx& z : all)
        if (std::abs(z) < 1.0 - 1e-10) x.push_back(z);
    if (x.size() != r) {
        char buf[128];
        std::snprintf(buf, sizeof buf, "qsys_mer1: found %d roots inside the unit disk, expected %d",
                      static_cast<int>(x.size()), static_cast<int>(r));
        throw NumericError(ErrorKind{"Root"}, buf);
    }
    std::vector<Cplx> ck(r);
    for (unsigned k = 0; k < r; ++k) {
        Cplx prod(1.0, 0.0);
        for (unsigned j = 0; j < r; ++j)
            if (j != k) prod *= (1.0 - x[j] / x[k]);
        ck[k] = (1.0 - rho) / prod;
    }
    const double rr = static_cast<double>(r);
    const double ES = rr / mu, ES2 = rr * (rr + 1.0) / (mu * mu), ES3 = rr * (rr + 1.0) * (rr + 2.0) / (mu * mu * mu);
    Mer1Result out;
    out.res = qsys_mg1_res(lambda, 1.0 / ES, std::sqrt(ES2 / (ES * ES) - 1.0), ES3);
 
    double xmax = 0.0, ckmax = 0.0;
    for (const Cplx& z : x) xmax = std::max(xmax, std::abs(z));
    for (const Cplx& c : ck) ckmax = std::max(ckmax, std::abs(c));
    long nlev;
    if (o.nmax < 0) {
        long nph = static_cast<long>(std::ceil(std::log(1e-14 * (1.0 - xmax) / ckmax) / std::log(xmax))) + r;
        nph = std::max(nph, static_cast<long>(r));
        nlev = static_cast<long>(std::ceil(static_cast<double>(nph) / rr));
    } else {
        nlev = o.nmax;
    }
    const long nph = nlev * static_cast<long>(r);
    std::vector<double> p(static_cast<std::size_t>(nph + 1));
    for (long n = 0; n <= nph; ++n) {
        Cplx s(0.0, 0.0);
        for (unsigned k = 0; k < r; ++k) s += ck[k] * std::pow(x[k], static_cast<double>(n));
        p[static_cast<std::size_t>(n)] = std::max(s.real(), 0.0);
    }
    std::vector<double> q(static_cast<std::size_t>(nlev + 1), 0.0);
    q[0] = p[0];
    for (long i = 1; i <= nlev; ++i)
        for (long j = (i - 1) * r + 1; j <= i * static_cast<long>(r); ++j) q[static_cast<std::size_t>(i)] += p[static_cast<std::size_t>(j)];
 
    out.res.set("queueLengthDist", q);
    out.res.set("phaseDist", p);
    out.roots = x;
    out.coeffs = ck;
    const std::function<double(double)> Wccdf = [=](double t) {
        if (t < 0.0) return 1.0;
        Cplx s(0.0, 0.0);
        for (unsigned k = 0; k < r; ++k) s += ck[k] * x[k] / (1.0 - x[k]) * std::exp(-mu * (1.0 - x[k]) * t);
        return s.real();
    };
    const std::function<double(double)> Tccdf = [=](double t) {
        if (t < 0.0) return 1.0;
        double s = 0.0;
        for (long n = 0; n <= nph; ++n)
            s += p[static_cast<std::size_t>(n)] * gammainc_upper(mu * t, static_cast<double>(n) + rr);
        return s;
    };
    if (!o.tPoints.empty()) {
        out.res.set("tPoints", o.tPoints);
        out.res.set("waitingTimeCdf", cdf_of(Wccdf, o.tPoints));
        out.res.set("sojournTimeCdf", cdf_of(Tccdf, o.tPoints));
    }
    if (!o.percentiles.empty()) {
        const double hi = 50.0 / (mu * (1.0 - xmax));
        out.res.set("percentiles", o.percentiles);
        out.res.set("waitingTimePercentile", pct_of(Wccdf, o.percentiles, hi));
        out.res.set("sojournTimePercentile", pct_of(Tccdf, o.percentiles, hi));
    }
    out.res.analyzer = "exact:M/E" + std::to_string(r) + "/1";
    return out;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MER1_H