qsys_mapmap1.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_MAPMAP1_H
#define LINE_API_QSYS_QSYS_MAPMAP1_H
 
/**
 * @file
 * @ingroup api_qsys
 * The MAP/MAP/1 FCFS queue: mean number in system, waiting time, sojourn time,
 * utilization and the queue-length distribution.
 *
 * ALGORITHM, AND HOW IT DIFFERS FROM THE MATLAB REFERENCE.
 * matlab/src/api/qsys/qsys_mapmap1.m obtains these quantities by calling
 * BUTools' MMAPPH1FCFS, which is not transcribed here; this port computes the
 * SAME quantities from the port's own
 * quasi-birth-death machinery (line/api/mam/qbd_mapmap1.h), which solves the
 * level-independent QBD whose level is the number in system and whose phase is
 * the pair (arrival phase, service phase). The algorithm is therefore
 * different, the quantity is the same.
 *
 * THE SERVICE MAP IS COLLAPSED TO A PHASE-TYPE LAW, as the reference does.
 * Before calling BUTools, qsys_mapmap1.m replaces the service MAP (D0,D1) by
 *     sigma = theta D1 / (theta D1 e),   S = D0,   theta (D0+D1) = 0,
 * i.e. by the stationary marginal of a service time (its mapToPh). That
 * embedded PH keeps the service-time MARGINAL exactly but DISCARDS the serial
 * correlation between consecutive service times, so the reference solves a
 * MAP/PH/1 queue in place of the MAP/MAP/1 queue it names. This port
 * follows the reference: the service MAP is collapsed the same way and the
 * resulting MAP/PH/1 queue is solved as a QBD, so the two agree also on a
 * correlated service MAP (case 3 below). A caller who wants the genuine
 * MAP/MAP/1 answer, with the service correlation honoured, has it in
 * line::mam::qbd_mapmap1 (line/api/mam/qbd_mapmap1.h).
 *
 * The higher moments the reference obtains from MMAPPH1FCFS are returned too:
 * queueLengthMoments (numQLMoms raw moments of the number in system, from
 * pi_0 and R in closed form, qbd_qlen_moment) and sojournTimeMoments
 * (numSTMoms raw moments, from the phase-type sojourn-time law of
 * mmapph1fcfs_stdistr_ph, the same age-process construction BUTools uses).
 *
 * MEASURED AGREEMENT (MATLAB R2025a, T = double). Metric order is
 * meanQueueLength / meanWaitingTime / meanSojournTime / utilization. Arrival
 * MAPs used below:
 *   P(2) = Poisson, D0 = [-2], D1 = [2]
 *   C    = correlated MMPP2, D0 = [-2.5 0.2; 0.1 -0.7], D1 = diag(2.3, 0.6),
 *          lambda = 7/6
 *
 * 1. M/M/1 collapse, qsys_mapmap1(P(2), service D0 = [-3], D1 = [3]). MATLAB
 *    1.999999999999999 / 0.6666666666666656 / 0.9999999999999989 /
 *    0.6666666666666666; port 2 / 0.666666666666668 / 1 / 0.666666666666667,
 *    the textbook values. Relative differences 5e-16, 1.8e-15, 1.1e-15.
 *
 * 2. Renewal service, where the reference's PH reduction is exact. Arrival C,
 *    service the Erlang-2 renewal MAP D0 = [-6 6; 0 -6], D1 = [0 0; 6 0].
 *    MATLAB 0.8364637529649164 / 0.3836355977794533 / 0.7169689311127866 /
 *    0.3888888888888888; port 0.836463752964917 / 0.383635597779453 /
 *    0.716968931112786 / 0.388888888888889. Relative differences below 1e-15
 *    on every metric: with a renewal service the two model the same system and
 *    the port reproduces the reference to machine precision.
 *
 * 3. Correlated service, where the reference's PH reduction is NOT exact.
 *    Arrival C, service MAP D0 = [-5 0.4; 0.2 -1.4], D1 = diag(4.6, 1.2)
 *    (lambda_s = 7/3, so rho = 1/2). MATLAB qsys_mapmap1 reports
 *    2.0305355320573 / 1.311887598906254 / 1.740459027477683 / 0.5. With the
 *    service collapsed to its PH marginal as the reference does, the port
 *    returns the same MAP/PH/1 answer. The genuine MAP/MAP/1 value, which
 *    LINE's own MATLAB qbd_mapmap1 gives as QN = 2.547554867631855, RN =
 *    2.183618457970162, is 20% away: that gap is the service correlation the
 *    reference's PH reduction throws away, and it is available through
 *    line::mam::qbd_mapmap1.
 * meanWaitingTime is meanSojournTime minus the mean service time 1/lambda_s,
 * and utilization is lambda_a/lambda_s, both matching the reference's own
 * definitions.
 *
 * ARITHMETIC. The function is gated on num_traits<T>::has_transcendental
 * because qbd_mapmap1 reaches R by cyclic reduction, a fixed-point iteration
 * driven to a tolerance that does not terminate in a finite number of field
 * operations; running it at exact rational arithmetic would produce a rational
 * with a denominator doubling every step and still not the exact R. See
 * qbd_r.h. Everything downstream of R here -- the mean, the level
 * probabilities, the Little's-law conversions -- is finite exact matrix
 * algebra and carries no additional error.
 */
 
#include <cstddef>
#include <vector>
 
#include "line/api/mam/map_moment.h"
#include "line/api/mam/mmap_lambda.h"
#include "line/api/mam/mmapph1fcfs.h"
#include "line/api/mam/qbd_mapmap1.h"
#include "line/api/mam/qbd_r.h"
#include "line/api/mc/ctmc_solve.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace qsys {
 
/**
 * Return value of the MAP/MAP/1 family (qsys_mapmap1, qsys_mapph1, qsys_phph1),
 * carrying the same quantities as the MATLAB result struct.
 */
template <class T>
struct MapMap1Result {
    T meanQueueLength;                ///< E[N], number in system
    T meanWaitingTime;                ///< E[Wq], time in queue
    T meanSojournTime;                ///< E[W] = E[Wq] + E[S]
    T utilization;                    ///< rho = lambda_a / lambda_s
    std::vector<T> queueLengthDist;   ///< P(N = n), n = 0, 1, ...
    std::vector<T> queueLengthMoments;   ///< raw moments E[N^k], k = 1..numQLMoms
    std::vector<T> sojournTimeMoments;   ///< raw moments E[W^k], k = 1..numSTMoms
};
 
namespace detail {
 
/**
 * The reference's mapToPh: the stationary marginal of a service time of the
 * MAP (D0, D1) as the PH (sigma, S) with sigma = theta D1 / (theta D1 e) and
 * S = D0; uniform sigma when the MAP never fires.
 */
template <class T>
void mapmap1_map_to_ph(const mam::Map<T>& service, std::vector<T>& sigma, Matrix<T>& S) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const std::size_t m = service.D0.rows();
    Matrix<T> Q(m, m, zero);
    for (std::size_t i = 0; i < m; ++i)
        for (std::size_t j = 0; j < m; ++j) Q(i, j) = service.D0(i, j) + service.D1(i, j);
    const std::vector<T> theta = mc::ctmc_solve(Q);
    sigma = vecmul(theta, service.D1);
    T total = zero;
    for (std::size_t j = 0; j < m; ++j) total += sigma[j];
    if (total > zero) {
        for (std::size_t j = 0; j < m; ++j) sigma[j] = sigma[j] / total;
    } else {
        sigma.assign(m, T(one / num_traits<T>::from_int(static_cast<long>(m))));
    }
    S = service.D0;
}
 
/** Raw moments k! alpha (-A)^-k e, k = 1..n, of the PH (alpha, A). */
template <class T>
std::vector<T> mapmap1_ph_moments(const std::vector<T>& alpha, const Matrix<T>& A, std::size_t n) {
    std::vector<T> out;
    if (n == 0) return out;
    const std::size_t m = A.rows();
    Matrix<T> negA(m, m);
    for (std::size_t i = 0; i < m; ++i)
        for (std::size_t j = 0; j < m; ++j) negA(i, j) = -A(i, j);
    const Matrix<T> iA = inverse(negA);
    std::vector<T> v = alpha;
    for (std::size_t k = 1; k <= n; ++k) {
        v = vecmul(v, iA);
        T s = num_traits<T>::from_int(0);
        for (std::size_t j = 0; j < m; ++j) s += v[j];
        out.push_back(T(num_factorial<T>(static_cast<unsigned>(k)) * s));
    }
    return out;
}
 
}  // namespace detail
 
/**
 * MAP/MAP/1 as the reference solves it: the service MAP collapsed to its PH
 * marginal (mapToPh), then the MAP/PH/1 queue by the exact QBD solution.
 * @param arrival    arrival MAP (C0, C1)
 * @param service    service MAP (D0, D1); its correlation is discarded, as in
 *                   the reference (see the file comment)
 * @param dist_size  how many entries of queueLengthDist to materialize
 *                   (the reference's numQLProbs)
 * @param numQLMoms  how many raw queue-length moments to return
 * @param numSTMoms  how many raw sojourn-time moments to return
 */
template <class T>
MapMap1Result<T> qsys_mapmap1(const mam::Map<T>& arrival, const mam::Map<T>& service,
                              std::size_t dist_size, std::size_t numQLMoms,
                              std::size_t numSTMoms) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mapmap1 requires transcendental arithmetic");
    if (dist_size == 0) throw InputError("qsys_mapmap1: dist_size must be positive");
    const T zero = num_traits<T>::from_int(0);
 
    // mapToPh, then the PH as its renewal MAP {S, (-S e) sigma}
    std::vector<T> sigma;
    Matrix<T> S;
    detail::mapmap1_map_to_ph(service, sigma, S);
    const std::size_t m = S.rows();
    mam::Map<T> service_ph;
    service_ph.D0 = S;
    service_ph.D1 = Matrix<T>(m, m, zero);
    for (std::size_t i = 0; i < m; ++i) {
        T exit = zero;
        for (std::size_t j = 0; j < m; ++j) exit -= S(i, j);
        for (std::size_t j = 0; j < m; ++j) service_ph.D1(i, j) = exit * sigma[j];
    }
 
    const mam::QbdMapMap1Result<T> q = mam::qbd_mapmap1(arrival, service_ph, T(zero), dist_size);
 
    // utilization from the arrival and service rates, rho = lambda / mu
    const T lambda_a = mam::map_lambda(arrival);
    const T mu = mam::map_lambda(service);
 
    // mean service time from the PH, sigma (-S)^-1 e
    const std::vector<T> svcMoms = detail::mapmap1_ph_moments(sigma, S, 1);
    const T mean_service = svcMoms[0];
 
    MapMap1Result<T> r;
    r.queueLengthMoments.clear();
    for (std::size_t k = 1; k <= numQLMoms; ++k)
        r.queueLengthMoments.push_back(mam::qbd_qlen_moment(q.pi0, q.R, static_cast<unsigned>(k)));
    r.sojournTimeMoments.clear();
    if (numSTMoms > 0) {
        mam::Mmap<T> mm;
        mm.D0 = arrival.D0;
        mm.D1 = arrival.D1;
        mm.Dc.assign(1, arrival.D1);
        std::vector<mam::PhService<T>> svc(1);
        svc[0].sigma = sigma;
        svc[0].S = S;
        const std::vector<mam::StDistrPh<T>> st = mam::mmapph1fcfs_stdistr_ph(mm, svc);
        r.sojournTimeMoments = detail::mapmap1_ph_moments(st[0].alpha, st[0].A, numSTMoms);
    }
 
    // meanQL = ncMoms(1), meanST = stMoms(1), each 0 when not requested
    r.meanQueueLength = r.queueLengthMoments.empty() ? zero : r.queueLengthMoments[0];
    r.meanSojournTime = r.sojournTimeMoments.empty() ? zero : r.sojournTimeMoments[0];
    const T wait = T(r.meanSojournTime - mean_service);
    r.meanWaitingTime = wait < zero ? zero : wait;
    r.utilization = lambda_a / mu;
    r.queueLengthDist.assign(q.pqueue.rows(), zero);
    for (std::size_t k = 0; k < q.pqueue.rows(); ++k) {
        T s = zero;
        for (std::size_t j = 0; j < q.pqueue.cols(); ++j) s += q.pqueue(k, j);
        r.queueLengthDist[k] = s;
    }
    return r;
}
 
/** qsys_mapmap1 with the reference's defaults numQLMoms = 3, numSTMoms = 3. */
template <class T>
MapMap1Result<T> qsys_mapmap1(const mam::Map<T>& arrival, const mam::Map<T>& service,
                              std::size_t dist_size) {
    return qsys_mapmap1(arrival, service, dist_size, static_cast<std::size_t>(3),
                        static_cast<std::size_t>(3));
}
 
/** qsys_mapmap1 with 100 materialized levels, the reference's numQLProbs. */
template <class T>
MapMap1Result<T> qsys_mapmap1(const mam::Map<T>& arrival, const mam::Map<T>& service) {
    return qsys_mapmap1(arrival, service, static_cast<std::size_t>(100));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MAPMAP1_H