qsys_mg1_sjf.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_SJF_H
#define LINE_API_QSYS_QSYS_MG1_SJF_H
 
/**
 * @file
 * @ingroup api_qsys
 * M/G/1 with non-preemptive shortest-job-first scheduling: per-class mean
 * response times.
 *
 * Port of matlab/src/api/qsys/qsys_mg1_sjf.m. A job of size x waits
 *
 *   W(x) = (lambda E[S^2]/2) / (1 - rho(x))^2,   rho(x) = lambda int_0^x t f(t) dt,
 *
 * (Phipps, 1956), and the class mean averages W(x) over the class's size law.
 * Each class's size law is the reference's two-moment fit (exponential,
 * balanced-means H2 above cs = 1, an Erlang mixture below), and both integrals
 * are the reference's trapezoidal rules on the same uniform grid, so the
 * numbers agree with MATLAB to rounding.
 *
 * ARITHMETIC: double, as the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/api/qsys/qsys_res.h"
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
namespace sjfdetail {
 
/** The reference's sjf_fit: exponential, H2 or Erlang mixture by cs^2. */
struct SjfFit {
    int type = 0;  // 0 exp, 1 h2, 2 erlmix
    double rate = 0, p = 0, r1 = 0, r2 = 0, rate_min = 0, rate_max = 0;
    int k = 0;
};
 
inline SjfFit sjf_fit(double mu, double cs) {
    SjfFit f;
    const double mean_x = 1.0 / mu;
    const double c2 = cs * cs;
    if (std::fabs(c2 - 1.0) < 1e-9) {
        f.type = 0;
        f.rate = f.rate_min = f.rate_max = mu;
    } else if (c2 > 1.0) {
        const double pr = 0.5 * (1.0 + std::sqrt((c2 - 1.0) / (c2 + 1.0)));
        f.type = 1;
        f.p = pr;
        f.r1 = 2.0 * pr * mu;
        f.r2 = 2.0 * (1.0 - pr) * mu;
        f.rate_min = std::min(f.r1, f.r2);
        f.rate_max = std::max(f.r1, f.r2);
    } else {
        const int k = static_cast<int>(std::ceil(1.0 / c2));
        const double pr = (1.0 / (1.0 + c2)) * (k * c2 - std::sqrt(k * (1.0 + c2) - static_cast<double>(k) * k * c2));
        f.type = 2;
        f.k = k;
        f.p = pr;
        f.rate = f.rate_min = f.rate_max = (k - pr) / mean_x;
    }
    return f;
}
 
/** Erlang(n, rate) density at x, in the reference's log-Poisson form. */
inline double erlang_pdf(int n, double rate, double x) {
    if (n <= 0) return 0.0;
    const double t = rate * x;
    const int m = n - 1;
    double logp;
    if (t <= 0.0) logp = (m == 0) ? 0.0 : -std::numeric_limits<double>::infinity();
    else logp = m * std::log(t) - t - std::lgamma(m + 1.0);
    return rate * std::exp(logp);
}
 
inline double sjf_pdf(const SjfFit& f, double x) {
    switch (f.type) {
        case 0: return f.rate * std::exp(-f.rate * x);
        case 1: return f.p * f.r1 * std::exp(-f.r1 * x) + (1.0 - f.p) * f.r2 * std::exp(-f.r2 * x);
        default: return f.p * erlang_pdf(f.k - 1, f.rate, x) + (1.0 - f.p) * erlang_pdf(f.k, f.rate, x);
    }
}
 
}  // namespace sjfdetail
 
/** MATLAB's [W, rho_total, res] outputs. */
struct Mg1SjfResult {
    std::vector<double> W; ///< per-class mean response time
    double rhohat = 0;     ///< Q/(1+Q), MATLAB's rho_total
    QsysRes res;           ///< meanSojournTime, meanWaitingTime, meanQueueLength, utilization, conservationResidual, analyzer
};
 
/**
 * @brief M/G/1 non-preemptive SJF.
 *
 * @param lambda per-class arrival rates
 * @param mu     per-class service rates
 * @param cs     per-class coefficients of variation of the service time
 */
inline Mg1SjfResult qsys_mg1_sjf(const std::vector<double>& lambda, const std::vector<double>& mu,
                                 const std::vector<double>& cs) {
    using namespace sjfdetail;
    const std::size_t K = lambda.size();
    if (mu.size() != K || cs.size() != K)
        throw InputError(ErrorKind("Argument"), "qsys_mg1_sjf: lambda, mu and cs must have the same length");
    for (std::size_t r = 0; r < K; ++r)
        if (!(lambda[r] > 0.0) || !(mu[r] > 0.0) || !(cs[r] >= 0.0))
            throw InputError(ErrorKind("Argument"), "qsys_mg1_sjf: lambda and mu must be positive, cs non-negative");
    double lambda_total = 0.0, rho_util = 0.0, lES2 = 0.0;
    for (std::size_t r = 0; r < K; ++r) {
        lambda_total += lambda[r];
        rho_util += lambda[r] / mu[r];
        lES2 += lambda[r] * (1.0 + cs[r] * cs[r]) / (mu[r] * mu[r]);
    }
    if (rho_util >= 1.0) resdetail::throw_unstable("qsys_mg1_sjf", "rho", rho_util);
    std::vector<SjfFit> fits(K);
    double rate_min = std::numeric_limits<double>::infinity(), rate_max = 0.0;
    for (std::size_t r = 0; r < K; ++r) {
        fits[r] = sjf_fit(mu[r], cs[r]);
        rate_min = std::min(rate_min, fits[r].rate_min);
        rate_max = std::max(rate_max, fits[r].rate_max);
    }
    const double xmax = 40.0 / rate_min;
    const long N = std::min(2000000L, std::max(20000L, static_cast<long>(std::ceil(200.0 * rate_max / rate_min))));
    const std::size_t np = static_cast<std::size_t>(N) + 1;
    std::vector<double> x(np), fmix(np, 0.0);
    // linspace(0, xmax, N+1), as MATLAB computes it
    for (std::size_t i = 0; i < np; ++i) x[i] = xmax * static_cast<double>(i) / static_cast<double>(N);
    x[np - 1] = xmax;
    std::vector<std::vector<double>> pdf(K, std::vector<double>(np));
    for (std::size_t r = 0; r < K; ++r) {
        const double p = lambda[r] / lambda_total;
        for (std::size_t i = 0; i < np; ++i) {
            pdf[r][i] = sjf_pdf(fits[r], x[i]);
            fmix[i] += p * pdf[r][i];
        }
    }
    // rho_x = lambda_total * cumtrapz(x, x .* fmix)
    std::vector<double> Wx(np);
    double cum = 0.0;
    for (std::size_t i = 0; i < np; ++i) {
        if (i > 0) cum += 0.5 * (x[i] - x[i - 1]) * (x[i] * fmix[i] + x[i - 1] * fmix[i - 1]);
        const double denom = std::max(1.0 - lambda_total * cum, 1e-12);
        Wx[i] = (lES2 / 2.0) / (denom * denom);
    }
    auto trapz = [&](const std::vector<double>& y) {
        double s = 0.0;
        for (std::size_t i = 1; i < np; ++i) s += 0.5 * (x[i] - x[i - 1]) * (y[i] + y[i - 1]);
        return s;
    };
    Mg1SjfResult out;
    std::vector<double> Wq(K), Qi(K), y(np);
    out.W.assign(K, 0.0);
    double Q = 0.0;
    for (std::size_t r = 0; r < K; ++r) {
        for (std::size_t i = 0; i < np; ++i) y[i] = Wx[i] * pdf[r][i];
        Wq[r] = trapz(y);
        out.W[r] = Wq[r] + 1.0 / mu[r];
        Qi[r] = lambda[r] * out.W[r];
        Q += Qi[r];
    }
    out.rhohat = Q / (1.0 + Q);
    for (std::size_t i = 0; i < np; ++i) y[i] = x[i] * fmix[i] * Wx[i];
    out.res.set("meanSojournTime", out.W);
    out.res.set("meanWaitingTime", Wq);
    out.res.set("meanQueueLength", Qi);
    out.res.set("utilization", rho_util);
    out.res.set("conservationResidual", lambda_total * trapz(y) - rho_util * lES2 / (2.0 * (1.0 - rho_util)));
    out.res.analyzer = "exact:M/G/1-SJF";
    return out;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MG1_SJF_H