solver_env_limit.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_ENV_SOLVER_ENV_LIMIT_H
#define LINE_SOLVERS_ENV_SOLVER_ENV_LIMIT_H
 
/**
 * @file
 * @ingroup line_solvers
 * The two CLOSED-FORM environment limits: `SolverENV.solveEnvLimit`, reached by
 * `options.method` in {`avg`, `dec`}.
 *
 * Both replace the fixed point with a single reading, and each is exact at one
 * end of the time-scale separation between the environment and the network:
 *
 *   `avg`, the FAST-environment limit. The environment switches so much faster
 *   than the network responds that the network only ever sees the AVERAGE of the
 *   modulated rates. One rate-averaged model is built -- every station-class
 *   rate that varies across stages replaced by its probEnv-weighted mean, as an
 *   exponential -- and solved once, in steady state. Exact as the stage-switch
 *   rate goes to infinity.
 *
 *   `dec`, the SLOW-environment limit, i.e. quasi-stationary decomposition. The
 *   environment stays in a stage long enough for the network to reach that
 *   stage's own steady state, so each stage is solved independently and the
 *   metrics are averaged with weights probEnv. Exact as the stage-switch rate
 *   goes to zero.
 *
 * NEITHER CARRIES ANYTHING ACROSS A SWITCH, which is what makes them closed
 * form and also what they give up: no entry state, no reset policy, no
 * transient. A reset policy declared on an arc is therefore inert here, exactly
 * as it is in the reference, whose `solveEnvLimit` never reads `resetFun`.
 *
 * WHAT IS AVERAGED, AND WHAT IS LEFT ALONE (`buildRateAveragedModel`). Only a
 * Source, a Queue or a Delay carries a rate to average. A station-class pair
 * that is disabled in any stage, or whose rate does not actually vary across
 * them, keeps its ORIGINAL distribution rather than being rewritten as an
 * exponential of its own mean -- so a non-modulated Erlang stays an Erlang, and
 * the base model is preserved exactly outside the modulated rates.
 *
 * A STAGE HOLDING A CACHE is solved, not refused. The hit, miss and delayed
 * fractions each stage reports are blended the way the reference's
 * `accumCacheMetric` does after 2026-09-13: as RATES weighted by that stage's
 * per-class arrival INTO THE CACHE, normalised once at the end. Rates and not
 * ratios, because a stage's hit ratio is conditional on arriving in that stage,
 * so a probEnv-weighted mean of ratios disagrees with the Sink-row hit
 * throughput -- `sum_e p_e*lambda_e*h_e` -- in the SAME node table. It is also
 * the convention `aggregate_cache_meanfield` and the statevec `cache_blend`
 * already hold, each with a comment saying so.
 *
 * WHICH SOLVER RUNS A STAGE. The reference hands `SolverENV` a FACTORY and calls
 * it for every stage, so "the stage solver is the calling solver" is literally
 * true there. The second constructor here takes an `EnvStageAvgFn`, the same
 * thing as a callable, and the no-factory path runs `stage_solver`: the fluid
 * analyzer, with `ctmc` the enumerated chain at `stage_cutoff` (an ensemble built
 * on SolverCTMC stages), with `mva` / `nc` the MVA or NC runner at `stage_mva` /
 * `stage_nc` (an ensemble built on SolverMVA / SolverNC stages). The callable is
 * what `map_env_approx` uses, so that the CALLING solver's own options reach
 * every stage.
 *
 * `occupation`, `SolverENV.solveEnvOccupation`, is the `dec` twin for an
 * environment that need not be ergodic. Each stage is solved in steady state
 * as `dec` solves it, but weighted by the time the environment is expected to
 * spend there from the initial stage law `EnvOptions::env_init` (default: stage
 * 1), so the environment may be ABSORBING. Over a finite horizon [0,T] the
 * occupation is int_0^T pi0 exp(A t) dt and the weights occupation/T, the
 * time-averaged stage law, which needs exponential arcs to integrate. With
 * T = Inf and absorbing stages the occupation is pi0_t (-A_tt)^-1 over the
 * transient stages and the weights are the absorption probabilities; with no
 * absorbing stage the weights are probEnv, as `dec`. The T = Inf cases read the
 * embedded chain and the mean holding times only, so they hold for a
 * semi-Markov environment too. A stage of weight zero is not solved.
 *
 * `occupation` is also the one limit that takes LAYERED stages, as the
 * reference does: each is run by its own `SolverLN` (`EnvOptions::lqn`) to the
 * fixed point and read over its LQN ELEMENTS (`SolverLN.getAvg`, an nidx x 1
 * table), not over its layers: a zero-capacity stage shares the elements with
 * the others while its layers do not (ln_zero_capacity drops the layers only a
 * trapped reference task reaches). The stages must then be all layered with the
 * same elements, which `Environment::set_lqn_stage` enforces.
 * The reference solves stages sharing one model object once; a stage here holds
 * its model by value, so there is no identity to share and each is solved.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <string>
#include <type_traits>
#include <vector>
 
#include "line/api/sn/sn_node_metrics.h"
#include "line/lang/qn/environment.h"
#include "line/lang/qn/network_struct.h"
#include "line/num/number.h"
#include "line/solvers/cache_metrics.h"
#include "line/solvers/ctmc/solver_ctmc_analyzer.h"
#include "line/solvers/env/solver_env.h"
#include "line/solvers/fluid/fluid_runner.h"
#include "line/solvers/fluid/solver_fluid.h"
#include "line/solvers/mva/solver_mva_runner.h"
#include "line/solvers/nc/solver_nc_runner.h"
#include "line/util/error.h"
#include "line/util/expm.h"
#include "line/util/linalg.h"
#include "line/util/matrix.h"
 
namespace line {
namespace env {
 
/**
 * What one stage solve reports back to the limits.
 *
 * The literal translation of what the reference's `[Qe,Ue,~,Te] = se.getAvg()`
 * plus `se.getAvgNode()` hand back, minus the node table: a Cache reports on the
 * NODE, so `cache` carries what `getAvgNode`'s resync would have written there.
 */
template <class T>
struct EnvStageAvg {
    Matrix<double> QN, UN, TN;
    solvers::CacheMetrics<T> cache;
};
 
/**
 * The stage solver, as a callable: the C++ spelling of the MATLAB function
 * handle `SolverENV(renv, @(m) SolverX(m, opts))` passes.
 *
 * It is how a caller binds ITS OWN options into every stage (`map_env_approx`);
 * the named stage solvers of `EnvOptions::stage_solver` cover the plain
 * ensemble. There is no include cycle in the way of either: the MVA and NC
 * runners reach no `env/` header (`ln/solver_ln.h` already includes both).
 */
template <class T>
using EnvStageAvgFn = std::function<EnvStageAvg<T>(const qn::NetworkStruct<T>&)>;
 
/** What a limit solve reports. */
struct EnvLimitSolution {
    /** Environment-averaged metrics, (nstations x nclasses). */
    Matrix<double> QN, UN, TN;
    /** The limit that ran: `avg`, `dec` or `occupation`. */
    std::string method;
    /**
     * `dec` and `occupation`: each stage's own steady state, before the blend.
     * A stage `occupation` does not solve (weight zero) stays all zeros.
     */
    std::vector<Matrix<double> > QStage, UStage, TStage;
    /**
     * `occupation` only (`getStageOccupation`): the expected time spent in each
     * stage over the horizon, Inf for a stage the environment ends in when the
     * horizon is unbounded, and the weight of each stage in the blend.
     */
    std::vector<double> occupation, occupation_weight;
    /** The stage probabilities the blend used: probEnv, or the occupation weights. */
    std::vector<double> prob_env;
    /**
     * The environment-blended cache surface, one entry per Cache node.
     *
     * EMPTY on a model with no cache. Within an entry every field is optional
     * and ABSENT MEANS NOT COMPUTED, never zero, exactly as `CacheNodeMetrics`
     * states: a class no stage measured stays NaN rather than reading as "never
     * hits". Kept `double`, since `init()` refuses `T != double` anyway.
     */
    solvers::CacheMetrics<double> cache;
};
 
/**
 * The limit solver.
 *
 * `envObj` must already carry a model per stage; `init()` is called here, as
 * `solveEnvLimit` calls `self.init()` before reading probEnv.
 */
template <class T>
class SolverEnvLimit {
public:
    SolverEnvLimit(Environment<T>& e, const EnvOptions& o) : envObj(e), opt(o) { init(); }
 
    /**
     * With an injected stage solver: `stage_fn` runs every stage, exactly as the
     * reference calls its factory. `opt.stage_solver` is then not consulted,
     * since the caller has already named the solver by handing one over.
     */
    SolverEnvLimit(Environment<T>& e, const EnvOptions& o, EnvStageAvgFn<T> stage_fn)
        : envObj(e), opt(o), stage_fn_(stage_fn) {
        init();
    }
 
    EnvLimitSolution solve() {
        if (occ_) return solve_occupation();
        return dec_ ? solve_dec() : solve_avg();
    }
 
    /**
     * `stageOccupation_`: the expected time spent in each stage over the
     * horizon and the weight of each stage, read off the initialised
     * environment without solving any stage. See the file header.
     */
    static void stage_occupation(const Environment<T>& env, const EnvOptions& o,
                                 std::vector<double>& occ, std::vector<double>& w) {
        const std::size_t E = env.nstages();
        const std::vector<double> pi0 = env_init_prob(o, E);
        const double inf = std::numeric_limits<double>::infinity();
        const double hz = env_horizon(o);
        occ.assign(E, 0.0);
        w.assign(E, 0.0);
        if (std::isfinite(hz)) {
            // [A I; 0 0] * T, whose exponential carries int_0^T exp(A t) dt in its
            // upper-right block; A is the generator of the exponential arcs.
            Matrix<double> X(2 * E, 2 * E, 0.0);
            for (std::size_t e = 0; e < E; ++e) {
                double out = 0.0;
                for (std::size_t h = 0; h < E; ++h) {
                    if (h == e) continue;
                    const EnvArc<T>& a = env.arc(e, h);
                    if (!a.enabled || a.dist.disabled) continue;
                    const double r = 1.0 / num_traits<T>::to_double(a.dist.mean);
                    X(e, h) = r * hz;
                    out += r;
                }
                X(e, e) = -out * hz;
                X(e, E + e) = hz;
            }
            const Matrix<double> Z = expm(X);
            for (std::size_t h = 0; h < E; ++h) {
                double s = 0.0;
                for (std::size_t e = 0; e < E; ++e) s += pi0[e] * Z(e, E + h);
                occ[h] = s;
                w[h] = s / hz;
            }
            return;
        }
        std::vector<char> absorb(E, 0);
        bool any = false;
        for (std::size_t e = 0; e < E; ++e) {
            absorb[e] = !(env.rate[e] > 0.0);
            any = any || absorb[e];
        }
        if (!any) {
            for (std::size_t e = 0; e < E; ++e) {
                w[e] = env.prob_env[e];
                occ[e] = w[e] == 0.0 ? 0.0 : inf;
            }
            return;
        }
        // The mean-equivalent generator diag(lambda) (Pemb - I), restricted to
        // the transient stages.
        std::vector<std::size_t> tr;
        for (std::size_t e = 0; e < E; ++e)
            if (!absorb[e]) tr.push_back(e);
        const std::size_t nt = tr.size();
        Matrix<double> N(nt, nt, 0.0);
        for (std::size_t a = 0; a < nt; ++a)
            for (std::size_t b = 0; b < nt; ++b)
                N(a, b) = -env.rate[tr[a]] * (env.pemb(tr[a], tr[b]) - (a == b ? 1.0 : 0.0));
        Matrix<double> Ninv;
        double rc = 0.0;
        try {
            Ninv = inverse(N);
            const double n1 = norm1(N), ni = norm1(Ninv);
            rc = (n1 > 0.0 && std::isfinite(ni)) ? 1.0 / (n1 * ni) : 0.0;
        } catch (const NumericError&) {
            rc = 0.0;
        }
        if (!(rc >= lang::GlobalConstants::FineTol))
            throw UnsupportedError(ErrorKind{"Structure"},
                "SolverENV: method 'occupation' with timespan(2) = Inf needs every non-absorbing "
                "stage to reach an absorbing one: some stages form a closed class. Set a finite "
                "horizon, timespan [0,T].");
        for (std::size_t b = 0; b < nt; ++b) {
            double s = 0.0;
            for (std::size_t a = 0; a < nt; ++a) s += pi0[tr[a]] * Ninv(a, b);
            occ[tr[b]] = s;
        }
        for (std::size_t h = 0; h < E; ++h) {
            if (!absorb[h]) continue;
            double s = pi0[h];
            for (std::size_t a = 0; a < nt; ++a)
                s += occ[tr[a]] * env.rate[tr[a]] * env.pemb(tr[a], h);
            w[h] = s;
            if (s > 0.0) occ[h] = inf;
        }
    }
 
    /**
     * `occupationRefusal_`: a finite horizon integrates exp(A t), which only an
     * environment of exponential arcs has; T = Inf reads the embedded chain and
     * mean holding times, which any semi-Markov stage has. The initial stage
     * law must be a distribution over the stages.
     */
    static void occupation_refusal(const Environment<T>& env, const EnvOptions& o) {
        const std::size_t E = env.nstages();
        if (std::isfinite(env_horizon(o))) {
            for (std::size_t e = 0; e < E; ++e) {
                if (env.is_kernel_stage(e))
                    throw UnsupportedError(ErrorKind{"Structure"},
                        "Stage " + std::to_string(e + 1) + " has kernel transitions, but method "
                        "'occupation' over a finite timespan needs exponential transitions to "
                        "integrate the stage law. Use timespan(2) = Inf.");
                for (std::size_t h = 0; h < E; ++h) {
                    const EnvArc<T>& a = env.arc(e, h);
                    if (h == e || !a.enabled || a.dist.disabled ||
                        a.dist.type == lang::ProcessType::EXP)
                        continue;
                    throw UnsupportedError(ErrorKind{"Structure"},
                        "The environment transition from stage " + std::to_string(e + 1) +
                        " to " + std::to_string(h + 1) + " is a " +
                        lang::process_to_text(a.dist.type) +
                        ": method 'occupation' over a finite timespan needs exponential "
                        "transitions. Expand it into exponential stages, or use timespan(2) = "
                        "Inf.");
                }
            }
        }
        const std::vector<double> pi0 = env_init_prob(o, E);
        double tot = 0.0;
        bool neg = false;
        for (double v : pi0) {
            tot += v;
            if (!(v >= 0.0)) neg = true;
        }
        if (pi0.size() != E || neg || std::fabs(tot - 1.0) > lang::GlobalConstants::CoarseTol)
            throw InputError("options.config.env_init must be a probability vector over the " +
                             std::to_string(E) + " stages of the environment.");
    }
 
    /** `envInitProb_`: `options.config.env_init`, default the first stage. */
    static std::vector<double> env_init_prob(const EnvOptions& o, std::size_t E) {
        if (!o.env_init.empty()) return o.env_init;
        std::vector<double> pi0(E, 0.0);
        pi0[0] = 1.0;
        return pi0;
    }
 
    /** `envHorizon_`: the length of SolverENV's timespan, Inf when unbounded. */
    static double env_horizon(const EnvOptions& o) {
        if (std::isfinite(o.env_timespan_start) && std::isfinite(o.env_timespan_end))
            return o.env_timespan_end - o.env_timespan_start;
        return std::numeric_limits<double>::infinity();
    }
 
private:
    /** The 1-norm, the column-sum maximum `rcond` is taken in. */
    static double norm1(const Matrix<double>& A) {
        double m = 0.0;
        for (std::size_t j = 0; j < A.cols(); ++j) {
            double s = 0.0;
            for (std::size_t i = 0; i < A.rows(); ++i) s += std::fabs(A(i, j));
            if (!std::isfinite(s)) return std::numeric_limits<double>::infinity();
            m = std::max(m, s);
        }
        return m;
    }
 
    void init() {
        if (opt.method != "avg" && opt.method != "dec" && opt.method != "occupation")
            throw UnsupportedError(ErrorKind{"Method"}, "SolverENV limit: '" + opt.method +
                                   "' is not a closed-form environment limit; the two are 'avg' "
                                   "(fast environment) and 'dec' (slow environment), and "
                                   "'occupation' is the occupation-weighted 'dec'");
        dec_ = (opt.method == "dec");
        occ_ = (opt.method == "occupation");
        if (!stage_fn_ && opt.stage_solver != "fluid" && opt.stage_solver != "ctmc" &&
            opt.stage_solver != "mva" && opt.stage_solver != "nc")
            throw UnsupportedError(ErrorKind{"Option"},
                "SolverENV limit: stage solver '" + opt.stage_solver +
                "' is not available; the limits solve each stage in STEADY STATE, by the fluid "
                "analyzer, the enumerated CTMC, MVA or NC. Hand a stage solver to the "
                "EnvStageAvgFn constructor to run another one");
        if (!std::is_same<T, double>::value)
            throw UnsupportedError(ErrorKind{"Arithmetic"},
                "SolverENV limit: a fluid stage integrates its drift with LSODA, which is double "
                "only; rerun with --arith double");
 
        // The two limits read the stage NETWORKS directly -- `dec` solves each
        // in steady state, `avg` builds one network at the probEnv-weighted
        // rates -- and a layered stage has no station rate table to weight. The
        // reference has no branch for it either: `solveEnvLimit` reaches for
        // `self.ensemble{1}.nodes`, which a LayeredNetwork does not carry.
        // `occupation` reads a layered stage in its block layout instead.
        if (!occ_) envObj.reject_lqn_stages(
            "SolverENV limit",
            "the fast/slow limits read a stage's station rates directly -- 'dec' solves each "
            "stage network in steady state and 'avg' builds one network at the probEnv-weighted "
            "rates -- and a layered model has no such rate table, only the layers SolverLN "
            "derives from it");
        if (occ_) {
            occupation_refusal(envObj, opt);
        } else {
            envObj.refuse_semi_markov(opt.method);
            envObj.refuse_absorbing(opt.method);
        }
        envObj.init();
        const std::size_t E = envObj.nstages();
        lqn_ = envObj.has_lqn_stages();
        if (lqn_) {
            for (std::size_t e = 0; e < E; ++e)
                if (!envObj.is_lqn(e))
                    throw InputError(
                        "SolverENV occupation: stage " + std::to_string(e + 1) + " is a Network and "
                        "another stage a LayeredNetwork; the metrics are blended entrywise, so "
                        "every stage must be read in one layout");
            M = envObj.stage(0).lqn_model.nidx;
            K = 1;
            return;
        }
        M = envObj.stage(0).model.nstations;
        K = envObj.stage(0).model.nclasses;
        for (std::size_t e = 1; e < E; ++e)
            if (envObj.stage(e).model.nstations != M || envObj.stage(e).model.nclasses != K)
                throw InputError(
                    "SolverENV limit: every stage must have the same stations and classes; the "
                    "metrics are blended entrywise across them");
    }
 
public:
    /**
     * The RATE accumulator behind the cache blend, `SolverENV.accumCacheMetric`. Public, since
     * the mean-field coupling's `substituteCacheBlend_` weights by the same rule.
     *
     * Keyed by the Cache node's NAME, never by its position: the struct's node
     * order is not the declaration order, and indexing by position is what once
     * made a host write cache results onto a Sink (see `cache_metrics.h`).
     */
    class CacheAccum {
    public:
        /** One stage's contribution, weighted by its arrival into each cache. */
        void add(const qn::NetworkStruct<T>& sn, const EnvStageAvg<T>& stage, double w) {
            const std::vector<double> arv_all = cache_arrival(sn, stage.TN);
            for (std::size_t c = 0; c < stage.cache.caches.size(); ++c) {
                const solvers::CacheNodeMetrics<T>& cm = stage.cache.caches[c];
                Entry& en = entry(cm);
                const std::size_t K = std::max(cm.hitprob.size(), cm.missprob.size());
                std::vector<double> arv(std::max(K, en.arv.size()), 0.0);
                for (std::size_t r = 0; r < arv.size(); ++r) {
                    const std::size_t ind = cm.node;  // 1-based node index
                    const std::size_t off = (ind >= 1 ? (ind - 1) * sn.nclasses + r : arv_all.size());
                    if (off < arv_all.size()) arv[r] = w * arv_all[off];
                }
                grow(en.arv, arv.size());
                for (std::size_t r = 0; r < arv.size(); ++r) en.arv[r] += arv[r];
                accum(en.hit, cm.hitprob, arv);
                accum(en.miss, cm.missprob, arv);
                accum(en.dhit, cm.delayedprob, arv);
                accum_rows(en.hitl, cm.hitproblist, arv);
            }
        }
 
        /** Divide the accumulated rates by the accumulated arrivals, ONCE. */
        solvers::CacheMetrics<double> finish() const {
            solvers::CacheMetrics<double> out;
            for (std::size_t c = 0; c < order_.size(); ++c) {
                const Entry& en = order_[c];
                solvers::CacheNodeMetrics<double> cm;
                cm.node = en.node;
                cm.name = en.name;
                cm.itemcap = en.itemcap;
                cm.nitems = en.nitems;
                cm.hitprob = divide(en.hit, en.arv);
                cm.missprob = divide(en.miss, en.arv);
                cm.delayedprob = divide(en.dhit, en.arv);
                if (en.hitl.rows() > 0) {
                    cm.hitproblist = Matrix<double>(en.hitl.rows(), en.hitl.cols(), nan_());
                    for (std::size_t r = 0; r < en.hitl.rows(); ++r) {
                        const double d = r < en.arv.size() ? en.arv[r] : 0.0;
                        for (std::size_t l = 0; l < en.hitl.cols(); ++l)
                            cm.hitproblist(r, l) = d > 0.0 ? en.hitl(r, l) / d : nan_();
                    }
                }
                out.caches.push_back(cm);
            }
            return out;
        }
 
    private:
        struct Entry {
            std::size_t node = 0;
            std::string name;
            std::vector<double> itemcap;
            std::size_t nitems = 0;
            std::vector<double> arv;
            /** EMPTY stays empty: a quantity no stage measured must not become 0. */
            std::vector<double> hit, miss, dhit;
            Matrix<double> hitl;
        };
 
        static double nan_() { return std::numeric_limits<double>::quiet_NaN(); }
 
        static void grow(std::vector<double>& v, std::size_t n) {
            if (v.size() < n) v.resize(n, 0.0);
        }
 
        Entry& entry(const solvers::CacheNodeMetrics<T>& cm) {
            for (std::size_t i = 0; i < order_.size(); ++i)
                if (order_[i].name == cm.name) return order_[i];
            Entry en;
            en.node = cm.node;
            en.name = cm.name;
            en.itemcap = cm.itemcap;
            en.nitems = cm.nitems;
            order_.push_back(en);
            return order_.back();
        }
 
        /** One quantity: an ABSENT ratio contributes nothing and leaves it absent. */
        static void accum(std::vector<double>& a, const std::vector<T>& ratio,
                          const std::vector<double>& arv) {
            if (ratio.empty()) return;
            grow(a, ratio.size());
            for (std::size_t r = 0; r < ratio.size() && r < arv.size(); ++r) {
                const double v = num_traits<T>::to_double(ratio[r]);
                if (std::isfinite(v)) a[r] += arv[r] * v;
            }
        }
 
        static void accum_rows(Matrix<double>& a, const Matrix<T>& hl,
                               const std::vector<double>& arv) {
            if (hl.rows() == 0 || hl.cols() == 0) return;
            if (a.rows() == 0) a = Matrix<double>(hl.rows(), hl.cols(), 0.0);
            for (std::size_t r = 0; r < hl.rows() && r < a.rows(); ++r) {
                const double w = r < arv.size() ? arv[r] : 0.0;
                for (std::size_t l = 0; l < hl.cols() && l < a.cols(); ++l) {
                    const double v = num_traits<T>::to_double(hl(r, l));
                    if (std::isfinite(v)) a(r, l) += w * v;
                }
            }
        }
 
        /**
         * A class with no arrival anywhere divides by nothing and STAYS NaN: its
         * ratio is undefined, and 0 would read as "never hits".
         */
        static std::vector<double> divide(const std::vector<double>& a,
                                          const std::vector<double>& d) {
            std::vector<double> r;
            if (a.empty()) return r;
            r.assign(a.size(), nan_());
            for (std::size_t i = 0; i < a.size(); ++i)
                if (i < d.size() && d[i] > 0.0) r[i] = a[i] / d[i];
            return r;
        }
 
        /**
         * The per-class arrival INTO EACH NODE, flattened node-major.
         *
         * The weight is the cache's own arrival and not the Source throughput,
         * because the latter is 0 on a closed cache model. The station table
         * handed in is all zeros on purpose: only the Cache rows are read, and
         * those come from the chain reference throughput and the node visits,
         * never from a station row.
         */
        static std::vector<double> cache_arrival(const qn::NetworkStruct<T>& sn,
                                                 const Matrix<double>& TN) {
            const std::size_t I = sn.nodes.size(), R = sn.nclasses, M = sn.nstations;
            std::vector<double> out(I * R, 0.0);
            if (TN.rows() == 0) return out;
            Matrix<T> TNt(TN.rows(), TN.cols(), num_traits<T>::from_int(0));
            for (std::size_t i = 0; i < TN.rows(); ++i)
                for (std::size_t r = 0; r < TN.cols(); ++r)
                    TNt(i, r) = num_traits<T>::from_double(TN(i, r));
            const Matrix<T> AN(M, R, num_traits<T>::from_int(0));
            const Matrix<T> ANn = api::sn_get_node_arvr_from_tput<T>(sn, TNt, AN);
            if (ANn.rows() != I) return out;
            for (std::size_t ind = 0; ind < I; ++ind)
                for (std::size_t r = 0; r < R; ++r) {
                    const double v = num_traits<T>::to_double(ANn(ind, r));
                    out[ind * R + r] = std::isfinite(v) ? v : 0.0;
                }
            return out;
        }
 
        std::vector<Entry> order_;
    };
 
    /**
     * One stage in steady state under the enumerated CTMC at `cutoff`, with its cache split,
     * which is exact (read off the chain). An ensemble built on SolverCTMC stages.
     */
    static EnvStageAvg<T> ctmc_stage_steady(const qn::NetworkStruct<T>& sn, double cutoff) {
        ctmc::CtmcOptions co;
        co.cutoff = cutoff;
        const mva::AvgResult<T> r = ctmc::solver_ctmc_run_analyzer<T>(sn, co);
        EnvStageAvg<T> out;
        out.QN = to_double(r.QN);
        out.UN = to_double(r.UN);
        out.TN = to_double(r.TN);
        out.cache = r.cache;
        return out;
    }
 
private:
    /** A stage table widened to double, which is what the limits blend in. */
    static Matrix<double> to_double(const Matrix<T>& A) {
        Matrix<double> out(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) out(i, j) = num_traits<T>::to_double(A(i, j));
        return out;
    }
 
    /** The slow-environment limit: each stage on its own, blended by probEnv. */
    EnvLimitSolution solve_dec() {
        const std::size_t E = envObj.nstages();
        EnvLimitSolution out;
        out.method = "dec";
        out.prob_env = envObj.prob_env;
        out.QN = Matrix<double>(M, K, 0.0);
        out.UN = Matrix<double>(M, K, 0.0);
        out.TN = Matrix<double>(M, K, 0.0);
        out.QStage.assign(E, Matrix<double>(M, K, 0.0));
        out.UStage.assign(E, Matrix<double>(M, K, 0.0));
        out.TStage.assign(E, Matrix<double>(M, K, 0.0));
        CacheAccum acc;
        for (std::size_t e = 0; e < E; ++e) {
            const EnvStageAvg<T> s = stage_steady_state(envObj.stage(e).model);
            const double p = envObj.prob_env[e];
            acc.add(envObj.stage(e).model, s, p);
            for (std::size_t i = 0; i < M; ++i)
                for (std::size_t r = 0; r < K; ++r) {
                    out.QStage[e](i, r) = s.QN(i, r);
                    out.UStage[e](i, r) = s.UN(i, r);
                    out.TStage[e](i, r) = s.TN(i, r);
                    out.QN(i, r) += p * s.QN(i, r);
                    out.UN(i, r) += p * s.UN(i, r);
                    out.TN(i, r) += p * s.TN(i, r);
                }
        }
        out.cache = acc.finish();
        return out;
    }
 
    /**
     * `solveEnvOccupation`: each stage of nonzero weight in steady state, as
     * `dec` solves it, blended by the occupation weights. A NaN or Inf cell is
     * an ABSENT metric (`local_absent_zero`) and adds nothing to the blend.
     */
    EnvLimitSolution solve_occupation() {
        const std::size_t E = envObj.nstages();
        EnvLimitSolution out;
        out.method = "occupation";
        stage_occupation(envObj, opt, out.occupation, out.occupation_weight);
        out.prob_env = out.occupation_weight;
        out.QN = Matrix<double>(M, K, 0.0);
        out.UN = Matrix<double>(M, K, 0.0);
        out.TN = Matrix<double>(M, K, 0.0);
        out.QStage.assign(E, Matrix<double>(M, K, 0.0));
        out.UStage.assign(E, Matrix<double>(M, K, 0.0));
        out.TStage.assign(E, Matrix<double>(M, K, 0.0));
        for (std::size_t e = 0; e < E; ++e) {
            const double w = out.occupation_weight[e];
            if (w == 0.0) continue;
            const EnvStageAvg<T> s = lqn_ ? lqn_stage_steady_state(e) : stage_steady_state(envObj.stage(e).model);
            for (std::size_t i = 0; i < M; ++i)
                for (std::size_t r = 0; r < K; ++r) {
                    out.QStage[e](i, r) = s.QN(i, r);
                    out.UStage[e](i, r) = s.UN(i, r);
                    out.TStage[e](i, r) = s.TN(i, r);
                    if (std::isfinite(s.QN(i, r))) out.QN(i, r) += w * s.QN(i, r);
                    if (std::isfinite(s.UN(i, r))) out.UN(i, r) += w * s.UN(i, r);
                    if (std::isfinite(s.TN(i, r))) out.TN(i, r) += w * s.TN(i, r);
                }
        }
        return out;
    }
 
    /**
     * A LAYERED stage in steady state: its own `SolverLN` run to the fixed point
     * and read over its LQN elements, as the reference reads `SolverLN.getAvg`
     * (row i is element i+1). An undefined cell is NaN, an absent metric that
     * `solve_occupation` leaves out of the blend.
     */
    EnvStageAvg<T> lqn_stage_steady_state(std::size_t e) const {
        const lqn::LqnStruct<T>& model = envObj.stage(e).lqn_model;
        if (model.nidx != M)
            throw InputError(
                "SolverENV occupation: stage " + std::to_string(e + 1) + " has " +
                std::to_string(model.nidx) + " LQN elements where stage 1 has " + std::to_string(M) +
                "; the layered stages are blended element by element, so they must have the same "
                "hosts, tasks, entries and activities");
        ln::SolverLN<T> s(model, opt.lqn);
        const ln::LnSolution<T> r = s.get_ensemble_avg();
        const double nan = std::numeric_limits<double>::quiet_NaN();
        EnvStageAvg<T> out;
        out.QN = Matrix<double>(M, 1, nan);
        out.UN = Matrix<double>(M, 1, nan);
        out.TN = Matrix<double>(M, 1, nan);
        for (std::size_t i = 1; i <= M; ++i) {
            if (r.defined_Q[i]) out.QN(i - 1, 0) = num_traits<T>::to_double(r.QN[i]);
            if (r.defined_U[i]) out.UN(i - 1, 0) = num_traits<T>::to_double(r.UN[i]);
            if (r.defined_T[i]) out.TN(i - 1, 0) = num_traits<T>::to_double(r.TN[i]);
        }
        return out;
    }
 
    /** The fast-environment limit: one rate-averaged model, solved once. */
    EnvLimitSolution solve_avg() {
        EnvLimitSolution out;
        out.method = "avg";
        out.prob_env = envObj.prob_env;
        const qn::NetworkStruct<T> avg = rate_averaged_model();
        const EnvStageAvg<T> s = stage_steady_state(avg);
        // One model, so the weight is the arrival itself and the normalisation
        // returns exactly what the single solve reported.
        CacheAccum acc;
        acc.add(avg, s, 1.0);
        out.cache = acc.finish();
        out.QN = Matrix<double>(M, K, 0.0);
        out.UN = Matrix<double>(M, K, 0.0);
        out.TN = Matrix<double>(M, K, 0.0);
        for (std::size_t i = 0; i < M; ++i)
            for (std::size_t r = 0; r < K; ++r) {
                out.QN(i, r) = s.QN(i, r);
                out.UN(i, r) = s.UN(i, r);
                out.TN(i, r) = s.TN(i, r);
            }
        return out;
    }
 
    /**
     * One stage, in steady state: the injected solver when there is one, else
     * the fluid analyzer, in the form that also REPORTS ITS CACHE SPLIT when the
     * stage holds a Cache.
     *
     * WHICH ENTRY POINT IS A MODEL QUESTION AND NOT A PREFERENCE, because the
     * two resolve `default` differently and therefore answer with DIFFERENT
     * CLOSURES. `solver_fluid` is `solver_fluid_analyzer.m`: it hands the method
     * to `fluid_dispatch`, which takes `default` to the matrix method (to
     * `closing` on a DPS model). `solver_fluid_run_analyzer` is the port of
     * `runAnalyzer`'s resolution and re-resolves `default` through the
     * minnormal/dae ladder. On a stage sitting at its saturation knee the two
     * disagree outright -- N/(Z+D) = 1/D is a continuum of first-order
     * equilibria, and which point a closure picks is the whole answer -- so
     * routing every stage through the runner silently moved the limits off the
     * matrix closure they are pinned against.
     *
     * A CACHE STAGE IS THE ONE THAT MUST TAKE THE RUNNER, and for the reference's
     * own reason rather than for a reporting convenience: `runAnalyzer` resolves
     * a Cache model to `rmf`, which lives in `fluid_cacheqn.h` and which
     * `fluid_dispatch` refuses BY NAME (it cannot call it without a cyclic
     * include). That route is also the one that reports the split, which is the
     * silence this file used to refuse a cache model over.
     */
    EnvStageAvg<T> stage_steady_state(const qn::NetworkStruct<T>& sn) const {
        if (stage_fn_) return stage_fn_(sn);
        if (opt.stage_solver == "ctmc") return ctmc_stage_steady(sn, opt.stage_cutoff);
        if (opt.stage_solver == "mva" || opt.stage_solver == "nc") {
            // `SolverENV(env, @SolverMVA)` / `@SolverNC`: the stage's own getAvg,
            // at the options the ensemble was built with, from no warm start.
            const mva::AvgResult<T> r =
                opt.stage_solver == "mva"
                    ? mva::solver_mva_run_analyzer<T>(sn, opt.stage_mva, Matrix<T>())
                    : nc::solver_nc_run_analyzer<T>(sn, opt.stage_nc);
            EnvStageAvg<T> out;
            out.QN = to_double(r.QN);
            out.UN = to_double(r.UN);
            out.TN = to_double(r.TN);
            out.cache = r.cache;
            return out;
        }
        fluid::FluidOptions fo = opt.stage;
        // A limit starts from nothing carried over -- there is no entry state to
        // seed with, which is the whole content of the approximation.
        fo.init_sol.clear();
        EnvStageAvg<T> out;
        const fluid::FluidSolution s =
            fluid::detail::fluid_has_cache(sn)
                ? fluid::solver_fluid_run_analyzer<T>(sn, fo, nullptr, nullptr, &out.cache)
                : fluid::solver_fluid<T>(sn, fo);
        out.QN = s.QN;
        out.UN = s.UN;
        out.TN = s.TN;
        return out;
    }
 
    /**
     * `buildRateAveragedModel`: stage 1's network with every MODULATED rate
     * replaced by its probEnv-weighted mean, as an exponential.
     *
     * The refresh chain is rerun because `rates` and everything derived from it
     * are read off the service table; editing the table alone would leave the
     * struct describing stage 1.
     */
    qn::NetworkStruct<T> rate_averaged_model() const {
        const std::size_t E = envObj.nstages();
        qn::NetworkStruct<T> sn = envObj.stage(0).model;
        std::vector<double> r(E, 0.0);
        for (std::size_t i = 0; i < M; ++i) {
            const lang::NodeType nt = sn.stations[i].nodetype;
            if (nt != lang::NodeType::Source && nt != lang::NodeType::Queue &&
                nt != lang::NodeType::Delay)
                continue;
            for (std::size_t k = 0; k < K; ++k) {
                bool ok = true;
                for (std::size_t e = 0; e < E && ok; ++e) {
                    const qn::NetworkStruct<T>& se = envObj.stage(e).model;
                    // The reference's `any(isnan(r)) || any(r<=0)`: a rate this
                    // port reports as disabled is MATLAB's NaN, and either way
                    // the pair is left as configured.
                    if (se.disabled[i][k]) {
                        ok = false;
                        break;
                    }
                    r[e] = num_traits<T>::to_double(se.rates(i, k));
                    if (!(r[e] > 0.0) || !std::isfinite(r[e])) ok = false;
                }
                if (!ok) continue;
                double lo = r[0], hi = r[0], avg = 0.0;
                for (std::size_t e = 0; e < E; ++e) {
                    lo = std::min(lo, r[e]);
                    hi = std::max(hi, r[e]);
                    avg += envObj.prob_env[e] * r[e];
                }
                // Not modulated: keep the original distribution, which may carry
                // a shape the exponential of its mean would throw away.
                if (hi - lo <= 1e-12 * std::max(1.0, hi)) continue;
                sn.set_service(i + 1, k + 1,
                               lang::Distrib<T>::exp_rate(num_traits<T>::from_double(avg)));
            }
        }
        sn.refresh_struct();
        return sn;
    }
 
    Environment<T>& envObj;
    EnvOptions opt;
    EnvStageAvgFn<T> stage_fn_;
    std::size_t M = 0, K = 0;
    bool dec_ = false;
    bool occ_ = false;
    bool lqn_ = false;  ///< `occupation` over layered stages
};
 
/** `solveEnvLimit` on the original stages. */
template <class T>
EnvLimitSolution solver_env_limit(Environment<T>& e, const EnvOptions& o) {
    SolverEnvLimit<T> s(e, o);
    return s.solve();
}
 
/**
 * `getStageOccupation`: method `occupation`'s per-stage expected occupation and
 * blend weight. It initialises the environment and solves no stage, so it is
 * the cheap reading; `solver_env_limit` reports the same two vectors.
 */
template <class T>
void stage_occupation(Environment<T>& e, const EnvOptions& o, std::vector<double>& occ,
                      std::vector<double>& w) {
    if (o.method != "occupation")
        throw InputError("getStageOccupation is defined for SolverENV method 'occupation'.");
    SolverEnvLimit<T>::occupation_refusal(e, o);
    e.init();
    SolverEnvLimit<T>::stage_occupation(e, o, occ, w);
}
 
/** `solveEnvLimit` with the CALLING solver running every stage. */
template <class T>
EnvLimitSolution solver_env_limit(Environment<T>& e, const EnvOptions& o,
                                  EnvStageAvgFn<T> stage_fn) {
    SolverEnvLimit<T> s(e, o, stage_fn);
    return s.solve();
}
 
}  // namespace env
}  // namespace line
 
#endif  // LINE_SOLVERS_ENV_SOLVER_ENV_LIMIT_H