qsys_mapmapc.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_MAPMAPC_H
#define LINE_API_QSYS_QSYS_MAPMAPC_H
 
/**
 * @file
 * @ingroup api_qsys
 * The MAP/MAP/c FCFS queue whose service MAP is CARRIED in service-start order.
 *
 * THE SERVICE LAW. The service MAP is one station-level process whose successive
 * draws are chained in service-START order, idle periods included: each job starts
 * in the phase the previous draw ENDED in, which is how JMT and both LDES engines
 * sample a MAP service and the law SolverCTMC follows through its pair form. A
 * renewal reading of the marginal (qsys_mapphc, or the single-fast-server
 * surrogate) discards the autocorrelation of the service times.
 *
 * THE STATE SPACE. With V = (-S0)^-1 S1, a draw started in h lands in j with
 * probability V(h,j); a busy server is a PAIR (i,j) of current phase i and landing
 * j, run under S0 conditioned on ending in j. A level holds the MULTISET of busy
 * pairs and the landing phase of the draw started last. `ldqbd_mphc(..., carry)`
 * builds those blocks; this header lifts them to a MAP arrival stream and solves
 * the level-independent QBD above level c matrix-geometrically, exactly as
 * qsys_mapphc does for a renewal service. At c = 1 it is the frozen-phase
 * MAP/MAP/1 queue.
 *
 * References:
 *   S. Asmussen and J.R. Moller, "Calculation of the steady state waiting time
 *   distribution in GI/PH/c and MAP/PH/c queues", Queueing Systems 37(1):9-29,
 *   2001.
 *   M.F. Neuts, "Matrix-geometric solutions in stochastic models", Johns Hopkins
 *   University Press, 1981.
 */
 
#include <cstddef>
#include <vector>
 
#include "line/api/mam/ldqbd_mphc.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/qbd_r.h"
#include "line/api/qsys/qsys_mapphc.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace qsys {
 
/**
 * MAP/MAP/c FCFS with the service phase carried in start order, exactly.
 *
 * @param arrival     arrival MAP (D0, D1) of order ma
 * @param service     service MAP (S0, S1) of order p
 * @param c           number of servers, c >= 1
 * @param dist_size   cap on the queue length probabilities materialized
 * @param num_w_moms  how many waiting-time moments to return (0: E[Wq] by Little)
 *
 * Returns the MapPhcResult struct; phaseCount is the repeating configuration count.
 */
template <class T>
MapPhcResult<T> qsys_mapmapc(const mam::Map<T>& arrival, const mam::Map<T>& service, unsigned c,
                             std::size_t dist_size = 500, std::size_t num_w_moms = 3) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mapmapc requires transcendental arithmetic");
    using mapphcdetail::madd2;
    using mapphcdetail::mkron;
 
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const Matrix<T>& D0 = arrival.D0;
    const Matrix<T>& D1 = arrival.D1;
    const std::size_t ma = D0.rows();
    const std::size_t p = service.D0.rows();
    if (D0.cols() != ma || D1.rows() != ma || D1.cols() != ma)
        throw InputError("qsys_mapmapc: D0 and D1 must be square and of equal order");
    if (service.D0.cols() != p || service.D1.rows() != p || service.D1.cols() != p)
        throw InputError("qsys_mapmapc: S0 and S1 must be square and of equal order");
    if (c < 1) throw InputError("qsys_mapmapc: c must be a positive integer");
    if (dist_size == 0) throw InputError("qsys_mapmapc: dist_size must be positive");
 
    const T lambda = mam::map_lambda(arrival);
    if (lambda <= zero) throw InputError("qsys_mapmapc: non-positive arrival rate");
    const T mean_service = one / mam::map_lambda(service);
    const T ct = num_traits<T>::from_int(static_cast<long>(c));
    const T rho = lambda * mean_service / ct;
    if (rho >= one) throw InputError("qsys_mapmapc: load rho must be strictly less than 1");
 
    // Carried-phase blocks at unit arrival rate; Q2[n] leaves level n, so DN_k = Q2[k]
    // for k <= c and Q2[c+1] is the repeating completion that hands the server on.
    const std::vector<T> alpha(p, one / num_traits<T>::from_int(static_cast<long>(p)));
    const std::vector<T> arr(c + 2, one);
    const mam::LdqbdMphcBlocks<T> B =
        mam::ldqbd_mphc(service.D0, service.D1, alpha, static_cast<double>(c), arr,
                        std::vector<T>(), true);
    std::vector<Matrix<T> > LOC(c + 1);
    for (unsigned k = 0; k <= c; ++k) {
        LOC[k] = B.Q1[k];
        for (std::size_t i = 0; i < LOC[k].rows(); ++i) LOC[k](i, i) += one;
    }
    const Matrix<T>& CDEP = B.Q2[c + 1];
    const std::size_t nc = CDEP.rows();
 
    const Matrix<T> Ima = eye<T>(ma);
    const Matrix<T> Inc = eye<T>(nc);
    const Matrix<T> A_up = mkron(D1, Inc);
    const Matrix<T> A_loc = madd2(mkron(D0, Inc), mkron(Ima, LOC[c]));
    const Matrix<T> A_dn = mkron(Ima, CDEP);
    const Matrix<T> R = mam::qbd_R_logred(A_dn, A_loc, A_up);
 
    const std::size_t n_op = nc * ma;
    Matrix<T> ImR(n_op, n_op, zero);
    for (std::size_t i = 0; i < n_op; ++i)
        for (std::size_t j = 0; j < n_op; ++j) ImR(i, j) = (i == j ? one : zero) - R(i, j);
    const Matrix<T> ImRinv = inverse(ImR);
    const std::vector<T> sum_geom = mulvec(ImRinv, ones<T>(n_op));
 
    std::vector<std::size_t> sz(c + 1), off(c + 2, 0);
    for (unsigned k = 0; k <= c; ++k) {
        sz[k] = ma * LOC[k].rows();
        off[k + 1] = off[k] + sz[k];
    }
    const std::size_t tot = off[c + 1];
    Matrix<T> Q(tot, tot, zero);
    const Matrix<T> RA_dn = matmul(R, A_dn);
    for (unsigned k = 0; k <= c; ++k) {
        const Matrix<T> Ick = eye<T>(LOC[k].rows());
        const Matrix<T> diag = (k < c) ? madd2(mkron(D0, Ick), mkron(Ima, LOC[k]))
                                       : madd2(A_loc, RA_dn);
        for (std::size_t i = 0; i < sz[k]; ++i)
            for (std::size_t j = 0; j < sz[k]; ++j) Q(off[k] + i, off[k] + j) += diag(i, j);
        if (k < c) {
            const Matrix<T> up = mkron(D1, B.Q0[k]);
            for (std::size_t i = 0; i < up.rows(); ++i)
                for (std::size_t j = 0; j < up.cols(); ++j) Q(off[k] + i, off[k + 1] + j) += up(i, j);
        }
        if (k > 0) {
            const Matrix<T> dn = mkron(Ima, B.Q2[k]);
            for (std::size_t i = 0; i < dn.rows(); ++i)
                for (std::size_t j = 0; j < dn.cols(); ++j) Q(off[k] + i, off[k - 1] + j) += dn(i, j);
        }
    }
 
    // pi Q = 0 with the normalization over boundary and tail replacing the last equation
    Matrix<T> M(tot, tot, zero);
    for (std::size_t i = 0; i < tot; ++i)
        for (std::size_t j = 0; j < tot; ++j) M(j, i) = Q(i, j);
    for (std::size_t col = 0; col < tot; ++col) M(tot - 1, col) = zero;
    for (std::size_t i = 0; i < off[c]; ++i) M(tot - 1, i) = one;
    for (std::size_t i = 0; i < sz[c]; ++i) M(tot - 1, off[c] + i) = sum_geom[i];
    std::vector<T> b(tot, zero);
    b[tot - 1] = one;
    const std::vector<T> pi_vec = line::solve(M, b);
 
    std::vector<T> pi_c(n_op);
    for (std::size_t i = 0; i < n_op; ++i) pi_c[i] = pi_vec[off[c] + i];
 
    std::vector<T> ql;
    T acc = zero;
    for (unsigned k = 0; k < c; ++k) {
        T s = zero;
        for (std::size_t i = 0; i < sz[k]; ++i) s += pi_vec[off[k] + i];
        ql.push_back(s);
        acc += s;
    }
    std::vector<T> tail = pi_c;
    {
        T s = zero;
        for (std::size_t i = 0; i < n_op; ++i) s += tail[i];
        ql.push_back(s);
        acc += s;
    }
    while (acc < one - num_traits<T>::from_double(1e-12) && ql.size() < dist_size) {
        std::vector<T> next(n_op, zero);
        for (std::size_t j = 0; j < n_op; ++j)
            for (std::size_t i = 0; i < n_op; ++i) next[j] += tail[i] * R(i, j);
        tail = next;
        T s = zero;
        for (std::size_t i = 0; i < n_op; ++i) s += tail[i];
        ql.push_back(s);
        acc += s;
    }
    // E[N] in closed form over the matrix-geometric tail, as in qsys_mapphc
    T meanQL = zero;
    for (unsigned k = 0; k < c; ++k)
        meanQL += num_traits<T>::from_int(static_cast<long>(k)) * ql[k];
    const std::vector<T> r_tail = mulvec(R, mulvec(ImRinv, sum_geom));
    for (std::size_t i = 0; i < n_op; ++i) meanQL += pi_c[i] * (ct * sum_geom[i] + r_tail[i]);
 
    // Waiting time: an arrival finding j waiting waits for j+1 completions of the
    // full-bank configuration process (LOC_c, CDEP), started at its arrival epoch.
    const Matrix<T> Lj = mkron(Ima, LOC[c]);
    const Matrix<T> Cj = mkron(Ima, CDEP);
    Matrix<T> G0 = matmul(ImRinv, mkron(D1, Inc));
    for (std::size_t i = 0; i < n_op; ++i)
        for (std::size_t j = 0; j < n_op; ++j) G0(i, j) = G0(i, j) / lambda;
    T probWait = zero;
    {
        const std::vector<T> v = mulvec(G0, ones<T>(n_op));
        for (std::size_t i = 0; i < n_op; ++i) probWait += pi_c[i] * v[i];
    }
    std::vector<T> wMoms;
    if (num_w_moms > 0) {
        // X(-Lj) - R X Cj = rhs, vectorized column-major
        const std::size_t n2 = n_op * n_op;
        Matrix<T> Kop(n2, n2, zero);
        for (std::size_t a = 0; a < n_op; ++a) {
            for (std::size_t bcol = 0; bcol < n_op; ++bcol) {
                for (std::size_t i = 0; i < n_op; ++i)
                    Kop(a * n_op + i, bcol * n_op + i) += -Lj(bcol, a);
                for (std::size_t i = 0; i < n_op; ++i)
                    for (std::size_t q = 0; q < n_op; ++q)
                        Kop(a * n_op + i, bcol * n_op + q) += -R(i, q) * Cj(bcol, a);
            }
        }
        Matrix<T> gPrev(n_op, n_op, zero);
        for (std::size_t k = 1; k <= num_w_moms; ++k) {
            std::vector<T> rhs(n2, zero);
            for (std::size_t j = 0; j < n_op; ++j)
                for (std::size_t i = 0; i < n_op; ++i)
                    rhs[j * n_op + i] = (k == 1)
                            ? G0(i, j)
                            : -num_traits<T>::from_int(static_cast<long>(k - 1)) * gPrev(i, j);
            const std::vector<T> gv = line::solve(Kop, rhs);
            Matrix<T> g(n_op, n_op, zero);
            for (std::size_t j = 0; j < n_op; ++j)
                for (std::size_t i = 0; i < n_op; ++i) g(i, j) = gv[j * n_op + i];
            const std::vector<T> v = mulvec(g, ones<T>(n_op));
            T mk = zero;
            for (std::size_t i = 0; i < n_op; ++i) mk += pi_c[i] * v[i];
            mk = mk * num_traits<T>::from_int(static_cast<long>(k));
            if (k % 2 == 0) mk = -mk;
            wMoms.push_back(mk);
            gPrev = g;
        }
    }
 
    MapPhcResult<T> r;
    r.meanQueueLength = meanQL;
    r.meanWaitingTime = wMoms.empty() ? meanQL / lambda - mean_service : wMoms[0];
    r.meanSojournTime = r.meanWaitingTime + mean_service;
    r.utilization = rho;
    r.queueLengthDist = ql;
    r.waitingTimeMoments = wMoms;
    r.probWait = probWait;
    r.phaseCount = nc;
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MAPMAPC_H