qsys_mmapmap1.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_MMAPMAP1_H
#define LINE_API_QSYS_QSYS_MMAPMAP1_H
 
/**
 * @file
 * @ingroup api_qsys
 * The MMAP[K]/MAP[K]/1 FCFS queue in which each class's service phase is CARRIED.
 *
 * THE SERVICE LAW. Each class owns ONE service process whose successive draws are
 * chained in the start order of THAT class: a class-k job starts in the phase the
 * previous class-k draw ended in, whatever was served in between. That is how JMT
 * and both LDES engines sample a MAP service, and the law SolverCTMC follows
 * through the per-class carried-phase local variable. A renewal class keeps its
 * marginal phase-type law, which is exact for it.
 *
 * THE REDUCTION. Under FCFS at one server the service-start order IS the arrival
 * order, so the landing phase of each class-k draw can be decided at the ARRIVAL of
 * its job. With V_k = (-S0)^-1 S1 a draw started in h ends in j with probability
 * V_k(h,j); conditioned on (h,j) it is phase type: S0 restricted to the phases i
 * with V_k(i,j) > 0, moved by the Doob transform S0(i,i')V_k(i',j)/V_k(i,j) and
 * started in h. The queue is a MMAP[K']/PH[K']/1 FCFS queue whose arrival
 * environment is (arrival phase, carried phase of every carried class) and whose
 * marks are (k,h,j); mmapph1fcfs (He, 2012) solves it exactly and the class
 * measures are the sums over each class's marks.
 *
 * References:
 *   Q. He, "Analysis of a continuous time SM[K]/PH[K]/1/FCFS queue: Age process,
 *   sojourn times, and queue lengths", Journal of Systems Science and Complexity,
 *   25(1):133-155, 2012.
 */
 
#include <cmath>
#include <cstddef>
#include <string>
#include <vector>
 
#include "line/api/mam/map_moment.h"
#include "line/api/mam/mmap_lambda.h"
#include "line/api/mam/mmapph1fcfs.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 {
 
/** Return value of qsys_mmapmap1, mirroring the MATLAB struct. */
template <class T>
struct MmapMap1Result {
    std::vector<T> meanQueueLength;  ///< E[N_k]
    std::vector<T> meanSojournTime;  ///< E[W_k], by Little's law
    std::vector<T> arrivalRate;      ///< class arrival rates
    std::vector<bool> carry;         ///< classes whose service phase is carried
    std::size_t markCount;           ///< marks of the lifted MMAP
};
 
namespace mmapmap1detail {
 
/** D1 = (-D0 e) alpha to within tolerance, as mam_is_renewal_map. */
template <class T>
bool is_renewal(const mam::Map<T>& m) {
    const std::size_t n = m.D0.rows();
    if (n == 1) return true;
    const std::vector<T> alpha = mam::map_pie(m);
    double err = 0.0, nrm = 0.0;
    for (std::size_t i = 0; i < n; ++i) {
        double ti = 0.0;
        for (std::size_t j = 0; j < n; ++j) ti += num_traits<T>::to_double(m.D1(i, j));
        for (std::size_t j = 0; j < n; ++j) {
            const double d1 = num_traits<T>::to_double(m.D1(i, j));
            const double e = d1 - ti * num_traits<T>::to_double(alpha[j]);
            err += e * e;
            nrm += d1 * d1;
        }
    }
    return std::sqrt(err) < 1e-9 * std::max(1.0, std::sqrt(nrm));
}
 
}  // namespace mmapmap1detail
 
/**
 * MMAP[K]/MAP[K]/1 FCFS with per-class carried service phases, exactly.
 *
 * @param arrival  marked arrival process (D0, D1, Dc[0..K-1])
 * @param svc      svc[k] = (S0, S1), the service MAP of class k
 * @param carry    per-class flags of the carried classes; empty for every arriving
 *                 class whose MAP is not renewal
 */
template <class T>
MmapMap1Result<T> qsys_mmapmap1(const mam::Mmap<T>& arrival, const std::vector<mam::Map<T> >& svc,
                                std::vector<bool> carry = std::vector<bool>()) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const std::size_t K = arrival.Dc.size();
    if (svc.size() != K)
        throw InputError("qsys_mmapmap1: svc must give one service MAP per arrival mark");
    const std::size_t ma = arrival.D0.rows();
    const std::vector<T> piA = mam::map_prob(mam::Map<T>{arrival.D0, arrival.D1});
    std::vector<T> lamk(K, zero);
    for (std::size_t k = 0; k < K; ++k)
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) lamk[k] += piA[i] * arrival.Dc[k](i, j);
    if (carry.empty()) {
        carry.assign(K, false);
        for (std::size_t k = 0; k < K; ++k) carry[k] = !mmapmap1detail::is_renewal(svc[k]);
    }
    if (carry.size() != K) throw InputError("qsys_mmapmap1: carry must give one flag per class");
    // A class that never arrives draws no service, so it has nothing to carry
    for (std::size_t k = 0; k < K; ++k) carry[k] = carry[k] && lamk[k] > zero;
 
    std::vector<std::size_t> pk(K, 1), cidx;
    std::size_t P = 1;
    for (std::size_t k = 0; k < K; ++k) {
        if (!carry[k]) continue;
        pk[k] = svc[k].D0.rows();
        cidx.push_back(k);
        P *= pk[k];
    }
 
    mam::Mmap<T> lifted;
    lifted.D0 = mam::kron(arrival.D0, eye<T>(P));
    lifted.D1 = Matrix<T>(ma * P, ma * P, zero);
    std::vector<mam::PhService<T> > ph;
    std::vector<std::size_t> markClass;
    const T tiny = num_traits<T>::from_double(1e-14);
    for (std::size_t k = 0; k < K; ++k) {
        const Matrix<T>& S0 = svc[k].D0;
        const Matrix<T>& S1 = svc[k].D1;
        if (!carry[k]) {
            lifted.Dc.push_back(mam::kron(arrival.Dc[k], eye<T>(P)));
            ph.push_back(mam::PhService<T>{mam::map_pie(svc[k]), S0});
            markClass.push_back(k);
            continue;
        }
        const std::size_t p = pk[k];
        Matrix<T> negS0(p, p, zero);
        for (std::size_t i = 0; i < p; ++i)
            for (std::size_t j = 0; j < p; ++j) negS0(i, j) = -S0(i, j);
        Matrix<T> V = matmul(inverse(negS0), S1);
        for (std::size_t i = 0; i < p; ++i)
            for (std::size_t j = 0; j < p; ++j)
                if (V(i, j) < tiny && V(i, j) > -tiny) V(i, j) = zero;
        for (std::size_t h = 0; h < p; ++h) {
            for (std::size_t j = 0; j < p; ++j) {
                if (!(V(h, j) > zero)) continue;
                // carried coordinate of class k moves h -> j, the others stay
                Matrix<T> E = eye<T>(1);
                for (std::size_t kk : cidx) {
                    Matrix<T> Ek;
                    if (kk == k) {
                        Ek = Matrix<T>(p, p, zero);
                        Ek(h, j) = V(h, j);
                    } else {
                        Ek = eye<T>(pk[kk]);
                    }
                    E = mam::kron(E, Ek);
                }
                std::vector<std::size_t> Jset;
                for (std::size_t i = 0; i < p; ++i)
                    if (V(i, j) > zero) Jset.push_back(i);
                const std::size_t nJ = Jset.size();
                Matrix<T> G(nJ, nJ, zero);
                std::vector<T> sig(nJ, zero);
                for (std::size_t a = 0; a < nJ; ++a) {
                    const std::size_t i = Jset[a];
                    G(a, a) = S0(i, i);
                    for (std::size_t bb = 0; bb < nJ; ++bb) {
                        if (bb == a) continue;
                        const std::size_t i2 = Jset[bb];
                        G(a, bb) = S0(i, i2) * V(i2, j) / V(i, j);
                    }
                    if (i == h) sig[a] = one;
                }
                lifted.Dc.push_back(mam::kron(arrival.Dc[k], E));
                ph.push_back(mam::PhService<T>{sig, G});
                markClass.push_back(k);
            }
        }
    }
    for (std::size_t m = 0; m < lifted.Dc.size(); ++m)
        for (std::size_t i = 0; i < ma * P; ++i)
            for (std::size_t j = 0; j < ma * P; ++j) lifted.D1(i, j) += lifted.Dc[m](i, j);
 
    const std::vector<T> Qm = mam::mmapph1fcfs_ncmean(lifted, ph);
    MmapMap1Result<T> r;
    r.meanQueueLength.assign(K, zero);
    r.meanSojournTime.assign(K, zero);
    for (std::size_t m = 0; m < markClass.size(); ++m) r.meanQueueLength[markClass[m]] += Qm[m];
    for (std::size_t k = 0; k < K; ++k)
        if (lamk[k] > zero) r.meanSojournTime[k] = r.meanQueueLength[k] / lamk[k];
    r.arrivalRate = lamk;
    r.carry = carry;
    r.markCount = markClass.size();
    return r;
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMAPMAP1_H