fluid_diffusion.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_FLUID_FLUID_DIFFUSION_H
#define LINE_SOLVERS_FLUID_FLUID_DIFFUSION_H
 
/**
 * @file
 * @ingroup line_solvers
 * The `diffusion` method: a port of `solver_fluid_diffusion.m`.
 *
 * WHAT IT COMPUTES. The fluid limit follows the MEAN drift and says nothing
 * about fluctuation. The diffusion approximation simulates the diffusion limit
 * by Euler-Maruyama. Every routing arc a -> b carries the flux phi_a(x) P(a,b),
 * with phi = mu*x at an INF station and mu*x/max(1,n_i) at a single server
 * (n_i the station total), plus its own Brownian noise sqrt(phi_a P(a,b) dt) Z
 * moved from a to b, so each step conserves the chain populations exactly. A
 * transfer that would drive a coordinate below zero is scaled down until it
 * lands there, which keeps the population exact without re-projecting.
 *
 * QN is the time average over the SECOND HALF of the trajectory (the first is
 * warm-up). TN is X_c * V(i,r) per chain, with X_c the mean departure rate
 * E[phi]/V at the chain's bottleneck (largest V/mu, single servers before INF
 * stations), so the throughputs are flow-balanced by construction; UN is TN/mu
 * at a single server and QN at an INF station.
 *
 * THE HORIZON IS NOT iter_max. `timestep` is the step and a finite
 * `timespan_end` sets the horizon, ceil(tend/dt) steps; otherwise 100000 steps
 * are taken. iter_max is the closing-iteration count every FLD run carries a
 * default for, and sizing the SDE with it once capped the run at ~200 steps.
 * That, together with a drift that served the whole queue at a single server
 * and unit noise per coordinate, made the method several-fold wrong.
 *
 * PARITY IS STATISTICAL, NOT EXACT. The trajectory is driven by pseudorandom
 * normals, and MATLAB's `randn` and this port's Mersenne Twister produce
 * different streams from the same seed. Two runs therefore agree in
 * distribution and not digit for digit -- the same situation the SSA solvers
 * are in across codebases. Tests must assert on averages and invariants
 * (population conservation, ordering), never on a specific trajectory.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <random>
#include <vector>
 
#include "line/api/mc/dtmc_stochcomp.h"
#include "line/lang/lang_types.h"
#include "line/lang/qn/network_struct.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace fluid {
 
/** Controls of the Euler-Maruyama trajectory. */
struct DiffusionOptions {
    std::size_t steps = 100000;     ///< number of steps (from timespan_end / dt when finite)
    double dt = 0.01;               ///< step size (the reference's timestep)
    unsigned long seed = 23000;     ///< RNG seed; see the header on parity
};
 
/** Steady-state averages of the diffusion trajectory. */
struct DiffusionResult {
    Matrix<double> QN, UN, RN, TN;
};
 
/**
 * Run the diffusion approximation of `sn`.
 *
 * Refuses by name every model shape the reference refuses: open classes, a
 * Source, a discipline outside {PS, FCFS, INF, SIRO}, and any finite
 * multiserver station.
 */
template <class T>
DiffusionResult fluid_diffusion(const qn::NetworkStruct<T>& sn, const DiffusionOptions& opt) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
    for (std::size_t r = 0; r < K; ++r)
        if (!std::isfinite(sn.classes[r].population))
            throw UnsupportedError(ErrorKind{"OpenClass"},
                "fluid diffusion: the method supports closed networks only; class '" +
                sn.classes[r].name + "' is open");
    for (std::size_t i = 0; i < M; ++i) {
        const lang::SchedStrategy sc = sn.stations[i].sched;
        if (sc == lang::SchedStrategy::EXT)
            throw UnsupportedError(ErrorKind{"OpenClass"}, "fluid diffusion: a Source is not supported (station '" +
                                   sn.stations[i].name + "'); the method is for closed networks");
        if (!(sc == lang::SchedStrategy::PS || sc == lang::SchedStrategy::FCFS ||
              sc == lang::SchedStrategy::INF || sc == lang::SchedStrategy::SIRO))
            throw UnsupportedError(ErrorKind{"SchedStrategy"},
                "fluid diffusion: scheduling at station '" + sn.stations[i].name +
                "' is outside the supported set {PS, FCFS, INF, SIRO}");
        const double c = sn.stations[i].nservers;
        if (std::isfinite(c) && c > 1.0)
            throw UnsupportedError(ErrorKind{"MultiServer"},
                "fluid diffusion: only single-server or infinite-server stations are supported; "
                "station '" + sn.stations[i].name + "' has more than one server");
    }
 
    // Service rate per (station, class); 0 marks "not served here".
    Matrix<double> mu(M, K, 0.0);
    std::vector<bool> inf_server(M, false);
    for (std::size_t i = 0; i < M; ++i) {
        inf_server[i] = !std::isfinite(sn.stations[i].nservers);
        for (std::size_t r = 0; r < K; ++r) {
            if (sn.disabled[i][r] || sn.service[i][r].D0.rows() == 0) continue;
            const double rate = num_traits<T>::to_double(sn.rates(i, r));
            if (rate > 0.0 && std::isfinite(rate)) mu(i, r) = rate;
        }
    }
 
    // Station-space routing, with the non-station stateful nodes folded out.
    const std::size_t S = sn.nof_stateful();
    std::vector<std::size_t> keep;
    keep.reserve(M * K);
    for (std::size_t i = 0; i < M; ++i) {
        const std::size_t isf = sn.stateful_of_station(i + 1) - 1;
        for (std::size_t r = 0; r < K; ++r) keep.push_back(isf * K + r);
    }
    Matrix<double> rt_full(S * K, S * K, 0.0);
    if (sn.rt.rows() == S * K)
        for (std::size_t a = 0; a < S * K; ++a)
            for (std::size_t b = 0; b < S * K; ++b)
                rt_full(a, b) = num_traits<T>::to_double(sn.rt(a, b));
    const Matrix<double> P = mc::dtmc_stochcomp(rt_full, keep);
 
    // Arcs between coordinates a = i*K + r; a self-loop moves no job.
    const std::size_t MK = M * K;
    std::vector<std::size_t> asrc, adst;
    std::vector<double> aprob;
    for (std::size_t a = 0; a < MK; ++a)
        for (std::size_t b = 0; b < MK; ++b)
            if (a != b && P(a, b) > 0.0) {
                asrc.push_back(a);
                adst.push_back(b);
                aprob.push_back(P(a, b));
            }
    const std::size_t narc = asrc.size();
 
    auto phi = [&](const std::vector<double>& x, std::vector<double>& f) {
        for (std::size_t i = 0; i < M; ++i) {
            double n = 0.0;
            for (std::size_t r = 0; r < K; ++r) n += x[i * K + r];
            const double den = inf_server[i] ? 1.0 : std::max(1.0, n);
            for (std::size_t r = 0; r < K; ++r) f[i * K + r] = mu(i, r) * x[i * K + r] / den;
        }
    };
 
    const double dt = opt.dt > 0.0 ? opt.dt : 0.01;
    const std::size_t steps = std::max<std::size_t>(2, opt.steps);
    std::mt19937_64 rng(opt.seed);
    std::normal_distribution<double> gauss(0.0, 1.0);
 
    // Each class starts spread evenly over the stations that serve it.
    std::vector<double> x(MK, 0.0), f(MK, 0.0);
    for (std::size_t r = 0; r < K; ++r) {
        std::size_t nsv = 0;
        for (std::size_t i = 0; i < M; ++i)
            if (mu(i, r) > 0.0) ++nsv;
        for (std::size_t i = 0; i < M; ++i)
            if (mu(i, r) > 0.0) x[i * K + r] = sn.classes[r].population / static_cast<double>(nsv);
    }
    phi(x, f);
 
    const std::size_t burn = steps / 2;
    std::vector<double> qsum(MK, 0.0), fsum(MK, 0.0), gain(MK), loss(MK), xn(MK), theta(MK);
    std::vector<std::size_t> loser(narc), gainer(narc);
    std::vector<double> amt(narc);
    for (std::size_t step = 1; step < steps; ++step) {
        for (std::size_t e = 0; e < narc; ++e) {
            const double fl = f[asrc[e]] * aprob[e] * dt;
            const double inc = fl + std::sqrt(fl) * gauss(rng);
            // a transfer drains its source when positive and its destination when negative
            loser[e] = inc > 0.0 ? asrc[e] : adst[e];
            gainer[e] = inc > 0.0 ? adst[e] : asrc[e];
            amt[e] = std::fabs(inc);
        }
        for (int it = 0; it < 100; ++it) {
            std::fill(gain.begin(), gain.end(), 0.0);
            std::fill(loss.begin(), loss.end(), 0.0);
            for (std::size_t e = 0; e < narc; ++e) {
                gain[gainer[e]] += amt[e];
                loss[loser[e]] += amt[e];
            }
            bool neg = false;
            for (std::size_t a = 0; a < MK; ++a) {
                xn[a] = x[a] + gain[a] - loss[a];
                theta[a] = 1.0;
                if (xn[a] < -lang::GlobalConstants::FineTol) {
                    theta[a] = (x[a] + gain[a]) / loss[a];
                    neg = true;
                }
            }
            if (!neg) break;
            for (std::size_t e = 0; e < narc; ++e) amt[e] *= theta[loser[e]];
        }
        for (std::size_t a = 0; a < MK; ++a) x[a] = std::max(xn[a], 0.0);
        phi(x, f);
        if (step >= burn)
            for (std::size_t a = 0; a < MK; ++a) {
                qsum[a] += x[a];
                fsum[a] += f[a];
            }
    }
    const double navg = static_cast<double>(steps - std::max<std::size_t>(burn, 1));
 
    DiffusionResult out;
    out.QN = Matrix<double>(M, K, 0.0);
    out.UN = Matrix<double>(M, K, 0.0);
    out.RN = Matrix<double>(M, K, 0.0);
    out.TN = Matrix<double>(M, K, 0.0);
    for (std::size_t a = 0; a < MK; ++a) out.QN(a / K, a % K) = qsum[a] / navg;
 
    // Per chain: the mean rate at its bottleneck, spread over the visits.
    for (std::size_t c = 0; c < sn.nchains; ++c) {
        const Matrix<T>& V = sn.visits[c];
        auto vis = [&](std::size_t i, std::size_t r) {
            return num_traits<T>::to_double(V(sn.stateful_of_station(i + 1) - 1, r));
        };
        bool found = false, bfin = false;
        std::size_t bi = 0, br = 0;
        double bD = -1.0;
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                if (!sn.chains[c][r] || mu(i, r) <= 0.0) continue;
                const double v = vis(i, r);
                if (!(v > 0.0)) continue;
                // a single server outranks every INF station: its rate is bounded by mu, so the
                // boundary distorts it least, and only there does TN = X*V keep sum_r U <= 1
                const bool fin = !inf_server[i];
                const double D = v / mu(i, r);
                if (!found || (fin && !bfin) || (fin == bfin && D > bD)) {
                    found = true;
                    bi = i;
                    br = r;
                    bfin = fin;
                    bD = D;
                }
            }
        if (!found) continue;
        const double Xc = (fsum[bi * K + br] / navg) / vis(bi, br);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r)
                if (sn.chains[c][r] && mu(i, r) > 0.0) out.TN(i, r) = Xc * vis(i, r);
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t r = 0; r < K; ++r) {
            out.UN(i, r) = inf_server[i] ? out.QN(i, r) : (mu(i, r) > 0.0 ? out.TN(i, r) / mu(i, r) : 0.0);
            // TN is zero only to the integrator's accuracy.
            if (out.TN(i, r) > lang::GlobalConstants::Zero) out.RN(i, r) = out.QN(i, r) / out.TN(i, r);
        }
    return out;
}
 
}  // namespace fluid
}  // namespace line
 
#endif  // LINE_SOLVERS_FLUID_FLUID_DIFFUSION_H