qsys_mmapgk1.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_MMAPGK1_H
#define LINE_API_QSYS_QSYS_MMAPGK1_H
 
/**
 * @file
 * @ingroup api_qsys
 * The MMAP[K]/G[K]/1 FCFS queue: K customer types with class-dependent GENERAL
 * service, fed by a marked Markovian arrival process.
 *
 * WHY THIS IS NOT MMAPPH1FCFS. That routine needs every type's service to be
 * PHASE TYPE, because it builds a QBD whose phase carries the service phase.
 * Here the service laws are arbitrary and may differ in family across types --
 * deterministic for one, uniform for another -- so no finite phase carries
 * them, and the analysis has to run through transforms instead.
 *
 * THE METHOD, which is He's, theorem for theorem. FCFS makes the actual waiting
 * time of a customer the WORKLOAD it finds on arrival, so everything follows
 * from the joint transform of workload and arrival phase,
 * f(s)_j = E[exp(-s V) 1{phase = j}], which by He's Theorem 4.1 (eq. 4.6)
 * satisfies
 *
 *     f(s) [ s I + D0 + sum_k Dk gk(s) ] = s v0,                          (*)
 *
 * with v0 the idle-phase vector, his y0. v0 needs NO root search: the matrix U
 * solving U = D0 + sum_k Dk Fk(U), Fk(U) = int exp(U t) dFk(t), is his
 * eq. (4.4), the generator of the underlying Markov process obtained by
 * EXCISING the busy periods, and eq. (4.5) with Theorem 4.2 give y0 Q = 0 and
 * y0 e = 1 - rho, i.e. v0 = (1 - rho) pi_U. The same vector is what the
 * analyticity of (*) forces, since for every left eigenpair (w, u) of U one has
 * w [D0 + sum_k Dk gk(-u) - u I] = 0, so the roots of the bracket in the closed
 * right half plane are exactly s = -u over the spectrum of U; the two agree to
 * 2.5e-13, and the stationary route is taken because it needs no complex
 * eigenvector.
 *
 * The per-type actual waiting time is the workload seen by a type-k arrival,
 * biased by that type's own arrival block, his Theorem 5.1 eq. (5.1) summed
 * over the post-arrival phase:
 *
 *     E[exp(-s Wk)] = f(s) Dk e / lambda_k.
 *
 * SCOPE. He allows an arrival to be a BATCH carrying a sequence of types, and
 * his Theorem 5.3 then multiplies the transform by prod_{i<n} f*_{h_i}(s), the
 * service of the customers ahead of the tagged one WITHIN its own batch. This
 * header covers the single-customer-per-arrival case, his Special case 3.3,
 * where that product is empty, which is exactly the MMAP convention LINE
 * carries.
 *
 * MOMENTS WITHOUT INVERSION. Differentiating (*) at s = 0 gives
 * sum_i C(j,i) f_i M_{j-i} = [j = 1] v0. M_0 = D is SINGULAR with right null
 * vector e, so each order fixes f_j only up to a multiple of theta, and that
 * multiple is what the NEXT order's solvability condition supplies. At j = 0
 * the same condition reads theta M_1 e = v0 e, i.e. 1 - rho = 1 - rho, which is
 * the identity that validates the whole setup.
 *
 * ARITHMETIC. Gated on num_traits<T>::has_transcendental: the U iteration runs
 * to a tolerance, and both the deterministic transform and the inversion need
 * exp. The CDF is the Abate-Whitt Euler sum, which evaluates the transform OFF
 * the real axis, so a small complex layer is carried here rather than in the
 * distribution interface.
 *
 * Reference:
 *   Qi-Ming He, "The versatility of MMAP[K] and the MMAP[K]/G[K]/1 queue",
 *   Queueing Systems 38(4):397-418, 2001.
 */
 
#include <algorithm>
#include <cmath>
#include <complex>
#include <functional>
#include <limits>
#include <string>
#include <cstddef>
#include <vector>
 
#include "line/api/mam/map_moment.h"
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/linalg.h"
#include "line/util/lu.h"
#include "line/util/matrix.h"
 
namespace line {
namespace qsys {
 
/** Return value of qsys_mmapgk1, mirroring the MATLAB struct. */
template <class T>
struct MmapGk1Result {
    std::vector<T> lambdas;                 ///< per-type arrival rates
    T arrivalRate;                          ///< sum of the per-type rates
    T utilization;                          ///< rho = sum_k lambda_k E[S_k]
    std::vector<T> idleVector;              ///< v0, summing to 1 - rho
    std::vector<std::vector<T> > waitMoments;  ///< [type][order], E[Wk^j]
    std::vector<T> meanWaitingTime;         ///< per-type E[Wq]
    std::vector<T> meanSojournTime;         ///< per-type E[Wq] + E[S]
    T meanQueueLength;                      ///< E[N] by Little over all types
    std::vector<std::vector<T> > waitCDF;   ///< [type][point], P(Wk <= t)
    std::vector<T> waitPoints;              ///< the requested points
};
 
/**
 * A service law given by its transform, the reference's struct form
 * `struct('lst', @(s) ..., 'moments', [m1 m2 ...])`: the LST at a complex
 * argument and its raw moments, of which at least numWMoms + 1 are needed.
 */
template <class T>
struct MmapGk1Transform {
    std::function<std::complex<double>(std::complex<double>)> lst;
    std::vector<T> moments;
};
 
/**
 * One marked type's service law, in any of the three forms the reference
 * accepts: a Distribution object, a (D0, D1) MAP pair, or a transform struct.
 */
template <class T>
struct MmapGk1Service {
    enum Kind { DISTRIB, MAP, TRANSFORM };
    Kind kind = DISTRIB;
    lang::Distrib<T> dist;
    mam::Map<T> map;
    MmapGk1Transform<T> transform;
    MmapGk1Service() {}
    MmapGk1Service(const lang::Distrib<T>& d) : kind(DISTRIB), dist(d) {}
    MmapGk1Service(const mam::Map<T>& m) : kind(MAP), map(m) {}
    MmapGk1Service(const MmapGk1Transform<T>& t) : kind(TRANSFORM), transform(t) {}
};
 
namespace mmapgk1detail {
 
typedef std::complex<double> cdbl;
 
/** Binomial coefficient; declared ahead of the complex block that uses it. */
inline double binom(std::size_t n, std::size_t k);
 
/** Left null vector of G normalized to sum one. */
template <class T>
std::vector<T> stat_left_null(const Matrix<T>& G) {
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    const std::size_t n = G.rows();
    Matrix<T> A(n, n, zero);
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) A(j, i) = G(i, j);
    for (std::size_t j = 0; j < n; ++j) A(n - 1, j) = one;
    std::vector<T> b(n, zero);
    b[n - 1] = one;
    return line::solve(A, b);
}
 
/** Midpoint nodes and true CDF increments over the support of the law. */
template <class T>
void stieltjes_nodes(const lang::Distrib<T>& law, std::vector<T>& x, std::vector<T>& w,
                     std::size_t quad_nodes = 120) {
    // nGrid = max(quadNodes, 20) * 20 in the reference
    const std::size_t n_grid = (quad_nodes > 20 ? quad_nodes : 20) * 20;
    const T zero = num_traits<T>::from_int(0);
    const double mean = num_traits<T>::to_double(lang::dist_moment(law, 1u));
    double hi = mean * 60.0;
    const double var = num_traits<T>::to_double((lang::dist_moment(law, 2u) - lang::dist_moment(law, 1u) * lang::dist_moment(law, 1u)));
    if (std::isfinite(var) && var > 0.0) hi = std::max(hi, mean + 12.0 * std::sqrt(var));
    const double step = hi / static_cast<double>(n_grid);
    x.assign(n_grid, zero);
    w.assign(n_grid, zero);
    T mass = zero;
    T prev = lang::dist_cdf(law, num_traits<T>::from_double(0.0));
    for (std::size_t i = 0; i < n_grid; ++i) {
        const T right = num_traits<T>::from_double((i + 1) * step);
        const T cur = lang::dist_cdf(law, right);
        x[i] = num_traits<T>::from_double((i + 0.5) * step);
        w[i] = cur - prev;
        mass += w[i];
        prev = cur;
    }
    if (mass > zero)
        for (std::size_t i = 0; i < n_grid; ++i) w[i] = w[i] / mass;
}
 
/**
 * Raw moment E[S^j]; dist_moment already evaluates every family exactly.
 *
 * A HEAVY TAIL HAS NO MOMENT of high enough order, and dist_moment says so by
 * THROWING (a Pareto of shape <= j). Here that is not an error but an answer:
 * the moment recursion consumes M_j only from order j onwards, so an infinite
 * M_3 leaves E[Wq] finite and makes E[Wq^2] infinite, which is the truth about
 * such a queue. Letting the throw escape would lose the finite moments too.
 */
template <class T>
T raw_moment(const lang::Distrib<T>& law, std::size_t j) {
    try {
        return lang::dist_moment(law, static_cast<unsigned>(j));
    } catch (const NumericError&) {
        return num_traits<T>::from_double(std::numeric_limits<double>::infinity());
    }
}
 
/** Raw moment E[S^j] of a service law in any of its three forms. */
template <class T>
T raw_moment(const MmapGk1Service<T>& svc, std::size_t j) {
    switch (svc.kind) {
        case MmapGk1Service<T>::MAP:
            return mam::map_moment(svc.map, static_cast<unsigned>(j));
        case MmapGk1Service<T>::TRANSFORM:
            if (j > svc.transform.moments.size())
                throw InputError("qsys_mmapgk1: The service law was given as a transform with " +
                                 std::to_string(svc.transform.moments.size()) +
                                 " moments, but order " + std::to_string(j) +
                                 " is needed. Supply at least numWMoms+1 moments.");
            return svc.transform.moments[j - 1];
        default:
            return raw_moment(svc.dist, j);
    }
}
 
/** Kronecker product; C++ has no shared templated kron. */
template <class T>
Matrix<T> gk_kron(const Matrix<T>& A, const Matrix<T>& B) {
    Matrix<T> C(A.rows() * B.rows(), A.cols() * B.cols(), num_traits<T>::from_int(0));
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j)
            for (std::size_t p = 0; p < B.rows(); ++p)
                for (std::size_t q = 0; q < B.cols(); ++q)
                    C(i * B.rows() + p, j * B.cols() + q) = A(i, j) * B(p, q);
    return C;
}
 
/** Gaussian elimination with partial pivoting over the complex field. */
inline std::vector<cdbl> complex_solve(std::vector<std::vector<cdbl> > A, std::vector<cdbl> b);
 
/**
 * The matrix transform int_0^inf exp(U t) dF(t) for a Markovian (D0, D1)
 * service law. The density is the SCALAR beta exp(St) s0, so the integral is
 * exact on the Kronecker sum: int exp(Ut) x exp(St) dt = -(U (+) S)^-1.
 */
template <class T>
Matrix<T> matrix_lst_map(const mam::Map<T>& rep, const Matrix<T>& U) {
    const T zero = num_traits<T>::from_int(0);
    const std::size_t n = U.rows();
    const Matrix<T>& S = rep.D0;
    const std::size_t ms = S.rows();
    const std::vector<T> beta = mam::map_pie(rep);
    std::vector<T> s0(ms, zero);
    for (std::size_t i = 0; i < ms; ++i) {
        T r = zero;
        for (std::size_t j = 0; j < ms; ++j) r += S(i, j);
        s0[i] = -r;
    }
    Matrix<T> Ims = eye<T>(ms);
    Matrix<T> KS = gk_kron(U, Ims);
    const Matrix<T> In = eye<T>(n);
    const Matrix<T> IkS = gk_kron(In, S);
    for (std::size_t i = 0; i < KS.rows(); ++i)
        for (std::size_t j = 0; j < KS.cols(); ++j) KS(i, j) += IkS(i, j);
    Matrix<T> betaM(1, ms, zero);
    for (std::size_t i = 0; i < ms; ++i) betaM(0, i) = beta[i];
    Matrix<T> s0M(ms, 1, zero);
    for (std::size_t i = 0; i < ms; ++i) s0M(i, 0) = s0[i];
    const Matrix<T> L = gk_kron(In, betaM);
    Matrix<T> R = gk_kron(In, s0M);
    Matrix<T> X = matmul(inverse(KS), R);
    for (std::size_t i = 0; i < X.rows(); ++i)
        for (std::size_t j = 0; j < X.cols(); ++j) X(i, j) = -X(i, j);
    return matmul(L, X);
}
 
/**
 * Complex eigen-decomposition U = V diag(ev) V^-1 of a small real matrix, the
 * reference's `[Vd, Dd] = eig(U)`. LAPACK-free: the eigenvalues are the roots
 * of the characteristic polynomial (Faddeev-LeVerrier, then Durand-Kerner),
 * the eigenvectors come from inverse iteration. `rcond` is the reciprocal
 * 1-norm condition number of V, which is what MATLAB's rcond(Vd) estimates;
 * it is 0 when V could not be inverted.
 */
inline void eig_decomp(const std::vector<std::vector<double> >& U, std::vector<cdbl>& ev,
                       std::vector<std::vector<cdbl> >& V, std::vector<std::vector<cdbl> >& Vinv,
                       double& rcond) {
    const std::size_t n = U.size();
    ev.assign(n, cdbl(0.0, 0.0));
    V.assign(n, std::vector<cdbl>(n, cdbl(0.0, 0.0)));
    Vinv.assign(n, std::vector<cdbl>(n, cdbl(0.0, 0.0)));
    rcond = 0.0;
    if (n == 0) return;
    // characteristic polynomial det(tau I - U) = tau^n + c1 tau^(n-1) + ... + cn
    std::vector<double> c(n + 1, 0.0);
    c[0] = 1.0;
    std::vector<std::vector<double> > Mk(n, std::vector<double>(n, 0.0));
    for (std::size_t i = 0; i < n; ++i) Mk[i][i] = 1.0;
    for (std::size_t k = 1; k <= n; ++k) {
        std::vector<std::vector<double> > TM(n, std::vector<double>(n, 0.0));
        double tr = 0.0;
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) {
                double acc = 0.0;
                for (std::size_t l = 0; l < n; ++l) acc += U[i][l] * Mk[l][j];
                TM[i][j] = acc;
                if (i == j) tr += acc;
            }
        c[k] = -tr / static_cast<double>(k);
        Mk = TM;
        for (std::size_t i = 0; i < n; ++i) Mk[i][i] += c[k];
    }
    // Durand-Kerner
    {
        const cdbl seed(0.4, 0.9);
        cdbl pw(1.0, 0.0);
        for (std::size_t i = 0; i < n; ++i) {
            ev[i] = pw;
            pw *= seed;
        }
        double scale = 0.0;
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) scale = std::max(scale, std::fabs(U[i][j]));
        if (scale > 0.0)
            for (std::size_t i = 0; i < n; ++i) ev[i] *= scale;
        for (int it = 0; it < 2000; ++it) {
            double move = 0.0;
            for (std::size_t i = 0; i < n; ++i) {
                cdbl den(1.0, 0.0);
                for (std::size_t j = 0; j < n; ++j)
                    if (j != i) den *= (ev[i] - ev[j]);
                if (std::abs(den) < 1e-300) continue;
                cdbl pv(0.0, 0.0);
                for (std::size_t q = 0; q <= n; ++q) pv = pv * ev[i] + c[q];
                const cdbl step = pv / den;
                ev[i] -= step;
                move = std::max(move, std::abs(step));
            }
            if (move < 1e-15 * std::max(1.0, scale)) break;
        }
    }
    // eigenvectors by inverse iteration on the slightly shifted eigenvalue
    double umax = 0.0;
    for (std::size_t i = 0; i < n; ++i) umax = std::max(umax, std::abs(ev[i]));
    for (std::size_t i = 0; i < n; ++i) {
        const cdbl shift = ev[i] + cdbl(1e-10 * std::max(1.0, umax), 1e-10 * std::max(1.0, umax));
        std::vector<cdbl> v(n);
        for (std::size_t q = 0; q < n; ++q) v[q] = cdbl(1.0 / std::sqrt(static_cast<double>(n)), 0.0);
        for (int it = 0; it < 3; ++it) {
            std::vector<std::vector<cdbl> > A(n, std::vector<cdbl>(n, cdbl(0.0, 0.0)));
            for (std::size_t p = 0; p < n; ++p) {
                for (std::size_t q = 0; q < n; ++q) A[p][q] = cdbl(U[p][q], 0.0);
                A[p][p] -= shift;
            }
            std::vector<cdbl> w = complex_solve(A, v);
            double nrm = 0.0;
            for (std::size_t q = 0; q < n; ++q) nrm += std::norm(w[q]);
            nrm = std::sqrt(nrm);
            if (!(nrm > 0.0) || !std::isfinite(nrm)) break;
            for (std::size_t q = 0; q < n; ++q) v[q] = w[q] / nrm;
        }
        for (std::size_t q = 0; q < n; ++q) V[q][i] = v[q];
    }
    // V^-1 column by column, and rcond = 1 / (||V||_1 ||V^-1||_1)
    double normV = 0.0, normW = 0.0;
    for (std::size_t j = 0; j < n; ++j) {
        std::vector<cdbl> ej(n, cdbl(0.0, 0.0));
        ej[j] = cdbl(1.0, 0.0);
        const std::vector<cdbl> col = complex_solve(V, ej);
        double cs = 0.0, vs = 0.0;
        for (std::size_t i = 0; i < n; ++i) {
            Vinv[i][j] = col[i];
            cs += std::abs(col[i]);
            vs += std::abs(V[i][j]);
        }
        normW = std::max(normW, cs);
        normV = std::max(normV, vs);
    }
    if (std::isfinite(normW) && normW > 0.0 && normV > 0.0) rcond = 1.0 / (normV * normW);
}
 
/** real(V diag(g) V^-1), the eigen route of the matrix transform. */
template <class T>
Matrix<T> matrix_lst_eig(std::size_t n, const std::vector<std::vector<cdbl> >& V,
                         const std::vector<std::vector<cdbl> >& Vinv, const std::vector<cdbl>& g) {
    Matrix<T> F(n, n, num_traits<T>::from_int(0));
    for (std::size_t i = 0; i < n; ++i)
        for (std::size_t j = 0; j < n; ++j) {
            cdbl acc(0.0, 0.0);
            for (std::size_t q = 0; q < n; ++q) acc += V[i][q] * g[q] * Vinv[q][j];
            F(i, j) = num_traits<T>::from_double(acc.real());
        }
    return F;
}
 
/** The matrix transform int_0^inf exp(U t) dF(t), for a Distribution object. */
template <class T>
Matrix<T> matrix_lst(const lang::Distrib<T>& law, const Matrix<T>& U, std::size_t quad_nodes = 120) {
    const T zero = num_traits<T>::from_int(0);
    const std::size_t n = U.rows();
    if (lang::process_is_markovian(law.type)) return matrix_lst_map(lang::dist_to_map(law), U);
    if (law.type == lang::ProcessType::DET) {
        Matrix<T> Ud = U;
        const T d = lang::dist_moment(law, 1u);
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) Ud(i, j) = Ud(i, j) * d;
        return line::expm(Ud);
    }
    // The reference's eigen-decomposition route: F = real(V diag(lst(-u)) V^-1)
    // when V is well conditioned and every transform value is finite, else the
    // Stieltjes sum over the CDF increments.
    {
        std::vector<std::vector<double> > Ud(n, std::vector<double>(n, 0.0));
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j) Ud[i][j] = num_traits<T>::to_double(U(i, j));
        std::vector<cdbl> ev;
        std::vector<std::vector<cdbl> > V, Vinv;
        double rcond = 0.0;
        eig_decomp(Ud, ev, V, Vinv, rcond);
        bool ok = rcond > 1e-12;
        std::vector<cdbl> g(n, cdbl(0.0, 0.0));
        if (ok) {
            for (std::size_t i = 0; i < n; ++i) {
                try {
                    g[i] = lang::dist_lst(law, -ev[i]);
                } catch (const Error&) {
                    ok = false;
                    break;
                }
                if (!std::isfinite(g[i].real()) || !std::isfinite(g[i].imag())) {
                    ok = false;
                    break;
                }
            }
        }
        if (ok) return matrix_lst_eig<T>(n, V, Vinv, g);
    }
    std::vector<T> x, w;
    stieltjes_nodes(law, x, w, quad_nodes);
    Matrix<T> F(n, n, zero);
    for (std::size_t i = 0; i < x.size(); ++i) {
        Matrix<T> Ut = U;
        for (std::size_t p = 0; p < n; ++p)
            for (std::size_t q = 0; q < n; ++q) Ut(p, q) = Ut(p, q) * x[i];
        const Matrix<T> E = line::expm(Ut);
        for (std::size_t p = 0; p < n; ++p)
            for (std::size_t q = 0; q < n; ++q) F(p, q) += w[i] * E(p, q);
    }
    return F;
}
 
/** The matrix transform for a service law in any of its three forms. */
template <class T>
Matrix<T> matrix_lst(const MmapGk1Service<T>& svc, const Matrix<T>& U, std::size_t quad_nodes) {
    switch (svc.kind) {
        case MmapGk1Service<T>::MAP:
            return matrix_lst_map(svc.map, U);
        case MmapGk1Service<T>::TRANSFORM: {
            // F = real(V diag(lst(-u)) V^-1), as the reference does for a struct law
            const std::size_t n = U.rows();
            std::vector<std::vector<double> > Ud(n, std::vector<double>(n, 0.0));
            for (std::size_t i = 0; i < n; ++i)
                for (std::size_t j = 0; j < n; ++j) Ud[i][j] = num_traits<T>::to_double(U(i, j));
            std::vector<cdbl> ev;
            std::vector<std::vector<cdbl> > V, Vinv;
            double rcond = 0.0;
            eig_decomp(Ud, ev, V, Vinv, rcond);
            std::vector<cdbl> g(n);
            for (std::size_t i = 0; i < n; ++i) g[i] = svc.transform.lst(-ev[i]);
            return matrix_lst_eig<T>(n, V, Vinv, g);
        }
        default:
            return matrix_lst(svc.dist, U, quad_nodes);
    }
}
 
/** Gaussian elimination with partial pivoting over the complex field. */
inline std::vector<cdbl> complex_solve(std::vector<std::vector<cdbl> > A, std::vector<cdbl> b) {
    const std::size_t n = b.size();
    for (std::size_t col = 0; col < n; ++col) {
        std::size_t piv = col;
        double best = std::abs(A[col][col]);
        for (std::size_t r = col + 1; r < n; ++r) {
            if (std::abs(A[r][col]) > best) { best = std::abs(A[r][col]); piv = r; }
        }
        if (piv != col) { std::swap(A[piv], A[col]); std::swap(b[piv], b[col]); }
        for (std::size_t r = col + 1; r < n; ++r) {
            const cdbl f = A[r][col] / A[col][col];
            for (std::size_t c = col; c < n; ++c) A[r][c] -= f * A[col][c];
            b[r] -= f * b[col];
        }
    }
    std::vector<cdbl> x(n);
    for (std::size_t row = n; row-- > 0;) {
        cdbl s = b[row];
        for (std::size_t c = row + 1; c < n; ++c) s -= A[row][c] * x[c];
        x[row] = s / A[row][row];
    }
    return x;
}
 
/**
 * Scalar Laplace-Stieltjes transform of a service law at a complex argument.
 *
 * One line, because the LANGUAGE LAYER owns this: dist_lst carries a complex
 * overload beside the real one, with the closed forms, the phase-type solve and
 * the CDF-increment fallback already tiered. A second implementation here would
 * be a second thing to keep true.
 */
template <class T>
cdbl scalar_lst(const lang::Distrib<T>& law, cdbl s) {
    return lang::dist_lst(law, s);
}
 
/** Scalar LST of a service law in any of its three forms. */
template <class T>
cdbl scalar_lst(const MmapGk1Service<T>& svc, cdbl s) {
    switch (svc.kind) {
        case MmapGk1Service<T>::MAP: {
            // alpha (s I - D0)^-1 (-D0 e), with alpha = map_pie
            const Matrix<T>& D0s = svc.map.D0;
            const std::size_t n = D0s.rows();
            const std::vector<T> alpha = mam::map_pie(svc.map);
            std::vector<std::vector<cdbl> > A(n, std::vector<cdbl>(n, cdbl(0.0, 0.0)));
            std::vector<cdbl> b(n, cdbl(0.0, 0.0));
            for (std::size_t i = 0; i < n; ++i) {
                double rs = 0.0;
                for (std::size_t j = 0; j < n; ++j) {
                    A[i][j] = cdbl(-num_traits<T>::to_double(D0s(i, j)), 0.0);
                    rs += num_traits<T>::to_double(D0s(i, j));
                }
                A[i][i] += s;
                b[i] = cdbl(-rs, 0.0);
            }
            const std::vector<cdbl> x = complex_solve(A, b);
            cdbl g(0.0, 0.0);
            for (std::size_t i = 0; i < n; ++i) g += num_traits<T>::to_double(alpha[i]) * x[i];
            return g;
        }
        case MmapGk1Service<T>::TRANSFORM:
            return svc.transform.lst(s);
        default:
            return lang::dist_lst(svc.dist, s);
    }
}
 
 
/** E[exp(-s Wk)] at a complex argument. */
template <class T>
cdbl wait_lst(const Matrix<T>& D0, const std::vector<Matrix<T> >& Dk,
              const std::vector<MmapGk1Service<T> >& svc, const std::vector<T>& v0,
              const std::vector<T>& lambdas, std::size_t k, cdbl s) {
    const std::size_t ma = D0.rows();
    std::vector<std::vector<cdbl> > M(ma, std::vector<cdbl>(ma, cdbl(0.0, 0.0)));
    for (std::size_t i = 0; i < ma; ++i)
        for (std::size_t j = 0; j < ma; ++j) M[i][j] = cdbl(num_traits<T>::to_double(D0(i, j)), 0.0);
    for (std::size_t i = 0; i < ma; ++i) M[i][i] += s;
    for (std::size_t q = 0; q < Dk.size(); ++q) {
        const cdbl g = scalar_lst(svc[q], s);
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j)
                M[i][j] += num_traits<T>::to_double(Dk[q](i, j)) * g;
    }
    // solve f M = s v0, i.e. M' f' = (s v0)'
    std::vector<std::vector<cdbl> > A(ma, std::vector<cdbl>(ma, cdbl(0.0, 0.0)));
    std::vector<cdbl> b(ma, cdbl(0.0, 0.0));
    for (std::size_t i = 0; i < ma; ++i) {
        for (std::size_t j = 0; j < ma; ++j) A[i][j] = M[j][i];
        b[i] = s * num_traits<T>::to_double(v0[i]);
    }
    const std::vector<cdbl> f = complex_solve(A, b);
    cdbl out(0.0, 0.0);
    for (std::size_t i = 0; i < ma; ++i) {
        double rowsum = 0.0;
        for (std::size_t j = 0; j < ma; ++j) rowsum += num_traits<T>::to_double(Dk[k](i, j));
        out += f[i] * rowsum;
    }
    return out / num_traits<T>::to_double(lambdas[k]);
}
 
/** Abate-Whitt Euler inversion of the type-k waiting time CDF. */
template <class T>
double euler_invert(const Matrix<T>& D0, const std::vector<Matrix<T> >& Dk,
                    const std::vector<MmapGk1Service<T> >& svc, const std::vector<T>& v0,
                    const std::vector<T>& lambdas, std::size_t k, double t) {
    if (t <= 0.0) return wait_lst(D0, Dk, svc, v0, lambdas, k, cdbl(1e12, 0.0)).real();
    const double A = 18.4;
    const std::size_t nE = 15, mE = 11;
    const double u = std::exp(A / 2) / t;
    const double x = A / (2 * t);
    std::vector<double> terms(nE + mE + 1, 0.0);
    terms[0] = wait_lst(D0, Dk, svc, v0, lambdas, k, cdbl(x, 0.0)).real() / x / 2.0;
    for (std::size_t j = 1; j <= nE + mE; ++j) {
        const cdbl s(x, M_PI * static_cast<double>(j) / t);
        terms[j] = ((j % 2 == 0) ? 1.0 : -1.0) * (wait_lst(D0, Dk, svc, v0, lambdas, k, s) / s).real();
    }
    std::vector<double> partial(terms.size(), 0.0);
    double run = 0.0;
    for (std::size_t j = 0; j < terms.size(); ++j) { run += terms[j]; partial[j] = run; }
    double F = 0.0;
    for (std::size_t j = 0; j <= mE; ++j) F += binom(mE, j) / std::pow(2.0, static_cast<double>(mE)) * partial[nE + j];
    F = u * F;
    return std::min(std::max(F, 0.0), 1.0);
}
 
 
inline double binom(std::size_t n, std::size_t k) {
    double r = 1.0;
    for (std::size_t i = 1; i <= k; ++i) r = r * static_cast<double>(n - k + i) / static_cast<double>(i);
    return std::floor(r + 0.5);
}
 
}  // namespace mmapgk1detail
 
/**
 * MMAP[K]/G[K]/1 FCFS, per type.
 *
 * @param MMAP       LINE convention {D0, D1, D^(1), ..., D^(K)}, D1 = sum_k D^(k)
 * @param svc        K service laws, one per marked type, each a Distribution
 *                   object, a (D0, D1) MAP pair or a transform struct
 * @param w_points   times at which to evaluate the per-type waiting time CDF
 * @param num_w_moms how many per-type waiting time moments to return
 * @param tol        fixed point tolerance on U
 * @param iter_max   fixed point iteration cap (the reference's 10000)
 * @param quad_nodes the reference's quadNodes (default 120): the Stieltjes
 *                   fallback uses max(quadNodes, 20) * 20 CDF increments
 */
template <class T>
MmapGk1Result<T> qsys_mmapgk1(const std::vector<Matrix<T> >& MMAP,
                              const std::vector<MmapGk1Service<T> >& svc,
                              const std::vector<T>& w_points, std::size_t num_w_moms,
                              double tol, std::size_t iter_max, std::size_t quad_nodes) {
    static_assert(num_traits<T>::has_transcendental,
                  "qsys_mmapgk1 requires transcendental arithmetic");
    using namespace mmapgk1detail;
    const T zero = num_traits<T>::from_int(0), one = num_traits<T>::from_int(1);
    if (MMAP.size() < 3) throw InputError("qsys_mmapgk1: the MMAP must carry a marked block");
    const std::size_t K = MMAP.size() - 2;
    if (svc.size() != K) throw InputError("qsys_mmapgk1: one service law per marked type");
    const Matrix<T>& D0 = MMAP[0];
    const std::size_t ma = D0.rows();
    std::vector<Matrix<T> > Dk;
    Matrix<T> Dsum = D0;
    for (std::size_t k = 0; k < K; ++k) {
        Dk.push_back(MMAP[k + 2]);
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) Dsum(i, j) += Dk[k](i, j);
    }
 
    const std::vector<T> theta = stat_left_null(Dsum);
    std::vector<T> lambdas(K, zero), mean_s(K, zero);
    T rho = zero;
    for (std::size_t k = 0; k < K; ++k) {
        T lam = zero;
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j) lam += theta[i] * Dk[k](i, j);
        lambdas[k] = lam;
        mean_s[k] = raw_moment(svc[k], 1);
        rho += lam * mean_s[k];
    }
    if (rho >= one) throw InputError("qsys_mmapgk1: load rho must be strictly less than 1");
 
    Matrix<T> U = D0;
    for (std::size_t it = 0; it < iter_max; ++it) {
        Matrix<T> Unew = D0;
        for (std::size_t k = 0; k < K; ++k) {
            const Matrix<T> Fk = matrix_lst(svc[k], U, quad_nodes);
            const Matrix<T> P = matmul(Dk[k], Fk);
            for (std::size_t i = 0; i < ma; ++i)
                for (std::size_t j = 0; j < ma; ++j) Unew(i, j) += P(i, j);
        }
        double diff = 0.0;
        for (std::size_t i = 0; i < ma; ++i)
            for (std::size_t j = 0; j < ma; ++j)
                diff = std::max(diff, std::fabs(num_traits<T>::to_double(Unew(i, j) - U(i, j))));
        U = Unew;
        if (diff <= tol) break;
    }
    // U e = 0 EXACTLY: a property of the fixed point, not of the iterate, and
    // the residue at the tolerance above would move the stationary solve.
    for (std::size_t i = 0; i < ma; ++i) {
        T rowsum = zero;
        for (std::size_t j = 0; j < ma; ++j) rowsum += U(i, j);
        U(i, i) -= rowsum;
    }
    std::vector<T> v0 = stat_left_null(U);
    for (std::size_t i = 0; i < ma; ++i) v0[i] = v0[i] * (one - rho);
 
    // Moment recursion, see the header comment
    std::vector<Matrix<T> > Mder;
    for (std::size_t j = 0; j <= num_w_moms + 1; ++j) {
        Matrix<T> Mj(ma, ma, zero);
        if (j == 0) {
            Mj = D0;
            for (std::size_t k = 0; k < K; ++k)
                for (std::size_t i = 0; i < ma; ++i)
                    for (std::size_t q = 0; q < ma; ++q) Mj(i, q) += Dk[k](i, q);
        } else {
            for (std::size_t k = 0; k < K; ++k) {
                const T mom = raw_moment(svc[k], j);
                const double sign = (j % 2 == 0) ? 1.0 : -1.0;
                for (std::size_t i = 0; i < ma; ++i)
                    for (std::size_t q = 0; q < ma; ++q)
                        Mj(i, q) += num_traits<T>::from_double(sign) * Dk[k](i, q) * mom;
            }
            if (j == 1)
                for (std::size_t i = 0; i < ma; ++i) Mj(i, i) += one;
        }
        Mder.push_back(Mj);
    }
    const std::vector<T> e(ma, one);
    T denom = zero;
    {
        const std::vector<T> v = mulvec(Mder[1], e);
        for (std::size_t i = 0; i < ma; ++i) denom += theta[i] * v[i];
    }
    std::vector<std::vector<T> > fder;
    fder.push_back(theta);
    // [M_0, e] with f_j^p e = 0 pins the particular solution; solved in the
    // least squares sense, since the system is one equation over-determined.
    Matrix<T> Abase(ma + 1, ma, zero);
    for (std::size_t i = 0; i < ma; ++i)
        for (std::size_t j = 0; j < ma; ++j) Abase(j, i) = Mder[0](i, j);
    for (std::size_t i = 0; i < ma; ++i) Abase(ma, i) = one;
    Matrix<T> AtA(ma, ma, zero);
    for (std::size_t i = 0; i < ma; ++i)
        for (std::size_t j = 0; j < ma; ++j) {
            T s = zero;
            for (std::size_t q = 0; q < ma + 1; ++q) s += Abase(q, i) * Abase(q, j);
            AtA(i, j) = s;
        }
    const Matrix<T> AtAinv = inverse(AtA);
    for (std::size_t j = 1; j <= num_w_moms; ++j) {
        std::vector<T> rhs(ma, zero);
        if (j == 1) for (std::size_t i = 0; i < ma; ++i) rhs[i] = v0[i];
        for (std::size_t i = 0; i < j; ++i) {
            const double c = binom(j, i);
            for (std::size_t q = 0; q < ma; ++q) {
                T acc = zero;
                for (std::size_t p = 0; p < ma; ++p) acc += fder[i][p] * Mder[j - i](p, q);
                rhs[q] -= num_traits<T>::from_double(c) * acc;
            }
        }
        std::vector<T> rhsAug(ma + 1, zero);
        for (std::size_t q = 0; q < ma; ++q) rhsAug[q] = rhs[q];
        std::vector<T> Atb(ma, zero);
        for (std::size_t i = 0; i < ma; ++i) {
            T s = zero;
            for (std::size_t q = 0; q < ma + 1; ++q) s += Abase(q, i) * rhsAug[q];
            Atb[i] = s;
        }
        const std::vector<T> fp = mulvec(AtAinv, Atb);
        T acc2 = zero;
        for (std::size_t i = 0; i < j; ++i) {
            const std::vector<T> v = mulvec(Mder[j + 1 - i], e);
            T inner = zero;
            for (std::size_t p = 0; p < ma; ++p) inner += fder[i][p] * v[p];
            acc2 += num_traits<T>::from_double(binom(j + 1, i)) * inner;
        }
        T fpM1e = zero;
        {
            const std::vector<T> v = mulvec(Mder[1], e);
            for (std::size_t p = 0; p < ma; ++p) fpM1e += fp[p] * v[p];
        }
        const T cfree = (-acc2 / num_traits<T>::from_int(static_cast<long>(j + 1)) - fpM1e) / denom;
        std::vector<T> fj(ma, zero);
        for (std::size_t p = 0; p < ma; ++p) fj[p] = fp[p] + cfree * theta[p];
        fder.push_back(fj);
    }
 
    MmapGk1Result<T> r;
    r.lambdas = lambdas;
    r.arrivalRate = zero;
    for (std::size_t k = 0; k < K; ++k) r.arrivalRate += lambdas[k];
    r.utilization = rho;
    r.idleVector = v0;
    r.waitMoments.assign(K, std::vector<T>(num_w_moms, zero));
    r.meanWaitingTime.assign(K, zero);
    r.meanSojournTime.assign(K, zero);
    r.meanQueueLength = zero;
    for (std::size_t k = 0; k < K; ++k) {
        for (std::size_t j = 1; j <= num_w_moms; ++j) {
            const std::vector<T> v = mulvec(Dk[k], e);
            T inner = zero;
            for (std::size_t p = 0; p < ma; ++p) inner += fder[j][p] * v[p];
            const double sign = (j % 2 == 0) ? 1.0 : -1.0;
            r.waitMoments[k][j - 1] = num_traits<T>::from_double(sign) * inner / lambdas[k];
        }
        r.meanWaitingTime[k] = num_w_moms >= 1 ? r.waitMoments[k][0] : zero;
        r.meanSojournTime[k] = r.meanWaitingTime[k] + mean_s[k];
        r.meanQueueLength += lambdas[k] * r.meanSojournTime[k];
    }
 
    r.waitPoints = w_points;
    if (!w_points.empty()) {
        r.waitCDF.assign(K, std::vector<T>(w_points.size(), zero));
        for (std::size_t k = 0; k < K; ++k) {
            for (std::size_t it = 0; it < w_points.size(); ++it) {
                const double t = num_traits<T>::to_double(w_points[it]);
                r.waitCDF[k][it] = num_traits<T>::from_double(
                    euler_invert(D0, Dk, svc, v0, lambdas, k, t));
            }
        }
    }
    return r;
}
 
 
/** qsys_mmapgk1 on service variants with the reference's quadNodes = 120. */
template <class T>
MmapGk1Result<T> qsys_mmapgk1(const std::vector<Matrix<T> >& MMAP,
                              const std::vector<MmapGk1Service<T> >& svc,
                              const std::vector<T>& w_points, std::size_t num_w_moms,
                              double tol, std::size_t iter_max) {
    return qsys_mmapgk1(MMAP, svc, w_points, num_w_moms, tol, iter_max,
                        static_cast<std::size_t>(120));
}
 
/** qsys_mmapgk1 on Distribution objects. */
template <class T>
MmapGk1Result<T> qsys_mmapgk1(const std::vector<Matrix<T> >& MMAP,
                              const std::vector<lang::Distrib<T> >& svc,
                              const std::vector<T>& w_points, std::size_t num_w_moms,
                              double tol, std::size_t iter_max) {
    std::vector<MmapGk1Service<T> > laws;
    laws.reserve(svc.size());
    for (std::size_t k = 0; k < svc.size(); ++k) laws.push_back(MmapGk1Service<T>(svc[k]));
    return qsys_mmapgk1(MMAP, laws, w_points, num_w_moms, tol, iter_max,
                        static_cast<std::size_t>(120));
}
 
/** qsys_mmapgk1 on Distribution objects with the reference's defaults. */
template <class T>
MmapGk1Result<T> qsys_mmapgk1(const std::vector<Matrix<T> >& MMAP,
                              const std::vector<lang::Distrib<T> >& svc) {
    return qsys_mmapgk1(MMAP, svc, std::vector<T>(), static_cast<std::size_t>(3), 1e-12,
                        static_cast<std::size_t>(10000));
}
 
/** qsys_mmapgk1 on service variants with the reference's defaults. */
template <class T>
MmapGk1Result<T> qsys_mmapgk1(const std::vector<Matrix<T> >& MMAP,
                              const std::vector<MmapGk1Service<T> >& svc) {
    return qsys_mmapgk1(MMAP, svc, std::vector<T>(), static_cast<std::size_t>(3), 1e-12,
                        static_cast<std::size_t>(10000), static_cast<std::size_t>(120));
}
 
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_MMAPGK1_H