polling_qsys_table.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_POLLING_POLLING_QSYS_TABLE_H
#define LINE_API_POLLING_POLLING_QSYS_TABLE_H
 
/**
 * @file
 * @ingroup api_polling
 * Exact mean waiting times of a polling system with a polling table.
 *
 * Templated port of matlab/src/api/polling/polling_qsys_table.m. The server
 * visits the queues in the order table[0], table[1], ..., table[M-1] and then
 * starts over, every queue appearing at least once and possibly several
 * times; all queues are served exhaustively or all gated.
 *
 * Method (buffer occupancy, Eisenberg 1972): the first and second moments of
 * the queue lengths at the polling instants obey an affine recursion over the
 * table positions. Composing the M steps gives a linear system of dimension
 * n + n^2 whose solution is the cycle-stationary moment vector, then
 * propagated to every position. The waiting time follows by length biasing:
 *
 *   exhaustive  W_i = lambda_i b2_i/(2(1-rho_i)) + sum_m E[X_i(X_i-1)]_m / (2 lambda_i sum_m E[X_i]_m)
 *   gated       W_i = (1+rho_i) sum_m E[X_i(X_i-1)]_m / (2 lambda_i sum_m E[X_i]_m)
 *
 * the sums running over the positions m of queue i. For table = 1..n this is
 * the classical cyclic system (Takagi 1988, eqs. 15 and 21), which the tests
 * use as the oracle together with the CTMC of the same station.
 *
 * Everything is rational in the input moments plus one linear solve, so the
 * port is instantiable at an exact field.
 */
 
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/api/polling/polling_qsys_1limited.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace polling {
 
/** Result of polling_qsys_table: waiting times and polling-instant moments. */
template <class T>
struct PollingTableResult {
    std::vector<T> W;                ///< mean waiting time per queue
    std::vector<std::vector<T> > EX; ///< EX[m][i] = E[X_i] at the polling instant of position m
    std::vector<Matrix<T> > EX2;     ///< EX2[m](i,j) = E[X_i X_j] at the polling instant of position m
};
 
/**
 * Polling table, exhaustive or gated.
 *
 * @param lambda  arrival rate per queue (Poisson)
 * @param b       mean service time per queue
 * @param b2      second raw moment of the service time per queue
 * @param table   queue indices, 1-based as in MATLAB, e.g. {1, 2, 1, 3}
 * @param s_mean  mean switchover of the leg out of each table position
 * @param s_var   variance of the switchover of the leg out of each position
 * @param gated   false for exhaustive, true for gated
 */
template <class T>
PollingTableResult<T> polling_qsys_table_full(const std::vector<T>& lambda, const std::vector<T>& b,
                                              const std::vector<T>& b2,
                                              const std::vector<std::size_t>& table,
                                              const std::vector<T>& s_mean,
                                              const std::vector<T>& s_var, bool gated) {
    static const char* who = "polling_qsys_table";
    const std::size_t n = lambda.size();
    const std::size_t M = table.size();
    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);
    if (n == 0 || b.size() != n || b2.size() != n)
        throw InputError(ErrorKind{"Argument"}, std::string(who) + ": lambda, b and b2 must have one entry per queue");
    if (s_mean.size() != M || s_var.size() != M)
        throw InputError(ErrorKind{"Argument"}, std::string(who) + ": polling_qsys_table needs one switchover per table position (" +
                         std::to_string(M) + ").");
    {
        std::vector<bool> seen(n, false);
        for (std::size_t p : table) {
            if (p < 1 || p > n)
                throw InputError(ErrorKind{"Argument"}, std::string(who) + ": The polling table must list every queue 1.." +
                                 std::to_string(n) + " at least once.");
            seen[p - 1] = true;
        }
        for (std::size_t i = 0; i < n; ++i)
            if (!seen[i])
                throw InputError(ErrorKind{"Argument"}, std::string(who) + ": The polling table must list every queue 1.." +
                                 std::to_string(n) + " at least once.");
    }
    std::vector<T> rho1(n), s2(M);
    T rho = zero;
    for (std::size_t i = 0; i < n; ++i) {
        if (lambda[i] <= zero)
            throw InputError(ErrorKind{"Argument"}, std::string(who) + ": polling_qsys_table needs a positive arrival rate at every queue.");
        rho1[i] = lambda[i] * b[i];
        rho += rho1[i];
    }
    if (rho >= one)
        throw NumericError(ErrorKind{"Instability"}, std::string(who) + ": The polling system is unstable: rho >= 1.");
    bool anysw = false;
    for (std::size_t m = 0; m < M; ++m) {
        s2[m] = s_var[m] + s_mean[m] * s_mean[m];
        if (s_mean[m] != zero) anysw = true;
    }
    if (!anysw)
        throw InputError(ErrorKind{"Argument"}, std::string(who) + ": polling_qsys_table needs a positive total switchover time per cycle.");
 
    // visit period moments per queue: a service time (gated) or a busy period (exhaustive)
    std::vector<T> y1(n), y2(n);
    for (std::size_t q = 0; q < n; ++q) {
        if (gated) {
            y1[q] = b[q];
            y2[q] = b2[q];
        } else {
            const T omr = one - rho1[q];
            y1[q] = b[q] / omr;
            y2[q] = b2[q] / (omr * omr * omr);
        }
    }
 
    const std::size_t N = n + n * n;
    // S is stored column-major after the n means, as MATLAB's reshape does
    const auto Sidx = [n](std::size_t i, std::size_t j) { return n + j * n + i; };
    // visit of queue q at position m followed by the switchover leg m
    const auto step = [&](const std::vector<T>& zi, std::size_t q, std::size_t m) {
        std::vector<T> lp(lambda);
        if (!gated) lp[q] = zero;
        const T muq = zi[q];
        const T Sqq = zi[Sidx(q, q)];
        const T EV = y1[q] * muq;
        const T EV2 = muq * (y2[q] - y1[q] * y1[q]) + Sqq * y1[q] * y1[q];
        std::vector<T> mu1(n);
        for (std::size_t i = 0; i < n; ++i) mu1[i] = (i == q ? zero : zi[i]) + lp[i] * y1[q] * muq;
        std::vector<T> zo(N, zero);
        for (std::size_t j = 0; j < n; ++j) {
            for (std::size_t i = 0; i < n; ++i) {
                const T DSD = (i != q && j != q) ? zi[Sidx(i, j)] : zero;
                const T DSq = (i != q) ? zi[Sidx(i, q)] : zero;
                const T SqD = (j != q) ? zi[Sidx(q, j)] : zero;
                T s1v = DSD + DSq * y1[q] * lp[j] + lp[i] * y1[q] * SqD + lp[i] * lp[j] * EV2;
                if (i == j) s1v += lp[i] * EV;
                T s2v = s1v + mu1[i] * lambda[j] * s_mean[m] + lambda[i] * mu1[j] * s_mean[m] +
                        lambda[i] * lambda[j] * s2[m];
                if (i == j) s2v += lambda[i] * s_mean[m];
                zo[Sidx(i, j)] = s2v;
            }
        }
        for (std::size_t i = 0; i < n; ++i) zo[i] = mu1[i] + lambda[i] * s_mean[m];
        return zo;
    };
 
    std::vector<Matrix<T> > A(M, Matrix<T>(N, N, zero));
    std::vector<std::vector<T> > c(M);
    for (std::size_t m = 0; m < M; ++m) {
        const std::size_t q = table[m] - 1;
        c[m] = step(std::vector<T>(N, zero), q, m);
        for (std::size_t k = 0; k < N; ++k) {
            std::vector<T> e(N, zero);
            e[k] = one;
            const std::vector<T> col = step(e, q, m);
            for (std::size_t r = 0; r < N; ++r) A[m](r, k) = col[r] - c[m][r];
        }
    }
    const auto matvec = [N, zero](const Matrix<T>& X, const std::vector<T>& v) {
        std::vector<T> out(N, zero);
        for (std::size_t r = 0; r < N; ++r)
            for (std::size_t k = 0; k < N; ++k) out[r] += X(r, k) * v[k];
        return out;
    };
    // cycle map from position 1 back to position 1
    Matrix<T> Acyc(N, N, zero);
    for (std::size_t r = 0; r < N; ++r) Acyc(r, r) = one;
    std::vector<T> ccyc(N, zero);
    for (std::size_t m = 0; m < M; ++m) {
        std::vector<T> nc = matvec(A[m], ccyc);
        for (std::size_t r = 0; r < N; ++r) nc[r] += c[m][r];
        ccyc = nc;
        Matrix<T> P(N, N, zero);
        for (std::size_t r = 0; r < N; ++r)
            for (std::size_t k = 0; k < N; ++k) {
                T s = zero;
                for (std::size_t l = 0; l < N; ++l) s += A[m](r, l) * Acyc(l, k);
                P(r, k) = s;
            }
        Acyc = P;
    }
    Matrix<T> IA(N, N, zero);
    for (std::size_t r = 0; r < N; ++r)
        for (std::size_t k = 0; k < N; ++k) IA(r, k) = (r == k ? one : zero) - Acyc(r, k);
    std::vector<std::vector<T> > z(M);
    z[0] = solve(IA, ccyc);
    for (std::size_t m = 0; m + 1 < M; ++m) {
        z[m + 1] = matvec(A[m], z[m]);
        for (std::size_t r = 0; r < N; ++r) z[m + 1][r] += c[m][r];
    }
 
    PollingTableResult<T> out;
    out.EX.assign(M, std::vector<T>(n, zero));
    out.EX2.assign(M, Matrix<T>(n, n, zero));
    for (std::size_t m = 0; m < M; ++m)
        for (std::size_t i = 0; i < n; ++i) {
            out.EX[m][i] = z[m][i];
            for (std::size_t j = 0; j < n; ++j)
                out.EX2[m](i, j) = (z[m][Sidx(i, j)] + z[m][Sidx(j, i)]) / two;
        }
    out.W.assign(n, zero);
    for (std::size_t i = 0; i < n; ++i) {
        T num = zero, den = zero;
        for (std::size_t m = 0; m < M; ++m) {
            if (table[m] != i + 1) continue;
            num += out.EX2[m](i, i) - out.EX[m][i];
            den += out.EX[m][i];
        }
        const T resid = num / (two * lambda[i] * den);
        if (gated)
            out.W[i] = (one + rho1[i]) * resid;
        else
            out.W[i] = lambda[i] * b2[i] / (two * (one - rho1[i])) + resid;
    }
    return out;
}
 
/** Mean waiting times only; see polling_qsys_table_full. */
template <class T>
std::vector<T> polling_qsys_table(const std::vector<T>& lambda, const std::vector<T>& b,
                                  const std::vector<T>& b2, const std::vector<std::size_t>& table,
                                  const std::vector<T>& s_mean, const std::vector<T>& s_var,
                                  bool gated) {
    return polling_qsys_table_full(lambda, b, b2, table, s_mean, s_var, gated).W;
}
 
}  // namespace polling
}  // namespace line
 
#endif  // LINE_API_POLLING_POLLING_QSYS_TABLE_H