sn_patience_handles.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_SN_SN_PATIENCE_HANDLES_H
#define LINE_API_SN_SN_PATIENCE_HANDLES_H
 
/**
 * @file
 * @ingroup api_sn
 * Port of matlab/src/api/sn/sn_patience_handles.m.
 *
 * The abandonment solvers need the patience law as FUNCTIONS -- a complementary
 * cdf and a hazard rate -- not as moments, because that is what the underlying
 * theory consumes: Whitt's engineering solution reads the hazard near the
 * origin, and the fluid models integrate the ccdf. This port reads the law from
 * `Station::patience`, which holds the distribution itself rather than the
 * (D0,D1) pair the MATLAB and Java structs carry.
 *
 * ARITHMETIC: transcendental. The ccdf of a phase-type law is a matrix
 * exponential.
 */
 
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
 
#include "line/api/mam/map_cdf.h"
#include "line/api/mam/map_pdf.h"
#include "line/api/qsys/qsys_mgisrgi_whitt.h"
#include "line/api/qsys/qsys_types.h"
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
#include "line/num/number.h"
 
namespace line {
namespace api {
 
/** The patience law of one station-class pair, in the forms the solvers consume. */
template <class T>
struct PatienceHandles {
    /** Whether the pair declares reneging at all; everything below is unset when false. */
    bool present = false;
    /** F^c(t) = P(patience > t). */
    std::function<T(const T&)> ccdf;
    /** The patience density. */
    std::function<T(const T&)> pdf;
    /** The hazard rate h = f/(1-F). */
    std::function<T(const T&)> hazard;
    /** Mean patience. */
    T mean = num_traits<T>::from_int(0);
    /** Whether the law is exponential, in which case the analysis is exact. */
    bool isExponential = false;
    /** The abandonment rate 1/mean. */
    T rate = num_traits<T>::from_int(0);
 
    /** This law as the `Patience` argument of the abandonment solvers. */
    qsys::Patience<T> as_patience() const {
        if (isExponential) return qsys::Patience<T>::exponential(rate);
        return qsys::Patience<T>::hazard(hazard);
    }
};
 
/**
 * Is `d` a Markovian law, i.e. one with a (D0,D1) pair the abandonment solvers
 * can read exactly? The same family `basic_detail::is_markovian_type` names
 * for SolverMAM's patience.
 */
template <class T>
bool patience_is_markovian(const lang::Distrib<T>& d) {
    switch (d.type) {
        case lang::ProcessType::EXP:
        case lang::ProcessType::ERLANG:
        case lang::ProcessType::HYPEREXP:
        case lang::ProcessType::PH:
        case lang::ProcessType::APH:
        case lang::ProcessType::ME:
        case lang::ProcessType::RAP:
        case lang::ProcessType::MAP:
        case lang::ProcessType::MMAP:
        case lang::ProcessType::MPH:
        case lang::ProcessType::COXIAN:
        case lang::ProcessType::COX2:
        case lang::ProcessType::MMPP2:
            return true;
        default:
            return false;
    }
}
 
/**
 * Port of `sn_patience_refusal.m`: the by-name refusal of a reneging patience
 * that is not a Markovian law at station `ist` (0-based), class `r`, or "" when
 * there is none. A Det, Gamma, Lognormal, Pareto, Uniform or Weibull patience is
 * REFUSED, never replaced by an exponential of the same mean (or a fitted PH);
 * SolverMAM, SolverMVA and SolverFLD raise this sentence in the gate and the run.
 */
template <class T>
std::string sn_patience_refusal(const qn::NetworkStruct<T>& sn, std::size_t ist, std::size_t r) {
    if (ist >= sn.nstations || r >= sn.nclasses) return std::string();
    const qn::Station<T>& st = sn.stations[ist];
    if (r >= st.impatience.size() || st.impatience[r] != lang::ImpatienceType::RENEGING) return std::string();
    if (r >= st.patience.size() || st.patience[r].disabled ||
        st.patience[r].type == lang::ProcessType::DISABLED)
        return std::string();  // reneging declared with no law at all: nothing to name
    if (patience_is_markovian(st.patience[r])) return std::string();
    return std::string("The patience of this class is ") + lang::process_to_text(st.patience[r].type) +
           ", not a Markovian law (the patience must be a Markovian law: Exp, Erlang, HyperExp, "
           "Coxian, APH or PH).";
}
 
/** `sn_patience_refusal` over every station-class pair: the first refusal, or "". */
template <class T>
std::string sn_patience_refusal(const qn::NetworkStruct<T>& sn) {
    for (std::size_t i = 0; i < sn.nstations; ++i)
        for (std::size_t r = 0; r < sn.nclasses; ++r) {
            std::string why = sn_patience_refusal(sn, i, r);
            if (!why.empty()) return why;
        }
    return std::string();
}
 
/**
 * Port of `sn_patience_exp_refusal.m`: the by-name refusal of a reneging
 * patience that is not EXPONENTIAL at station `ist` (0-based), class `r`, or ""
 * when there is none. The stricter twin of `sn_patience_refusal`, worded in the
 * same shape, for the analyses that read the patience as one memoryless rate:
 * SolverCTMC and SolverSSA, whose RENEGE event keeps no patience phase per
 * waiting job, and SolverMVA method 'erlanga', the M/M/s/r+M chain. Callers
 * append their own redirect.
 */
template <class T>
std::string sn_patience_exp_refusal(const qn::NetworkStruct<T>& sn, std::size_t ist, std::size_t r) {
    if (ist >= sn.stations.size()) return std::string();
    const qn::Station<T>& st = sn.stations[ist];
    if (r >= st.impatience.size() || st.impatience[r] != lang::ImpatienceType::RENEGING) return std::string();
    if (r >= st.patience.size() || st.patience[r].disabled || st.patience[r].type == lang::ProcessType::EXP)
        return std::string();
    return std::string("The patience of this class is ") + lang::process_to_text(st.patience[r].type) +
           ", not an exponential law (the patience must be an exponential law: Exp).";
}
 
/** `sn_patience_exp_refusal` over every station-class pair: the first refusal, or "". */
template <class T>
std::string sn_patience_exp_refusal(const qn::NetworkStruct<T>& sn) {
    for (std::size_t i = 0; i < sn.stations.size(); ++i)
        for (std::size_t r = 0; r < sn.stations[i].impatience.size(); ++r) {
            std::string why = sn_patience_exp_refusal(sn, i, r);
            if (!why.empty()) return why;
        }
    return std::string();
}
 
/**
 * Build the patience handles of station `ist` (0-based), class `r`.
 *
 * Returns a handle set with `present == false` when the pair declares no
 * Markovian reneging patience; a non-Markovian one is named by
 * `sn_patience_refusal` (before 2026-10-22 it was read through dist_to_map).
 */
template <class T>
PatienceHandles<T> sn_patience_handles(const qn::NetworkStruct<T>& sn, std::size_t ist,
                                       std::size_t r) {
    PatienceHandles<T> h;
    if (ist >= sn.nstations || r >= sn.nclasses) return h;
    const qn::Station<T>& st = sn.stations[ist];
    if (r >= st.impatience.size() || st.impatience[r] != lang::ImpatienceType::RENEGING) return h;
    if (r >= st.patience.size() || st.patience[r].disabled) return h;
    if (!patience_is_markovian(st.patience[r])) return h;
 
    const lang::Distrib<T>& d = st.patience[r];
    const T zero = num_traits<T>::from_int(0);
    const T one = num_traits<T>::from_int(1);
    (void)one;
    h.present = true;
    h.mean = d.mean;
    h.rate = d.mean > zero ? T(one / d.mean) : zero;
    h.isExponential = (d.type == lang::ProcessType::EXP);
 
    // WHAT IS AVAILABLE UNDER EXACT ARITHMETIC, and what is not. Whether the
    // pair reneges, at what rate and under which family are all read off the
    // struct, so `present`, `mean`, `rate` and `isExponential` are answered for
    // every T -- which is what lets the method list and the method resolution
    // run under Rational. The ccdf and the hazard are a matrix exponential, so
    // they exist only where the arithmetic has one; the abandonment methods
    // themselves are refused by name before they are ever called.
    if constexpr (num_traits<T>::has_transcendental) {
        if (h.isExponential) {
            const T theta = h.rate;
            h.ccdf = [theta](const T& t) { return qsys::detail::num_exp(T(-theta * t)); };
            h.pdf = [theta](const T& t) { return T(theta * qsys::detail::num_exp(T(-theta * t))); };
            h.hazard = [theta](const T&) { return theta; };
            return h;
        }
 
        const mam::Map<T> m = lang::dist_to_map(d);
        h.ccdf = [m, one](const T& t) {
            std::vector<T> pts(1, t);
            return T(one - mam::map_cdf(m, pts)[0]);
        };
        h.pdf = [m](const T& t) {
            std::vector<T> pts(1, t);
            return mam::map_pdf(m, pts)[0];
        };
        const T asymptotic = h.rate;
        std::function<T(const T&)> ccdf = h.ccdf;
        std::function<T(const T&)> pdf = h.pdf;
        h.hazard = [ccdf, pdf, asymptotic](const T& t) {
            // h = f/(1-F). Past the point where the ccdf underflows the hazard
            // is the asymptotic decay rate, and returning that is better
            // conditioned than dividing two zeros.
            const T c = ccdf(t);
            if (num_traits<T>::to_double(c) <= 1e-300) return asymptotic;
            return T(pdf(t) / c);
        };
    }
    return h;
}
 
}  // namespace api
}  // namespace line
 
#endif  // LINE_API_SN_SN_PATIENCE_HANDLES_H