qsys_res.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_RES_H
#define LINE_API_QSYS_QSYS_RES_H
 
/**
 * @file
 * @ingroup api_qsys
 * The `res` struct the exact queueing-system functions return, and the numeric
 * helpers those functions share.
 *
 * MATLAB's exact closed forms (qsys_mm1, qsys_mmk, qsys_mg1, qsys_gm1,
 * qsys_mmck, qsys_mg1_prio and the finite-source, priority, busy-period and
 * renewal-arrival functions) return a struct whose field set depends on the
 * options passed (tPoints, percentiles, ES3, lambda, ...). `QsysRes` holds the
 * same fields in the same order, by the same names, so the dispatch layer can
 * encode them one for one; `analyzer` is MATLAB's `res.analyzer` string.
 *
 * The distributional measures (incomplete gamma tails, Bessel densities,
 * log-space distributions, bisection percentiles) are evaluated in double
 * precision, as in the reference.
 */
 
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <functional>
#include <limits>
#include <string>
#include <utility>
#include <vector>
 
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>
 
#include "line/util/error.h"
 
namespace line {
namespace qsys {
 
/** One measure of a `res` struct: a scalar or a row vector. */
template <class T = double>
struct QsysResFieldT {
    std::string name;
    std::vector<T> value;
    bool scalar = true;
};
 
/**
 * MATLAB's `res` struct: named measures in insertion order, plus `res.analyzer`.
 * Templated on the arithmetic so that the closed forms which are pure field
 * arithmetic keep their exact Rational values; `QsysRes` is the double one.
 */
template <class T = double>
struct QsysResT {
    std::vector<QsysResFieldT<T>> fields;
    std::string analyzer;
 
    /** Sets (or overwrites) a scalar field. */
    void set(const std::string& name, const T& v) { put(name, std::vector<T>{v}, true); }
    /** Sets (or overwrites) a vector field. */
    void set(const std::string& name, const std::vector<T>& v) { put(name, v, false); }
    bool has(const std::string& name) const { return find(name) != nullptr; }
    /** The scalar field `name`; throws when absent. */
    T scalar(const std::string& name) const {
        const QsysResFieldT<T>* f = find(name);
        if (f == nullptr || f->value.empty()) throw Error("QsysRes: no field '" + name + "'");
        return f->value[0];
    }
    /** The vector field `name`; throws when absent. */
    const std::vector<T>& vec(const std::string& name) const {
        const QsysResFieldT<T>* f = find(name);
        if (f == nullptr) throw Error("QsysRes: no field '" + name + "'");
        return f->value;
    }
 
  private:
    const QsysResFieldT<T>* find(const std::string& name) const {
        for (const QsysResFieldT<T>& f : fields)
            if (f.name == name) return &f;
        return nullptr;
    }
    void put(const std::string& name, const std::vector<T>& v, bool sc) {
        for (QsysResFieldT<T>& f : fields)
            if (f.name == name) {
                f.value = v;
                f.scalar = sc;
                return;
            }
        fields.push_back(QsysResFieldT<T>{name, v, sc});
    }
};
 
using QsysRes = QsysResT<double>;
 
/**
 * The name-value options shared by the exact functions: `tPoints`,
 * `percentiles` and `nmax` (a negative nmax means "not given", MATLAB's []).
 */
struct QsysResOptions {
    std::vector<double> tPoints;
    std::vector<double> percentiles;
    long nmax = -1;
};
 
namespace resdetail {
 
/** MATLAB's gammaln. */
inline double gammaln(double x) { return std::lgamma(x); }
 
/** MATLAB's gammainc(x, a, 'upper'), the regularized upper tail Q(a, x). */
inline double gammainc_upper(double x, double a) {
    if (x <= 0.0) return 1.0;
    if (std::isinf(x)) return 0.0;
    return boost::math::gamma_q(a, x);
}
 
/** MATLAB's gammainc(x, a, 'lower'), the regularized lower tail P(a, x). */
inline double gammainc_lower(double x, double a) {
    if (x <= 0.0) return 0.0;
    if (std::isinf(x)) return 1.0;
    return boost::math::gamma_p(a, x);
}
 
/**
 * MATLAB's besseli(1, x, 1) = I_1(x) exp(-x), x >= 0. Past the range where
 * I_1 itself is finite the Hankel asymptotic expansion of the scaled function
 * is used (six terms, relative error far below 1e-15 at x > 600).
 */
inline double besseli1_scaled(double x) {
    if (x < 600.0) return boost::math::cyl_bessel_i(1.0, x) * std::exp(-x);
    const double mu4 = 4.0;  // 4 nu^2
    double term = 1.0, sum = 1.0;
    for (int k = 1; k <= 6; ++k) {
        const double odd = 2.0 * k - 1.0;
        term *= -(mu4 - odd * odd) / (k * 8.0 * x);
        sum += term;
    }
    return sum / std::sqrt(2.0 * M_PI * x);
}
 
/** The 0:nmax power law c*r.^(0:nmax). */
inline std::vector<double> geometric(double c, double r, long nmax) {
    std::vector<double> v(static_cast<std::size_t>(nmax + 1));
    for (long n = 0; n <= nmax; ++n) v[static_cast<std::size_t>(n)] = c * std::pow(r, static_cast<double>(n));
    return v;
}
 
/**
 * The bisection inverse of a complementary CDF shared by the reference's
 * functions (ccdf_inverse / mmck_ccdf_inverse / mmckk_ccdf_inverse): the time at
 * which `ccdf` falls to `target`, bracketing upward from `hi`.
 */
inline double ccdf_inverse(const std::function<double(double)>& ccdf, double target, double hi) {
    if (ccdf(0.0) <= target) return 0.0;
    double lo = 0.0;
    while (ccdf(hi) > target) hi = 2.0 * hi;
    for (int it = 1; it <= 200; ++it) {
        const double mid = (lo + hi) / 2.0;
        if (ccdf(mid) > target) lo = mid;
        else hi = mid;
        if (hi - lo <= 1e-12 * std::max(1.0, hi)) break;
    }
    return (lo + hi) / 2.0;
}
 
/** sum(x) */
inline double vsum(const std::vector<double>& x) {
    double s = 0.0;
    for (double v : x) s += v;
    return s;
}
 
/** exp(logp - logsumexp(lognorm)), the reference's log-space normalization. */
inline std::vector<double> lognormalize(const std::vector<double>& logp, const std::vector<double>& lognorm) {
    double mx = -std::numeric_limits<double>::infinity();
    for (double v : lognorm) mx = std::max(mx, v);
    double s = 0.0;
    for (double v : lognorm) s += std::exp(v - mx);
    const double ln = mx + std::log(s);
    std::vector<double> p(logp.size());
    for (std::size_t i = 0; i < logp.size(); ++i) p[i] = std::exp(logp[i] - ln);
    return p;
}
 
/** Fills `out` with 1 - ccdf(t) at each t, as `1 - arrayfun(ccdf, t)`. */
inline std::vector<double> cdf_of(const std::function<double(double)>& ccdf, const std::vector<double>& t) {
    std::vector<double> v(t.size());
    for (std::size_t i = 0; i < t.size(); ++i) v[i] = 1.0 - ccdf(t[i]);
    return v;
}
 
/** The percentiles by bisection, `arrayfun(@(r) ccdf_inverse(ccdf, 1 - r/100, hi), pct)`. */
inline std::vector<double> pct_of(const std::function<double(double)>& ccdf, const std::vector<double>& pct,
                                  double hi) {
    std::vector<double> v(pct.size());
    for (std::size_t i = 0; i < pct.size(); ++i) v[i] = ccdf_inverse(ccdf, 1.0 - pct[i] / 100.0, hi);
    return v;
}
 
/** The reference's `LINE:input:Unstable` refusal. */
inline void require_stable(double rho, const char* fn) {
    if (rho >= 1.0) {
        char buf[64];
        std::snprintf(buf, sizeof buf, "%g", rho);
        throw InputError(ErrorKind("Unstable"), std::string(fn) + ": unstable system: rho = " + buf + " >= 1");
    }
}
 
/** A `LINE:input:Unstable` refusal "fn: unstable system: <what> = %g >= 1". */
[[noreturn]] inline void throw_unstable(const char* fn, const char* what, double v) {
    char buf[64];
    std::snprintf(buf, sizeof buf, "%g", v);
    throw InputError(ErrorKind("Unstable"), std::string(fn) + ": unstable system: " + what + " = " + buf + " >= 1");
}
 
/** MATLAB's "unknown option 'name'" refusal, as LINE:input:Argument or as a bare LINE:Error. */
[[noreturn]] inline void throw_unknown_option(const char* fn, const std::string& name, bool argument_id) {
    const std::string msg = std::string(fn) + ": unknown option '" + name + "'";
    if (argument_id) throw InputError(ErrorKind("Argument"), msg);
    throw Error(msg);
}
 
/**
 * The reference's nchoosek_r: exp(gammaln) rounded to the integer it is while
 * that integer is exactly representable (below 2^52).
 */
inline double nchoosek_r(double n, double k) {
    const double v = std::exp(std::lgamma(n + 1.0) - std::lgamma(k + 1.0) - std::lgamma(n - k + 1.0));
    return v < 4503599627370496.0 ? std::round(v) : v;
}
 
/** A percentile vector checked against (0, 100), MATLAB's LINE:Error refusal. */
inline void check_percentiles(const std::vector<double>& pct, const char* fn) {
    for (double p : pct)
        if (p <= 0.0 || p >= 100.0) throw Error(std::string(fn) + ": percentiles must lie in (0, 100)");
}
 
}  // namespace resdetail
}  // namespace qsys
}  // namespace line
 
#endif  // LINE_API_QSYS_QSYS_RES_H