solver_mam_reneging.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_MAM_SOLVER_MAM_RENEGING_H
#define LINE_SOLVERS_MAM_SOLVER_MAM_RENEGING_H
 
/**
 * @file
 * @ingroup line_solvers
 * The reneging route of `solver_mam_retrial.m`: the MAP/M/s+G queue whose
 * waiting jobs abandon after a generally distributed patience, after O. Gursoy,
 * K. A. Mehr and N. Akar, "The MAP/M/s + G Call Center Model with Generally
 * Distributed Patience Times". Gate: `mam_reneging_applicable.m`.
 *
 * THE FLUID IS THE VIRTUAL WAITING TIME. The background chain is (busy servers,
 * MAP phase), s+1 blocks of MAP phases; the last block has drift +1 and stands
 * for the work an admitted arrival adds, every other block drift -1. The
 * patience CDF is quantized into `MamOptions::mapmsg_quantization` regimes over
 * [0, T], T the first mean * 1.25^k at which the CDF reaches 1 - 1e-4. An arrival
 * that finds the virtual wait in regime k abandons with the patience CDF at the
 * regime's representative wait (the midpoint, and T for the tail [T, inf)),
 * which is the abandonment of a FCFS queue because the virtual wait IS the wait
 * that arrival would see.
 *
 * MATLAB solves the fluid queue with the MAPMsG library's `MRMFQSolver`; this
 * port, like the JAR's, solves the same queue with `mfq_multiregime` (Kankaya
 * and Akar, the method MRMFQSolver implements) and reads the per-regime state
 * masses off its distribution function: the mass on regime d is
 * P(X < B_d) - P(X <= B_{d-1}) and the atom at zero is P(X <= 0). On M/M/3+M,
 * MAP/M/3+M and M/M/4+Erlang-2 this construction and the native-python MRMFQ
 * port coincide to 1e-14.
 *
 * THE PATIENCE MUST BE MARKOVIAN, and that is refused by name rather than
 * fitted. The regime grid, the horizon and E[min(v, P)] are all read off a
 * (D0, D1) pair. A Det, Gamma, Lognormal, Pareto, Uniform or Weibull patience
 * carries none in the C++ model layer, and the JAR's `refreshImpatience` stores
 * none either. MATLAB stores `getProcess()` for these laws, which is the raw
 * parameter list and not a generator, and the reneging route then fails inside
 * `schur` on NaN, so there is no reference behaviour to reproduce.
 *
 * DOUBLE ONLY. `mfq_multiregime` is double only (an ordered real Schur form,
 * LAPACK), so an exact or multiprecision instantiation refuses by name instead
 * of carrying a double-accurate answer under a wider type.
 */
 
#include <cstddef>
#include <string>
#include <type_traits>
#include <vector>
 
#include "line/api/mam/map_cdf.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/mfq_multiregime.h"
#include "line/lang/distribution.h"
#include "line/lang/qn/network_struct.h"
#include "line/solvers/mam/mam_types.h"
#include "line/solvers/mam/solver_mam_basic.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 mam {
 
/** What `mam_reneging_applicable.m` returns: the station it found, or why it found none. */
struct MamRenegingInfo {
    bool ok = false;
    std::size_t source = 0;   ///< 1-based station index of the Source
    std::size_t station = 0;  ///< 1-based station index of the reneging queue
    std::size_t cls = 1;      ///< 1-based class index; the reference is single-class
    int nservers = 0;         ///< s
    std::string why;          ///< the reference's `errorMsg`, empty when ok
};
 
namespace reneging_detail {
 
/** Does station `ist` (1-based) declare RENEGING for class `r` (1-based)? */
template <class T>
bool declares(const qn::NetworkStruct<T>& L, std::size_t ist, std::size_t r) {
    const qn::Station<T>& st = L.stations[ist - 1];
    return r - 1 < st.impatience.size() && st.impatience[r - 1] == lang::ImpatienceType::RENEGING;
}
 
/** The patience law at (ist, r) when it is a usable (D0, D1) pair; null otherwise. */
template <class T>
const lang::Distrib<T>* markovian_patience(const qn::NetworkStruct<T>& L, std::size_t ist,
                                          std::size_t r) {
    if (!declares(L, ist, r)) return nullptr;
    const qn::Station<T>& st = L.stations[ist - 1];
    if (r - 1 >= st.patience.size()) return nullptr;
    const lang::Distrib<T>& p = st.patience[r - 1];
    if (p.disabled || !p.has_map()) return nullptr;
    if (p.type == lang::ProcessType::IMMEDIATE || p.type == lang::ProcessType::DISABLED ||
        !basic_detail::is_markovian_type(p.type))
        return nullptr;
    return &p;
}
 
/** Matrix<T> to Matrix<double>. */
template <class T>
Matrix<double> to_double(const Matrix<T>& A) {
    Matrix<double> B(A.rows(), A.cols(), 0.0);
    for (std::size_t i = 0; i < A.rows(); ++i)
        for (std::size_t j = 0; j < A.cols(); ++j) B(i, j) = num_traits<T>::to_double(A(i, j));
    return B;
}
 
inline void add_block(Matrix<double>& A, std::size_t r0, std::size_t c0, const Matrix<double>& B,
                      double scale) {
    for (std::size_t i = 0; i < B.rows(); ++i)
        for (std::size_t j = 0; j < B.cols(); ++j) A(r0 + i, c0 + j) += scale * B(i, j);
}
 
/** The station metrics of the MAP/M/s+G queue. */
struct RenegingMetrics {
    double tput = 0.0, util = 0.0, qlen = 0.0, respt = 0.0;
};
 
/**
 * `solveReneging` + `convertPatienceToRegimes` of `solver_mam_retrial.m`, as
 * fixed in 2c0848d76f (stationary arrival rate, the regime grid, the QLen
 * formula), in double.
 *
 * @param C, D   the arrival MAP
 * @param mu     the exponential service rate
 * @param s      the server count
 * @param P0, P1 the patience as a (D0, D1) pair
 * @param Qz     the regime count, at least 2
 */
inline RenegingMetrics solve(const Matrix<double>& C, const Matrix<double>& D, double mu, int s,
                             const Matrix<double>& P0, const Matrix<double>& P1, std::size_t Qz) {
    Map<double> pat;
    pat.D0 = P0;
    pat.D1 = P1;
 
    // ---- the regimes of the patience CDF ------------------------------------
    double Th = map_mean(pat);
    while (map_cdf(pat, std::vector<double>(1, Th))[0] < 1.0 - 1e-4) Th *= 1.25;
    std::vector<double> BL(Qz);
    for (std::size_t k = 0; k < Qz; ++k)
        BL[k] = Th * static_cast<double>(k) / static_cast<double>(Qz - 1);  // linspace(0, T, Qz)
    BL[Qz - 1] = Th;
    std::vector<double> vrep(Qz);
    for (std::size_t k = 0; k + 1 < Qz; ++k) vrep[k] = (BL[k] + BL[k + 1]) / 2.0;
    vrep[Qz - 1] = BL[Qz - 1];
    const std::vector<double> ga = map_cdf(pat, vrep);  // ga[k] is the reference's ga(k+2)
 
    // ---- the background generators -----------------------------------------
    const std::size_t m = C.rows();
    const std::size_t S = static_cast<std::size_t>(s);
    const std::size_t n = (S + 1) * m;
    const Matrix<double> I = eye<double>(m);
    Matrix<double> Qy0(n, n, 0.0);
    for (std::size_t row = 1; row <= S + 1; ++row) {
        const std::size_t rs = (row - 1) * m;
        const double busy = static_cast<double>(row - 1) * mu;
        if (row == 1) {
            add_block(Qy0, rs, rs, C, 1.0);
            add_block(Qy0, rs, rs + m, D, 1.0);
        } else if (row == S + 1) {
            add_block(Qy0, rs, rs, I, -busy);
            add_block(Qy0, rs, rs - m, I, busy);
        } else {
            add_block(Qy0, rs, rs, C, 1.0);
            add_block(Qy0, rs, rs, I, -busy);
            add_block(Qy0, rs, rs - m, I, busy);
            add_block(Qy0, rs, rs + m, D, 1.0);
        }
    }
    // Above level zero every server is busy: an arrival is admitted to the
    // all-busy block with probability 1 - ga[k] and abandons otherwise.
    std::vector<Matrix<double>> Qy;
    Qy.reserve(Qz);
    for (std::size_t k = 0; k < Qz; ++k) {
        Matrix<double> Qk(n, n, 0.0);
        const std::size_t rs = (S - 1) * m;
        add_block(Qk, rs, rs, D, ga[k]);
        add_block(Qk, rs, rs, C, 1.0);
        add_block(Qk, rs, rs + m, D, 1.0 - ga[k]);
        const std::size_t rs2 = S * m;
        add_block(Qk, rs2, rs2, I, -static_cast<double>(S) * mu);
        add_block(Qk, rs2, rs2 - m, I, static_cast<double>(S) * mu);
        Qy.push_back(Qk);
    }
    std::vector<double> Rd(n, -1.0);
    for (std::size_t i = n - m; i < n; ++i) Rd[i] = 1.0;
    std::vector<std::vector<double>> R(Qz, Rd), Rt(Qz + 1, Rd);
    std::vector<Matrix<double>> Qt;
    Qt.reserve(Qz + 1);
    Qt.push_back(Qy0);
    for (std::size_t k = 0; k < Qz; ++k) Qt.push_back(Qy[k]);
    // thresholds: B(1) dropped and 1e7 appended, as the reference's solveReneging
    std::vector<double> B(Qz);
    for (std::size_t k = 0; k + 1 < Qz; ++k) B[k] = BL[k + 1];
    B[Qz - 1] = 1e7;
    std::vector<double> cdfpts(Qz + 1, 0.0);
    for (std::size_t k = 0; k < Qz; ++k) cdfpts[k + 1] = B[k];
    const MultiRegimeResult res =
        mfq_multiregime(Qy, R, Qt, Rt, B, std::vector<double>(), cdfpts);
    const std::vector<std::vector<double>>& F = res.cdf;    // P(X < p)
    const std::vector<std::vector<double>>& Fm = res.cdfm;  // P(X <= p)
 
    // ---- the per-regime masses, weighted by the per-phase arrival rate -----
    std::vector<double> lmap(m, 0.0);
    for (std::size_t i = 0; i < m; ++i)
        for (std::size_t j = 0; j < m; ++j) lmap[i] += D(i, j);
    // E[min(v, P)] = alpha (-S)^-1 (I - expm(S v)) e for the patience P = (alpha, S = D0)
    const std::vector<double> alpha = map_pie(pat);
    const std::size_t np = P0.rows();
    Matrix<double> negS(np, np, 0.0);
    for (std::size_t i = 0; i < np; ++i)
        for (std::size_t j = 0; j < np; ++j) negS(i, j) = -P0(i, j);
    std::vector<double> meanMin(Qz, 0.0);
    for (std::size_t k = 0; k < Qz; ++k) {
        const Matrix<double> E = expm(P0, vrep[k]);
        std::vector<double> rhs(np, 0.0);
        for (std::size_t i = 0; i < np; ++i) {
            rhs[i] = 1.0;
            for (std::size_t j = 0; j < np; ++j) rhs[i] -= E(i, j);
        }
        const std::vector<double> x = solve(negS, rhs);
        for (std::size_t i = 0; i < np; ++i) meanMin[k] += alpha[i] * x[i];
    }
    const std::size_t norm = S * m;
    double total = 0.0, abandon = 0.0, queued = 0.0;
    for (std::size_t i = 0; i < n; ++i) total += Fm[0][i] * lmap[i % m];
    for (std::size_t d = 0; d < Qz; ++d)
        for (std::size_t i = 0; i < norm; ++i) {
            const double mass = (F[d + 1][i] - Fm[d][i]) * lmap[i % m];
            total += mass;
            abandon += ga[d] * mass;
            queued += meanMin[d] * mass;
        }
    const double abandonProb = abandon / total;
    const double expectedQueued = queued / total;
 
    // pi * D1 * e, the stationary arrival rate, not the per-phase rates added up
    Map<double> arv;
    arv.D0 = C;
    arv.D1 = D;
    const double lambda = map_lambda(arv);
    // Abandoning customers occupy the waiting room too, so QLen = busy servers +
    // lambda * E[min(V, P)] (Little on the waiting room) and RespT = QLen / Tput.
    RenegingMetrics out;
    out.tput = lambda * (1.0 - abandonProb);
    out.util = out.tput / (mu * static_cast<double>(S));
    out.qlen = out.tput / mu + lambda * expectedQueued;
    out.respt = out.qlen / out.tput;
    return out;
}
 
}  // namespace reneging_detail
 
/**
 * True when some station-class pair declares RENEGING, with or without a
 * patience law this analyzer can read. A RENEGING pair whose law is not
 * Markovian must still be seen by the gates: read off `mam_has_reneging_patience`
 * alone, such a model would fall through to dec.source and be solved with the
 * patience ignored.
 */
template <class T>
bool mam_declares_reneging(const qn::NetworkStruct<T>& L) {
    for (std::size_t i = 1; i <= L.nstations; ++i)
        for (std::size_t r = 1; r <= L.nclasses; ++r)
            if (reneging_detail::declares(L, i, r)) return true;
    return false;
}
 
/** `mam_has_reneging_patience.m`: RENEGING together with a usable patience law. */
template <class T>
bool mam_has_reneging_patience(const qn::NetworkStruct<T>& L) {
    for (std::size_t i = 1; i <= L.nstations; ++i)
        for (std::size_t r = 1; r <= L.nclasses; ++r)
            if (reneging_detail::markovian_patience(L, i, r) != nullptr) return true;
    return false;
}
 
/**
 * Port of `mam_reneging_applicable.m`: is this the MAP/M/s+G shape? Returns
 * rather than throws, since the dispatch uses it as a branch predicate; the
 * messages are the reference's.
 */
template <class T>
MamRenegingInfo mam_reneging_detect(const qn::NetworkStruct<T>& L) {
    MamRenegingInfo info;
    for (std::size_t k = 0; k < L.nclasses; ++k)
        if (std::isfinite(L.classes[k].population)) {
            info.why = "MAPMsG requires open queueing model.";
            return info;
        }
    if (L.nclasses > 1) {
        info.why = "MAPMsG currently supports single class only.";
        return info;
    }
    info.cls = 1;
    for (std::size_t i = 1; i <= L.nstations; ++i) {
        const qn::NodeType nt = L.stations[i - 1].nodetype;
        if (nt == qn::NodeType::Source) {
            info.source = i;
        } else if (nt == qn::NodeType::Queue) {
            if (info.station == 0) {
                info.station = i;
            } else {
                info.station = 0;
                info.why = "MAPMsG requires single queue station.";
                return info;
            }
        }
    }
    if (info.source == 0) {
        info.why = "No Source node found.";
        return info;
    }
    if (info.station == 0) {
        info.why = "No Queue node found.";
        return info;
    }
    const std::size_t q = info.station, src = info.source, r = info.cls;
    const qn::Station<T>& st = L.stations[q - 1];
    if (st.impatience.empty()) {
        info.why = "No patience/impatience configuration found.";
        return info;
    }
    if (!reneging_detail::declares(L, q, r)) {
        info.why = "Queue does not have reneging configured.";
        return info;
    }
    if (reneging_detail::markovian_patience(L, q, r) == nullptr) {
        // a declared law that is not Markovian is named (MATLAB mam_reneging_applicable, JAR renegingApplicable)
        const bool named = r - 1 < st.patience.size() && !st.patience[r - 1].disabled &&
                           st.patience[r - 1].type != lang::ProcessType::DISABLED;
        info.why = named ? std::string("The patience of this class is ") +
                               lang::process_to_text(st.patience[r - 1].type) +
                               ", not a Markovian law (the patience must be a Markovian law: "
                               "Exp, Erlang, HyperExp, Coxian, APH or PH)."
                         : std::string("No patience distribution for this class (the patience must be a "
                                       "Markovian law: Exp, Erlang, HyperExp, Coxian, APH or PH).");
        return info;
    }
    if (st.sched != lang::SchedStrategy::FCFS) {
        info.why = "MAPMsG requires FCFS scheduling.";
        return info;
    }
    const lang::Distrib<T>& svc = L.service[q - 1][r - 1];
    if (L.disabled[q - 1][r - 1] || svc.disabled) {
        info.why = "Invalid service process.";
        return info;
    }
    if (!svc.has_map() || svc.D0.rows() != 1) {
        info.why = "MAPMsG requires exponential service (single-phase).";
        return info;
    }
    const lang::Distrib<T>& arv = L.service[src - 1][r - 1];
    if (L.disabled[src - 1][r - 1] || arv.disabled || !arv.has_map() ||
        !basic_detail::is_markovian_type(arv.type)) {
        info.why = "Invalid arrival process.";
        return info;
    }
    const double ns = L.stations[q - 1].nservers;
    if (!std::isfinite(ns) || ns < 1.0) {
        info.why = "MAPMsG requires a finite, positive server count.";
        return info;
    }
    info.nservers = static_cast<int>(ns);
    info.ok = true;
    return info;
}
 
/**
 * The sentence `supportsModelMethod('default')` refuses a declared reneging
 * with when the model is not the MAP/M/s+G shape, or combines it with
 * setup/delay-off; empty otherwise. dec.source would serve the patience as if
 * it were absent, so a 'default' run must not fall through to it.
 */
template <class T>
std::string mam_reneging_default_refusal(const qn::NetworkStruct<T>& L) {
    if (!mam_declares_reneging(L)) return std::string();
    const MamRenegingInfo info = mam_reneging_detect(L);
    if (!info.ok)
        return "SolverMAM: this model declares a reneging patience outside the MAP/M/s+G shape "
               "(one Source, one FCFS Queue with exponential service, one open class) that the "
               "default method solves through the reneging analyzer; dec.source would ignore the "
               "patience. " + info.why + " Use SolverCTMC, SolverSSA, SolverJMT or SolverLDES.";
    if (!L.setupparam.empty())
        return "SolverMAM: this model combines setup/delay-off times with a reneging topology; only "
               "the dec.source decomposition honours setup times, and the default method routes "
               "such a topology to an analyzer that does not read them. Use SolverJMT or SolverLDES.";
    return std::string();
}
 
/**
 * The reneging route of `solver_mam_retrial.m`.
 *
 * @param L    the refreshed struct
 * @param opt  the MAM options; `mapmsg_quantization` is the regime count
 * @param info a successful `mam_reneging_detect(L)`
 */
template <class T>
mva::MvaSolution<T> solver_mam_reneging(const qn::NetworkStruct<T>& L, const MamOptions& opt,
                                        const MamRenegingInfo& info) {
    if (!info.ok)
        throw UnsupportedError(ErrorKind{"Reneging"},
                               "SolverMAM: not a MAP/M/s+G reneging queue: " + info.why);
    if (opt.mapmsg_quantization < 2)
        throw InputError("SolverMAM: options.config.mapmsg_quantization must be at least 2, got " +
                         std::to_string(opt.mapmsg_quantization) +
                         "; one regime leaves no grid between 0 and the patience horizon");
    if constexpr (!std::is_same<T, double>::value) {
        throw UnsupportedError(ErrorKind{"Arithmetic"},
            "SolverMAM: the MAP/M/s+G reneging analyzer solves a multi-regime fluid queue "
            "through mfq_multiregime, whose ordered real Schur basis is double only; use --arith "
            "double");
    } else {
        const std::size_t q = info.station, src = info.source, r = info.cls;
        const lang::Distrib<T>& arv = L.service[src - 1][r - 1];
        const lang::Distrib<T>& svc = L.service[q - 1][r - 1];
        const lang::Distrib<T>* pat = reneging_detail::markovian_patience(L, q, r);
        const double mu = -svc.D0(0, 0);
        const reneging_detail::RenegingMetrics met =
            reneging_detail::solve(arv.D0, arv.D1, mu, info.nservers, pat->D0, pat->D1,
                                   opt.mapmsg_quantization);
 
        mva::MvaSolution<T> out;
        const std::size_t M = L.nstations, K = L.nclasses;
        out.Q = Matrix<T>(M, K, 0.0);
        out.U = Matrix<T>(M, K, 0.0);
        out.R = Matrix<T>(M, K, 0.0);
        out.Tp = Matrix<T>(M, K, 0.0);
        out.C.assign(K, 0.0);
        out.X.assign(K, 0.0);
        out.Q(q - 1, r - 1) = met.qlen;
        out.U(q - 1, r - 1) = met.util;
        out.R(q - 1, r - 1) = met.respt;
        out.Tp(q - 1, r - 1) = met.tput;
        out.X[r - 1] = met.tput;
        out.C[r - 1] = met.respt;
        // the reference's totiter for this route is the regime count
        out.iter = static_cast<int>(opt.mapmsg_quantization);
        return out;
    }
}
 
}  // namespace mam
}  // namespace line
 
#endif  // LINE_SOLVERS_MAM_SOLVER_MAM_RENEGING_H