qsys_gimc.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_GIMC_H
#define LINE_API_QSYS_QSYS_GIMC_H
 
/**
 * @file
 * @ingroup api_qsys
 * Exact analysis of the GI/M/c queue (Takacs, 1962).
 *
 * Port of matlab/src/api/qsys/qsys_gimc.m. With A(s) the interarrival LST,
 * omega in (0,1) is the root of omega = A(c mu (1 - omega)). An arrival finds
 * n >= c-1 jobs with probability D omega^(n-c), and fewer with the alternating
 * binomial sums of the U_n built from g_j = A(j mu); D normalizes. Waiting is
 * then exponential, P(W > t) = PW exp(-c mu (1 - omega) t), PW = D/(1-omega),
 * and the time-average queue law follows from the arrival law by rate
 * balance. MATLAB's refusals when omega cannot be bracketed (LINE:numerical:Root)
 * and when the alternating sums lose accuracy (LINE:numerical:Cancellation)
 * are kept.
 *
 * ARITHMETIC: double, as the reference. MATLAB's fzero on [0, 1-d] is a
 * bracketing root search; it is reproduced here by bisection run to the
 * resolution of the double grid, which brackets the same simple root.
 */
 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/lang/distribution.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/**
 * @brief Exact GI/M/c queue, interarrival law given by its LST.
 *
 * @param lst    interarrival LST at real arguments
 * @param lambda arrival rate (MATLAB's 'lambda' option, required with a handle)
 * @param mu     service rate of one server
 * @param c      number of servers
 * @param o      tPoints, percentiles and nmax options
 */
inline QsysRes qsys_gimc(const std::function<double(double)>& lst, double lambda, double mu, unsigned c,
                         const QsysResOptions& o = QsysResOptions()) {
    using namespace resdetail;
    check_percentiles(o.percentiles, "qsys_gimc");
    if (!lst) throw Error("qsys_gimc: the interarrival law must be a Distribution or an LST function handle");
    if (!(lambda > 0.0) || !std::isfinite(lambda))
        throw Error("qsys_gimc: with an LST handle the arrival rate must be given with 'lambda'");
    if (!(mu > 0.0)) throw Error("qsys_gimc: mu must be a positive scalar");
    if (c < 1) throw Error("qsys_gimc: c must be a positive integer");
    const double A = lambda / mu;
    if (A >= c) throw_unstable("qsys_gimc", "lambda/(c mu)", A / c);
 
    const std::function<double(double)> g = [&](double w) { return lst(c * mu * (1.0 - w)) - w; };
    double d = 1e-3;
    while (g(1.0 - d) >= 0.0 && d > 1e-15) d /= 10.0;
    if (g(1.0 - d) >= 0.0) throw NumericError(ErrorKind{"Root"}, "qsys_gimc: could not bracket the root omega in (0,1)");
    double lo = 0.0, hi = 1.0 - d;
    for (int it = 0; it < 2000; ++it) {
        const double mid = 0.5 * (lo + hi);
        if (mid <= lo || mid >= hi) break;
        const double gm = g(mid);
        if (gm == 0.0) {
            lo = hi = mid;
            break;
        }
        if (gm > 0.0) lo = mid;
        else hi = mid;
    }
    const double omega = std::fabs(g(lo)) <= std::fabs(g(hi)) ? lo : hi;
 
    std::vector<double> gj(c), Cj(c), terms(c);
    for (unsigned j = 1; j <= c; ++j) gj[j - 1] = lst(j * mu);
    double cp = 1.0;
    for (unsigned j = 1; j <= c; ++j) {
        cp *= gj[j - 1] / (1.0 - gj[j - 1]);
        Cj[j - 1] = cp;
    }
    double tsum = 0.0;
    for (unsigned j = 1; j <= c; ++j) {
        terms[j - 1] = nchoosek_r(c, j) * (c * (1.0 - gj[j - 1]) - j) /
                       (Cj[j - 1] * (1.0 - gj[j - 1]) * (c * (1.0 - omega) - j));
        tsum += terms[j - 1];
    }
    const double D = 1.0 / (1.0 / (1.0 - omega) + tsum);
    std::vector<double> Cn(c + 1);
    Cn[0] = 1.0;
    for (unsigned j = 1; j <= c; ++j) Cn[j] = Cj[j - 1];
    std::vector<double> Un(c);
    for (unsigned nn = 0; nn < c; ++nn) {
        double s = 0.0;
        for (unsigned j = nn; j < c; ++j) s += terms[j];
        Un[nn] = D * Cn[nn] * s;
    }
    long nmax = o.nmax;
    if (nmax < 0)
        nmax = static_cast<long>(c) +
               static_cast<long>(std::max(1.0, std::ceil(std::log(1e-14 * (1.0 - omega)) /
                                                         std::log(std::max(omega, std::numeric_limits<double>::min())))));
    nmax = std::max<long>(nmax, c);
    std::vector<double> Pi(static_cast<std::size_t>(nmax + 1));
    for (long nn = 0; nn <= nmax; ++nn) {
        if (nn >= static_cast<long>(c) - 1) {
            Pi[nn] = D * std::pow(omega, static_cast<double>(nn - static_cast<long>(c)));
        } else {
            double s = 0.0;
            for (long ii = nn; ii < static_cast<long>(c); ++ii)
                s += (((ii - nn) % 2) ? -1.0 : 1.0) * nchoosek_r(ii, nn) * Un[ii];
            Pi[nn] = s;
        }
    }
    double PiTot = D / omega / (1.0 - omega);
    for (long nn = 0; nn + 1 < static_cast<long>(c); ++nn) PiTot += Pi[nn];
    bool neg = false;
    for (long nn = 0; nn < std::max<long>(static_cast<long>(c) - 1, 1); ++nn)
        if (Pi[nn] < -1e-10) neg = true;
    if (std::fabs(PiTot - 1.0) > 1e-8 || neg) {
        char buf[256];
        std::snprintf(buf, sizeof buf,
                      "qsys_gimc: the arrival distribution does not normalize (sum = %.12g): the alternating "
                      "sums lost accuracy at c = %u",
                      PiTot, c);
        throw NumericError(ErrorKind{"Cancellation"}, buf);
    }
    for (double& v : Pi) v = std::max(v, 0.0);
 
    std::vector<double> Pn(Pi.size(), 0.0);
    for (long nn = 1; nn <= nmax; ++nn) Pn[nn] = A * Pi[nn - 1] / std::min<long>(nn, c);
    double sj = 0.0;
    for (unsigned j = 1; j + 1 <= c; ++j) sj += Pi[j - 1] * (1.0 / j - 1.0 / c);
    Pn[0] = 1.0 - A / c - A * sj;
 
    const double PW = D / (1.0 - omega);
    const double theta = c * mu * (1.0 - omega);
    const double EW = PW / theta;
    const double EW2 = 2.0 * PW / (theta * theta);
    const double VW = EW2 - EW * EW;
    const double tailc = (A / c) * D;
    const double om1 = 1.0 - omega;
    const double EQ = tailc / (om1 * om1);
    const double EQ2 = tailc * (1.0 + omega) / (om1 * om1 * om1);
    const double m = c + 1.0;
    double EN2 = 0.0;
    for (unsigned n0 = 0; n0 <= c; ++n0) EN2 += static_cast<double>(n0) * n0 * Pn[n0];
    EN2 += tailc * (m * m / om1 + 2.0 * m * omega / (om1 * om1) + omega * (1.0 + omega) / (om1 * om1 * om1));
    const double EN = EQ + A;
 
    QsysRes r;
    r.set("meanQueueLength", EN);
    r.set("varQueueLength", EN2 - EN * EN);
    r.set("meanQueueLengthQ", EQ);
    r.set("varQueueLengthQ", EQ2 - EQ * EQ);
    r.set("meanWaitingTime", EW);
    r.set("varWaitingTime", VW);
    r.set("meanSojournTime", EW + 1.0 / mu);
    r.set("varSojournTime", VW + 1.0 / (mu * mu));
    r.set("probWait", PW);
    r.set("condMeanWaitingTime", 1.0 / theta);
    r.set("throughput", lambda);
    r.set("utilization", A / c);
    r.set("omega", omega);
    r.set("queueLengthDist", Pn);
    r.set("arrivalDist", Pi);
    const std::function<double(double)> Wccdf = [PW, theta](double t) {
        return t < 0.0 ? 1.0 : PW * std::exp(-theta * t);
    };
    const std::function<double(double)> Tccdf = [PW, theta, mu](double t) {
        if (t < 0.0) return 1.0;
        if (std::fabs(theta - mu) > 1e-12 * mu) {
            const double G = PW * mu / (mu - theta);
            return (1.0 - G) * std::exp(-mu * t) + G * std::exp(-theta * t);
        }
        return (1.0 + PW * mu * t) * std::exp(-mu * t);
    };
    if (!o.tPoints.empty()) {
        r.set("tPoints", o.tPoints);
        r.set("waitingTimeCdf", cdf_of(Wccdf, o.tPoints));
        r.set("sojournTimeCdf", cdf_of(Tccdf, o.tPoints));
    }
    if (!o.percentiles.empty()) {
        r.set("percentiles", o.percentiles);
        std::vector<double> wp(o.percentiles.size());
        for (std::size_t i = 0; i < wp.size(); ++i)
            wp[i] = std::max(0.0, std::log(PW / (1.0 - o.percentiles[i] / 100.0)) / theta);
        r.set("waitingTimePercentile", wp);
        r.set("sojournTimePercentile", pct_of(Tccdf, o.percentiles, 50.0 / std::min(mu, theta)));
    }
    r.analyzer = "exact:GI/M/" + std::to_string(c);
    return r;
}
 
/**
 * @brief Exact GI/M/c queue, interarrival law given as a Distribution (its LST,
 *        and its rate 1/mean unless `lambda` > 0 is given, as MATLAB's
 *        LA.evalLST / LA.getMean).
 */
inline QsysRes qsys_gimc(const lang::Distrib<double>& LA, double mu, unsigned c,
                         const QsysResOptions& o = QsysResOptions(), double lambda = -1.0) {
    const lang::Distrib<double> d = LA;
    if (!(lambda > 0.0)) lambda = 1.0 / lang::dist_moment(d, 1);
    return qsys_gimc([d](double s) { return lang::dist_lst(d, s); }, lambda, mu, c, o);
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_GIMC_H