pfqn_lldfun.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_API_PFQN_LLDFUN_H
#define LINE_API_PFQN_LLDFUN_H
 
/**
 * @file
 * @ingroup api_pfqn
 * AMVA-QD limited-load-dependence function.
 *
 * Templated port of matlab/src/api/pfqn/pfqn_lldfun.m, cross-checked against
 * jar/src/main/java/jline/api/pfqn/ld/Pfqn_lldfun.java. Returns, for every
 * station i, the reciprocal service-capacity multiplier at the (fractional)
 * queue length n(i), combining
 *
 *  - the multiserver term 1/softmin(n_i, c_i, alpha), a smooth approximation
 *    of 1/min(n_i, c_i), with alpha = 20 as in the reference. A delay station
 *    (c_i infinite) contributes 1; the caller handles the 1/n_i itself.
 *  - the generic load-dependence term 1/alpha_i(n_i), obtained by CLAMPED
 *    LINEAR interpolation of the lattice lldscaling(i, 1..smax) at
 *    min(max(n_i,1), smax). Linear is exact for the piecewise-linear
 *    min(1:N, c) lattice that load dependence is overwhelmingly used to
 *    express, whereas a cubic spline overshoots between knots and extrapolates
 *    to negative rates past the lattice. A station whose lattice row is
 *    constant is skipped, exactly as `range(lldscaling(i,:)) > 0` does.
 *
 * softmin. The reference rewrites the literal weighted average
 * (x e^{-ax} + y e^{-ay})/(e^{-ax} + e^{-ay}) as lo + gap w/(1+w) with
 * w = e^{-a gap}, so the exponent argument is never positive and the limit
 * w -> 0 returns min(x,y) exactly; that form is reproduced here, including its
 * cutoff at a*gap > 745 where exp underflows to zero in double.
 *
 * Arithmetic: TRANSCENDENTAL-GATED, and here on a genuine transcendental, not
 * on a tolerance: softmin evaluates exp(-alpha |n_i - c_i|), which is not an
 * element of the field generated by the inputs for any alpha and gap of
 * interest. There is no exact formulation of a soft minimum.
 */
 
#include <cmath>
#include <cstddef>
#include <limits>
#include <vector>
 
#include "line/num/number.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace pfqn {
 
namespace detail {
 
/** matlab/src/util/softmin.m, in the overflow-free form. */
template <class T>
T softmin(const T& x, const T& y, double alpha) {
    const T lo = x < y ? x : y;
    const T hi = x < y ? y : x;
    const T gap = hi - lo;
    // exp(-t) underflows to exactly 0 for t > 745.13; beyond that softmin is min.
    if (!(num_traits<T>::to_double(gap) < 745.0 / alpha)) return lo;
    using std::exp;
    const T w = exp(T(-num_traits<T>::from_double(alpha) * gap));
    return lo + gap * w / (num_traits<T>::from_int(1) + w);
}
 
}  // namespace detail
 
/**
 * @brief AMVA-QD limited-load-dependence function.
 *
 * @param n          (M) queue lengths, possibly fractional
 * @param lldscaling (M x smax) rate lattice; empty for none
 * @param nservers   (M) server counts; a non-finite entry marks a delay
 *                   station. Empty to skip the multiserver term entirely,
 *                   which is what omitting the argument does in MATLAB.
 * @return (M) reciprocals of the effective rate multipliers
 */
template <class T>
std::vector<T> pfqn_lldfun(const std::vector<T>& n, const Matrix<T>& lldscaling,
                           const std::vector<double>& nservers) {
    static_assert(num_traits<T>::has_transcendental,
                  "pfqn_lldfun requires transcendental arithmetic");
 
    const std::size_t M = n.size();
    const T one = num_traits<T>::from_int(1);
    std::vector<T> r(M, one);
    if (!nservers.empty() && nservers.size() != M)
        throw InputError("pfqn_lldfun: server-count vector has the wrong station count");
    if (!lldscaling.empty() && lldscaling.rows() != M)
        throw InputError("pfqn_lldfun: rate lattice has the wrong station count");
    const std::size_t smax = lldscaling.empty() ? 0 : lldscaling.cols();
    const double alpha = 20.0;
 
    for (std::size_t i = 0; i < M; ++i) {
        if (!nservers.empty()) {
            if (!std::isfinite(nservers[i])) {
                // Delay station: the 1/n_i is applied by the caller.
                r[i] = one;
            } else {
                const T c = num_traits<T>::from_double(nservers[i]);
                const T sm = detail::softmin(n[i], c, alpha);
                if (sm == num_traits<T>::from_int(0)) {
                    // The soft minimum collapses only when both arguments do.
                    const T m = n[i] < c ? n[i] : c;
                    if (m == num_traits<T>::from_int(0))
                        throw NumericError("pfqn_lldfun: zero queue length and zero server count");
                    r[i] = r[i] / m;
                } else {
                    r[i] = r[i] / sm;
                }
            }
        }
        if (smax > 0) {
            T lo = lldscaling(i, 0), hi = lldscaling(i, 0);
            for (std::size_t k = 1; k < smax; ++k) {
                if (lldscaling(i, k) < lo) lo = lldscaling(i, k);
                if (lldscaling(i, k) > hi) hi = lldscaling(i, k);
            }
            // A CONSTANT ROW IS A NO-OP ONLY WHEN IT IS ONE. The test here used to be
            // `hi > lo`, which skipped alpha(n) = c for EVERY c: a station declaring a
            // uniform rate multiplier was then solved at its UNSCALED rate, while the
            // exact recursions (solver_mvald, pfqn_mvaldmx) and SolverCTMC applied it, so
            // AMVA contradicted them on the same model -- silently, and by the whole
            // factor c. Only a row of ones divides by 1.
            if (hi != one || lo != one) {
                if (smax == 1) {
                    // One column IS the constant alpha(n): nothing to interpolate, and the
                    // two-knot arithmetic below has no second knot to read.
                    if (lldscaling(i, 0) == num_traits<T>::from_int(0))
                        throw NumericError("pfqn_lldfun: the rate multiplier is zero");
                    r[i] = r[i] / lldscaling(i, 0);
                    continue;
                }
                // Clamp into [1, smax] and interpolate linearly between knots.
                T x = n[i];
                const T xlo = one, xhi = num_traits<T>::from_int(static_cast<long>(smax));
                if (x < xlo) x = xlo;
                if (x > xhi) x = xhi;
                const double xd = num_traits<T>::to_double(x);
                std::size_t k0 = static_cast<std::size_t>(xd) - 1;
                if (k0 + 1 >= smax) k0 = smax - 2;
                const T t = x - num_traits<T>::from_int(static_cast<long>(k0) + 1);
                const T val =
                    lldscaling(i, k0) + t * (lldscaling(i, k0 + 1) - lldscaling(i, k0));
                if (val == num_traits<T>::from_int(0))
                    throw NumericError("pfqn_lldfun: the interpolated rate multiplier is zero");
                r[i] = r[i] / val;
            }
        }
    }
    return r;
}
 
/** Overload without the multiserver term, matching the two-argument MATLAB call. */
template <class T>
std::vector<T> pfqn_lldfun(const std::vector<T>& n, const Matrix<T>& lldscaling) {
    return pfqn_lldfun(n, lldscaling, std::vector<double>());
}
 
}  // namespace pfqn
}  // namespace line
 
#endif  // LINE_API_PFQN_LLDFUN_H