environment.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_LANG_QN_ENVIRONMENT_H
#define LINE_LANG_QN_ENVIRONMENT_H
 
/**
 * @file
 * @ingroup line_lang
 * A random environment: a port of `matlab/src/lang/Environment.m`, restricted
 * to what `SolverENV` reads out of it.
 *
 * WHAT IT MODELS. A queueing network whose PARAMETERS change over time because
 * the world around it does: a server that breaks down and is repaired, a
 * workload that has a day phase and a night phase, a system that degrades and
 * is reset. Each of those is a STAGE, holding its own complete network, and the
 * environment is a semi-Markov process over the stages. The point is that the
 * network never restarts empty at a switch -- the jobs in it at the moment of
 * the switch are carried into the next stage, which is what couples the stages
 * and what makes this more than solving each stage separately.
 *
 * THE HOLDING TIME IS A COMPETING RISK, and this is the part worth reading the
 * reference for. Every enabled transition e -> h has its own distribution, and
 * they all run at once: the stage ends when the FIRST of them fires, and which
 * one fired decides the next stage. `init()` therefore superposes the outgoing
 * transitions of e by a Kronecker sum into a single marked process
 * `hold_time[e]`, whose per-mark rates give the embedded jump chain `Pemb`.
 * The collapse that follows -- each block's row sums moved into its FIRST
 * column -- is the reference's, and it makes the superposed process restart in
 * phase one after every jump, i.e. renewal at each stage entry.
 *
 * WHAT init() PRODUCES, and what the solver consumes:
 *   proc[e][h]   the e -> h transition process, whose CDF weights the stage
 *                transient when computing the metrics AT a switch to h
 *   hold_time[e] the superposed holding time of stage e, whose CDF weights the
 *                transient when computing the metrics OVER the whole stage
 *   prob_env[e]  the stationary probability of being in stage e
 *   prob_orig    prob_orig(h, e) = P(the previous stage was h | now entering e)
 *
 * A SEMI-MARKOV STAGE has no such MAP form: a stage with a non-Markovian arc
 * (Det, Gamma, Weibull, ...) or with KERNEL arcs (`add_kernel_transition`, a
 * jump probability and a conditional sojourn law per arc). `init()` then reads
 * its embedded jump row and mean holding time off the arc CDFs by
 * `mc::smp_stage_quad`, leaves `hold_time[e]` and `proc[e][h]` as the disabled
 * pair, and `hold_time_cdf` / `hold_time_mean` evaluate the holding-time law
 * from the arcs. Only SolverENV method `smp` solves such an environment.
 *
 * A DISABLED transition is the 1 x 1 zero pair, exactly as in the reference:
 * krons(A, 0) leaves A unchanged, so a disabled arc costs nothing and needs no
 * special case anywhere below.
 *
 * NODE BREAKDOWNS are the one stage pattern the reference gives a name to
 * (`addNodeBreakdown` / `addNodeRepair` / `addNodeFailureRepair`): the UP stage
 * holds the base network, the DOWN_<node> stage holds the same network with one
 * node's service replaced by its degraded distribution, and the two arcs
 * between them are the time to failure and the time to repair. It is a macro
 * over `set_stage` and `add_transition` and nothing more, EXCEPT for the pair of
 * queue-length reset policies it attaches, which are the only part of a
 * breakdown that the expanded stages cannot express -- hence `NodeFailure`,
 * which records them beside the stages so an environment read back from
 * `model.json` is the same model it was written from.
 */
 
#include <cctype>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <functional>
#include <limits>
#include <string>
#include <vector>
 
#include "line/api/mam/map_cdf.h"
#include "line/api/mam/map_moment.h"
#include "line/api/mam/mmap_lambda.h"
#include "line/api/mc/ctmc_courtois.h"
#include "line/api/mc/smp_stage_cif.h"
#include "line/api/mc/smp_stage_quad.h"
#include "line/lang/distribution.h"
#include "line/lang/lang_types.h"
#include "line/lang/lqn/lqn_struct.h"
#include "line/lang/prior.h"
#include "line/lang/qn/network_struct.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace env {
 
/**
 * The reset policy of a transition, `resetFun` in the reference.
 *
 * It maps the mean queue lengths at the moment of the switch onto the mean
 * queue lengths the next stage starts from. Identity means the jobs are simply
 * carried over; zero means the buffer is flushed on the switch.
 */
using ResetMarginal = std::function<Matrix<double>(const Matrix<double>&)>;
 
/**
 * The two NAMED reset policies of `Environment.resolveResetPolicy`, which are
 * the only ones the JSON interchange can carry (a function handle is written as
 * `custom` and warned about by both writers, never reloaded).
 *
 * `keep` resolves to the EMPTY function rather than to an explicit identity,
 * because empty is how this port spells identity everywhere a reset is read:
 * `SolverEnv::post` skips the call, and the compression's macro-arc fold
 * compares resets only by whether one is PRESENT, so an explicit identity on
 * one arc and nothing on another would be refused as a disagreement although
 * the two mean the same thing.
 */
inline ResetMarginal env_reset_policy(const std::string& name) {
    std::string low;
    low.reserve(name.size());
    for (std::size_t i = 0; i < name.size(); ++i)
        low.push_back(static_cast<char>(std::tolower(static_cast<unsigned char>(name[i]))));
    if (low == "keep") return ResetMarginal();
    if (low == "clear")
        return [](const Matrix<double>& q) { return Matrix<double>(q.rows(), q.cols(), 0.0); };
    throw InputError("Environment: unknown reset policy '" + name +
                     "'. The serializable policies are 'keep' (carry the queue lengths across "
                     "the switch) and 'clear' (empty the queues); a custom policy is a function "
                     "and is installed through add_transition, not by name");
}
 
/**
 * `Environment.nodeFailures{k}`: the declarative record of one node breakdown.
 *
 * The expanded stages and transitions already carry the STRUCTURE of a
 * breakdown; what only this record carries is the pair of reset policies, which
 * decide what the next stage starts from and therefore change the numbers. It
 * is kept on the environment for the same reason the reference keeps it: so an
 * environment read from `model.json` in its expanded form can be given back its
 * policies, and so it serializes out again as the same model.
 */
template <class T>
struct NodeFailure {
    std::string node;            ///< the node that breaks down
    lang::Distrib<T> breakdown;  ///< time to failure, the UP -> DOWN transition
    lang::Distrib<T> down_service;  ///< the node's service while it is down
    lang::Distrib<T> repair;     ///< time to repair, the DOWN -> UP transition
    bool has_repair = false;     ///< false for a breakdown with no repair arc
    std::string breakdown_reset = "keep";
    std::string repair_reset;    ///< empty when there is no repair
};
 
/**
 * One stage: a name, a category, and the model in force while it lasts.
 *
 * THE MODEL IS ONE OF TWO KINDS, and which one is the stage's own property
 * rather than the environment's: a flat `qn::NetworkStruct`, or a
 * `lqn::LqnStruct` -- a LayeredNetwork, whose stations and classes are those of
 * the LAYERS SolverLN builds out of it. The reference draws no distinction at
 * this level either (`Environment.addStage` stores whatever model it is handed
 * and `SolverENV` branches on its class), and the two fields are kept apart
 * rather than unified because nothing about them is shared below the name: a
 * flat stage is integrated directly, a layered one through a whole fixed point
 * over its layers.
 *
 * `has_model` and `has_lqn` are mutually exclusive; `set_stage` and
 * `set_lqn_stage` each clear the other, so a stage re-declared as the other
 * kind cannot leave a stale twin behind for a later consumer to read.
 */
template <class T>
struct EnvStage {
    std::string name;
    std::string type;  ///< the stage category, informational only
    qn::NetworkStruct<T> model;
    bool has_model = false;
    /** The layered model, when this stage holds a LayeredNetwork. */
    lqn::LqnStruct<T> lqn_model;
    bool has_lqn = false;
    /** True when the stage carries a model of either kind. */
    bool declared() const { return has_model || has_lqn; }
};
 
/**
 * `resetEnvRatesFun` in the reference: the state-dependent environment rate.
 *
 * It is given the arc's CURRENT transition distribution and the exit metrics of
 * the stage the arc leaves -- mean queue lengths, utilizations and throughputs,
 * averaged over when that arc fires -- and returns the distribution the arc
 * should carry next. It is what makes the environment process depend on the
 * network it modulates, and `method = "statedep"` is what applies it.
 */
template <class T>
using ResetEnvRates =
    std::function<lang::Distrib<T>(const lang::Distrib<T>&, const Matrix<double>&,
                                   const Matrix<double>&, const Matrix<double>&)>;
 
/** One arc of the environment process. */
template <class T>
struct EnvArc {
    bool enabled = false;
    lang::Distrib<T> dist;   ///< the e -> h transition time
    ResetMarginal reset;     ///< empty means the identity
    ResetEnvRates<T> reset_rates;  ///< empty means the rate does not depend on the state
};
 
/**
 * The DOWN stage network of `addNodeBreakdown`: the base network with ONE
 * node's service replaced by its degraded distribution, for EVERY class.
 *
 * Every class, and not only the ones that were enabled there, is what the
 * reference does (`for c = 1:length(classes), nodes{nodeIdx}.setService(...)`),
 * so a class that was disabled at the node while it was up is served at the
 * degraded rate while it is down. The whole refresh chain is rerun afterwards
 * because `rates`, `scv` and the chain-derived tables are all read off the
 * service table; editing the table alone would leave the struct describing the
 * UP stage and the solver reading the DOWN one.
 */
template <class T>
qn::NetworkStruct<T> env_degraded_model(const qn::NetworkStruct<T>& base,
                                        const std::string& node_name,
                                        const lang::Distrib<T>& down_service) {
    qn::NetworkStruct<T> sn = base;
    std::size_t node = 0;
    for (std::size_t i = 0; i < sn.nodes.size(); ++i)
        if (sn.nodes[i].name == node_name) {
            node = i + 1;
            break;
        }
    if (node == 0)
        throw InputError("Environment: node '" + node_name +
                         "' is not in the base model, so it has no service to degrade");
    const std::size_t ist = sn.nodes[node - 1].station;
    if (ist == 0)
        throw InputError("Environment: node '" + node_name +
                         "' is not a station and carries no service distribution, so it cannot "
                         "break down into a degraded service");
    lang::Distrib<T> d = down_service;
    if (d.is_prior())
        lang::prior_refresh_moments(d);
    else
        lang::dist_refresh_moments(d);
    for (std::size_t r = 1; r <= sn.classes.size(); ++r) sn.set_service(ist, r, d);
    sn.refresh_struct();
    return sn;
}
 
template <class T>
class Environment {
public:
    Environment(const std::string& nm, std::size_t nstages)
        : name_(nm), stages_(nstages), arcs_(nstages, std::vector<EnvArc<T>>(nstages)),
          kernel_prob_(nstages, nstages, std::numeric_limits<double>::quiet_NaN()) {
        if (nstages == 0) throw InputError("Environment: a random environment needs a stage");
    }
 
    std::size_t nstages() const { return stages_.size(); }
    const std::string& name() const { return name_; }
 
    /** `addStage`: name the stage and give it its network. */
    void set_stage(std::size_t e, const std::string& nm, const std::string& type,
                   const qn::NetworkStruct<T>& model) {
        check(e);
        stages_[e].name = nm;
        stages_[e].type = type;
        stages_[e].model = model;
        stages_[e].has_model = true;
        stages_[e].has_lqn = false;
        stages_[e].lqn_model = lqn::LqnStruct<T>();
    }
 
    /**
     * `addStage` with a LayeredNetwork: the stage holds a LAYERED model.
     *
     * The environment itself does nothing with the difference -- probEnv,
     * probOrig and the holding times are read off the ARCS and never off a
     * stage model -- so the whole content of this overload is that the stage
     * records which kind of model it carries and `SolverEnv` runs the matching
     * stage solver. What the two kinds must still agree on is the (station,
     * class) SHAPE of the metrics being blended, and for a layered stage that
     * shape is the block-diagonal union of its layers; `SolverEnv::init` is
     * where the shapes are compared, because only there is a SolverLN built and
     * the layer blocks known.
     */
    void set_lqn_stage(std::size_t e, const std::string& nm, const std::string& type,
                       const lqn::LqnStruct<T>& model) {
        check(e);
        // compare LQN elements, not layers: a zero-capacity stage drops the layers
        // only its trapped reference tasks reach (SolverLN, ln_zero_capacity)
        for (std::size_t h = 0; h < stages_.size(); ++h)
            if (h != e && stages_[h].has_lqn &&
                (stages_[h].lqn_model.nidx != model.nidx || stages_[h].lqn_model.tshift != model.tshift))
                throw InputError(ErrorKind{"Model"},
                    "Unsupported feature. Random environment stages must map to layered networks "
                    "with the same hosts, tasks, entries and activities.");
        stages_[e].name = nm;
        stages_[e].type = type;
        stages_[e].lqn_model = model;
        stages_[e].has_lqn = true;
        stages_[e].has_model = false;
        stages_[e].model = qn::NetworkStruct<T>();
    }
 
    /** True when stage `e` holds a LayeredNetwork rather than a flat network. */
    bool is_lqn(std::size_t e) const {
        check(e);
        return stages_[e].has_lqn;
    }
 
    /** True when ANY stage holds a LayeredNetwork. */
    bool has_lqn_stages() const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (stages_[e].has_lqn) return true;
        return false;
    }
 
    /**
     * Refuse an environment carrying a LayeredNetwork stage, by name.
     *
     * Shared by every consumer that reads `stage(e).model` directly -- the
     * state-vector coupling, the closed-form limits, the compression -- so that
     * each says the same thing about the same gap rather than reading an EMPTY
     * NetworkStruct and reporting a confident answer about a model that is not
     * there. `who` names the caller and `why` says what about it needs a flat
     * stage.
     */
    void reject_lqn_stages(const std::string& who, const std::string& why) const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (stages_[e].has_lqn)
                throw UnsupportedError(ErrorKind{"LayeredNetwork"},
                    who + ": stage " + std::to_string(e + 1) + " ('" + stages_[e].name +
                    "') holds a LayeredNetwork, and " + why +
                    ". The mean-field coupling (method 'default' / 'meanfield') is the one that "
                    "solves a layered stage, through SolverLN");
    }
 
    /**
     * `addTransition`: enable e -> h with a distribution and a reset policy.
     *
     * A COMPETING arc: `d` is a latent clock started on entry to the stage, and
     * the first of the stage's clocks to ring picks the destination. Any law is
     * accepted; a non-Markovian one makes the environment semi-Markov, which
     * only SolverENV method `smp` solves.
     */
    void add_transition(std::size_t e, std::size_t h, const lang::Distrib<T>& d,
                        const ResetMarginal& reset = ResetMarginal()) {
        check(e);
        check(h);
        arcs_[e][h].enabled = true;
        arcs_[e][h].dist = d;
        arcs_[e][h].reset = reset;
        kernel_prob_(e, h) = std::numeric_limits<double>::quiet_NaN();
    }
 
    /**
     * `addKernelTransition`: a KERNEL arc of Q(e,h,t) = prob * F(t).
     *
     * On entry to e the next stage is h with probability `prob`, and the sojourn
     * is then distributed as `d`: the sojourn-time parameterization (Approach I
     * of Asanjarani et al.), beside `add_transition`'s competing clocks (Approach
     * II). A stage's arcs must be all kernel or all competing, and the
     * probabilities of a kernel stage must sum to one (`stage_arcs`). Only
     * SolverENV method `smp` solves an environment with kernel arcs.
     */
    void add_kernel_transition(std::size_t e, std::size_t h, double prob,
                               const lang::Distrib<T>& d,
                               const ResetMarginal& reset = ResetMarginal()) {
        if (!(prob > 0.0 && prob <= 1.0))
            throw InputError("Environment: the probability of a kernel transition must be a "
                             "scalar in (0,1].");
        add_transition(e, h, d, reset);
        kernel_prob_(e, h) = prob;
    }
 
    /** `kernelProb(e,h)`: the kernel arc's probability, NaN for a competing arc or none. */
    double kernel_prob(std::size_t e, std::size_t h) const {
        check(e);
        check(h);
        return kernel_prob_(e, h);
    }
 
    /** `isKernelStage`: true if stage e's arcs are kernel arcs. */
    bool is_kernel_stage(std::size_t e) const {
        check(e);
        for (std::size_t h = 0; h < stages_.size(); ++h)
            if (!std::isnan(kernel_prob_(e, h))) return true;
        return false;
    }
 
    /**
     * `isSemiMarkovStage`: true if stage e's holding time has no Markovian (MAP)
     * form, i.e. it has kernel arcs or a competing arc whose law is not
     * Markovian (the reference's `isa(arc, 'Markovian')`).
     */
    bool is_semi_markov_stage(std::size_t e) const {
        if (is_kernel_stage(e)) return true;
        for (std::size_t h = 0; h < stages_.size(); ++h) {
            const EnvArc<T>& a = arcs_[e][h];
            if (a.enabled && !a.dist.disabled && !arc_is_markovian(a.dist)) return true;
        }
        return false;
    }
 
    /**
     * `isAbsorbingStage`: true if stage e has no outgoing transition, so once
     * entered the environment stays there with an infinite holding time. A
     * self arc does not count as a way out.
     */
    bool is_absorbing_stage(std::size_t e) const {
        check(e);
        if (is_kernel_stage(e)) return false;
        for (std::size_t h = 0; h < stages_.size(); ++h)
            if (h != e && arcs_[e][h].enabled && !arcs_[e][h].dist.disabled) return false;
        return true;
    }
 
    /**
     * Refuse a stage with no model. Runs ahead of the structural refusals, since
     * an undeclared stage also has no arcs and would otherwise be reported as absorbing.
     */
    void require_declared() const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (!stages_[e].declared())
                throw InputError("Environment: stage " + std::to_string(e + 1) + " has no model");
    }
 
    /** True if some stage is absorbing, i.e. `init()` left probEnv empty. */
    bool is_absorbing() const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (is_absorbing_stage(e)) return true;
        return false;
    }
 
    /**
     * Refuse an absorbing environment for a coupling that reads probEnv: its
     * limit depends on the initial stage, so it has no single steady state.
     */
    void refuse_absorbing(const std::string& method) const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (is_absorbing_stage(e))
                throw UnsupportedError(ErrorKind{"Structure"},
                    "Stage " + std::to_string(e + 1) + " ('" + stages_[e].name +
                    "') has no outgoing transition, so the environment is absorbing and has no "
                    "single steady state: its limit depends on the initial stage. Method '" +
                    method + "' reads the stationary stage probabilities; use method "
                    "'occupation', with options.config.env_init naming the initial stage law.");
    }
 
    /** `isSemiMarkov`: true if any stage is semi-Markov. */
    bool is_semi_markov() const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (is_semi_markov_stage(e)) return true;
        return false;
    }
 
    /**
     * `getStageArcs`: the transitions out of stage e in the form
     * `mc::smp_stage_cif` takes. `arcs[h]` is the law of e -> h, DISABLED where
     * there is none (and for h == e), and `probs` is empty for competing arcs or
     * holds the kernel probabilities, validated as the reference does.
     */
    void stage_arcs(std::size_t e, std::vector<lang::Distrib<T>>& arcs,
                    std::vector<double>& probs) const {
        check(e);
        const std::size_t E = stages_.size();
        arcs.assign(E, lang::Distrib<T>());
        for (std::size_t h = 0; h < E; ++h)
            if (h != e && arcs_[e][h].enabled && !arcs_[e][h].dist.disabled)
                arcs[h] = arcs_[e][h].dist;
        probs.clear();
        if (!is_kernel_stage(e)) return;
        probs.assign(E, 0.0);
        double tot = 0.0;
        for (std::size_t h = 0; h < E; ++h) {
            if (!arcs[h].disabled) {
                if (std::isnan(kernel_prob_(e, h)))
                    throw InputError("Environment: stage " + std::to_string(e + 1) +
                                     " mixes kernel and competing transitions (the one to stage " +
                                     std::to_string(h + 1) +
                                     " was added by addTransition). Declare all of a stage's "
                                     "transitions with addKernelTransition, or none.");
                probs[h] = kernel_prob_(e, h);
                tot += probs[h];
            } else if (!std::isnan(kernel_prob_(e, h))) {
                throw InputError("Environment: the kernel transition " + std::to_string(e + 1) +
                                 "->" + std::to_string(h + 1) + " has no law.");
            }
        }
        if (std::fabs(tot - 1.0) > 1e-8) {
            char buf[64];
            std::snprintf(buf, sizeof buf, "%.12g", tot);
            throw InputError("Environment: the kernel probabilities out of stage " +
                             std::to_string(e + 1) + " sum to " + buf + ", not 1.");
        }
    }
 
    /**
     * The `methodRefusal` rule for every coupling but `smp`: refuse a kernel
     * stage, or a competing arc that is not Markovian, with the reference's
     * sentences. Those couplings read each stage's holding time as a (D0,D1)
     * race of its transitions, which neither has.
     */
    void refuse_semi_markov(const std::string& method) const {
        if (method == "smp") return;
        const std::size_t E = stages_.size();
        for (std::size_t e = 0; e < E; ++e) {
            if (is_kernel_stage(e))
                throw UnsupportedError(ErrorKind{"Structure"},
                    "Stage " + std::to_string(e + 1) +
                    " has kernel transitions (addKernelTransition), which only method 'smp' of "
                    "SolverENV solves: the other couplings read each stage's holding time as a "
                    "Markovian (D0,D1) race of its transitions.");
            for (std::size_t h = 0; h < E; ++h) {
                const EnvArc<T>& a = arcs_[e][h];
                if (!a.enabled || a.dist.disabled || arc_is_markovian(a.dist)) continue;
                throw UnsupportedError(ErrorKind{"Structure"},
                    "The distribution of 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) + ", which is not Markovian: method '" +
                    method +
                    "' reads every arc as a (D0,D1) process to build the holding-time laws. Use "
                    "method 'smp', which solves the semi-Markov environment, or fit the arc to a "
                    "phase-type law, e.g. APH.fitMeanAndSCV.");
            }
        }
    }
 
    const EnvStage<T>& stage(std::size_t e) const {
        check(e);
        return stages_[e];
    }
    const EnvArc<T>& arc(std::size_t e, std::size_t h) const {
        check(e);
        check(h);
        return arcs_[e][h];
    }
 
    /**
     * `resetEnvRatesFun{e,h}`: make the e -> h transition depend on the state
     * the stage is left in. Only `method = "statedep"` reads it.
     */
    void set_env_rate_reset(std::size_t e, std::size_t h, const ResetEnvRates<T>& f) {
        check(e);
        check(h);
        if (!arcs_[e][h].enabled)
            throw InputError("Environment: no transition " + stages_[e].name + " -> " +
                             stages_[h].name +
                             " is declared, so its rate cannot be made state dependent");
        arcs_[e][h].reset_rates = f;
    }
 
    /**
     * Replace the distribution of an arc that is already declared.
     *
     * This is what the state-dependent method writes back each iteration; every
     * other caller declares the arc once through `add_transition`. `init()` must
     * be rerun afterwards, because the superposed holding times and the stage
     * probabilities are all derived from these distributions.
     */
    void set_transition_dist(std::size_t e, std::size_t h, const lang::Distrib<T>& d) {
        check(e);
        check(h);
        if (!arcs_[e][h].enabled)
            throw InputError("Environment: no transition " + stages_[e].name + " -> " +
                             stages_[h].name + " is declared, so its distribution cannot be set");
        arcs_[e][h].dist = d;
    }
 
    /** The index of the stage called `nm`, or `nstages()` when there is none. */
    std::size_t find_stage(const std::string& nm) const {
        for (std::size_t e = 0; e < stages_.size(); ++e)
            if (stages_[e].has_model && stages_[e].name == nm) return e;
        return stages_.size();
    }
 
    /**
     * Install a reset policy on an arc that is already declared.
     *
     * The arc must exist: a reset on a disabled arc is a policy for a switch
     * the environment cannot make, and it would sit there reporting nothing.
     * The reference reaches this through `setBreakdownResetPolicy`, which
     * likewise errors when the stages it names are absent.
     */
    void set_reset(std::size_t e, std::size_t h, const ResetMarginal& reset) {
        check(e);
        check(h);
        if (!arcs_[e][h].enabled)
            throw InputError("Environment: no transition " + stages_[e].name + " -> " +
                             stages_[h].name +
                             " is declared, so it cannot be given a reset policy");
        arcs_[e][h].reset = reset;
    }
 
    // ---- node breakdown and repair ---------------------------------------
 
    /** `nodeFailures`, the declarative record of the breakdowns declared here. */
    const std::vector<NodeFailure<T>>& node_failures() const { return node_failures_; }
 
    /** `findNodeFailure`: the descriptor for `nm`, or `node_failures().size()`. */
    std::size_t find_node_failure(const std::string& nm) const {
        for (std::size_t i = 0; i < node_failures_.size(); ++i)
            if (node_failures_[i].node == nm) return i;
        return node_failures_.size();
    }
 
    /** The name `addNodeBreakdown` gives the stage in which `nm` is down. */
    static std::string down_stage_name(const std::string& nm) { return "DOWN_" + nm; }
 
    /**
     * Port of `addNodeBreakdown`, on a FIXED stage count.
     *
     * The reference grows its stage graph as breakdowns are declared; this
     * environment is sized at construction, exactly as `set_stage` is, so the
     * caller says which slot is UP and which is the DOWN stage of this node.
     * Everything else is the reference's: the UP stage holds the base model and
     * is named `UP`, the DOWN stage holds the degraded copy and is named
     * `DOWN_<node>`, and the UP -> DOWN arc carries the breakdown time and the
     * breakdown reset policy.
     */
    void add_node_breakdown(std::size_t up, std::size_t down, const qn::NetworkStruct<T>& base,
                            const std::string& node_name, const lang::Distrib<T>& breakdown,
                            const lang::Distrib<T>& down_service,
                            const std::string& reset_policy = "keep") {
        check(up);
        check(down);
        if (up == down)
            throw InputError("Environment: the UP and DOWN stages of node '" + node_name +
                             "' must be different stages");
        if (!stages_[up].has_model) set_stage(up, "UP", "operational", base);
        set_stage(down, down_stage_name(node_name), "failed",
                  env_degraded_model(base, node_name, down_service));
        add_transition(up, down, breakdown, env_reset_policy(reset_policy));
 
        NodeFailure<T> nf;
        nf.node = node_name;
        nf.breakdown = breakdown;
        nf.down_service = down_service;
        nf.breakdown_reset = reset_policy;
        record_node_failure(nf);
    }
 
    /** Port of `addNodeRepair`: the DOWN_<node> -> UP arc and its policy. */
    void add_node_repair(const std::string& node_name, const lang::Distrib<T>& repair,
                         const std::string& reset_policy = "keep") {
        const std::size_t down = find_stage(down_stage_name(node_name));
        const std::size_t up = find_stage("UP");
        if (down == stages_.size())
            throw InputError("Environment: stage '" + down_stage_name(node_name) +
                             "' is not defined; declare the breakdown of node '" + node_name +
                             "' before its repair");
        if (up == stages_.size())
            throw InputError("Environment: no UP stage is defined, so node '" + node_name +
                             "' has nothing to be repaired into");
        add_transition(down, up, repair, env_reset_policy(reset_policy));
 
        const std::size_t idx = find_node_failure(node_name);
        if (idx == node_failures_.size())
            throw InputError("Environment: no breakdown is recorded for node '" + node_name +
                             "', so its repair would describe a failure that was never declared");
        node_failures_[idx].repair = repair;
        node_failures_[idx].has_repair = true;
        node_failures_[idx].repair_reset = reset_policy;
    }
 
    /** `addNodeFailureRepair`: the two calls above, in order. */
    void add_node_failure_repair(std::size_t up, std::size_t down,
                                 const qn::NetworkStruct<T>& base, const std::string& node_name,
                                 const lang::Distrib<T>& breakdown, const lang::Distrib<T>& repair,
                                 const lang::Distrib<T>& down_service,
                                 const std::string& breakdown_reset = "keep",
                                 const std::string& repair_reset = "keep") {
        add_node_breakdown(up, down, base, node_name, breakdown, down_service, breakdown_reset);
        add_node_repair(node_name, repair, repair_reset);
    }
 
    /**
     * Port of `registerNodeFailure`: attach a breakdown descriptor, and its
     * reset policies, to stages that ALREADY exist.
     *
     * This is the read path of an environment saved in its expanded form: the
     * `UP` and `DOWN_<node>` stages and their two arcs came off the wire, and
     * the only thing the wire could not carry is the pair of policies, which is
     * what this installs.
     */
    void register_node_failure(const std::string& node_name, const lang::Distrib<T>& breakdown,
                               const lang::Distrib<T>& down_service, bool has_repair,
                               const lang::Distrib<T>& repair,
                               const std::string& breakdown_reset,
                               const std::string& repair_reset) {
        const std::size_t up = find_stage("UP");
        const std::size_t down = find_stage(down_stage_name(node_name));
        if (up == stages_.size())
            throw InputError("Environment: cannot register a node failure on '" + node_name +
                             "': no UP stage is defined in this environment");
        if (down == stages_.size())
            throw InputError("Environment: cannot register a node failure on '" + node_name +
                             "': no '" + down_stage_name(node_name) +
                             "' stage is defined in this environment");
        set_reset(up, down, env_reset_policy(breakdown_reset));
 
        NodeFailure<T> nf;
        nf.node = node_name;
        nf.breakdown = breakdown;
        nf.down_service = down_service;
        nf.breakdown_reset = breakdown_reset;
        if (has_repair) {
            set_reset(down, up, env_reset_policy(repair_reset));
            nf.repair = repair;
            nf.has_repair = true;
            nf.repair_reset = repair_reset;
        }
        record_node_failure(nf);
    }
 
    /**
     * Port of `Environment.init()`.
     *
     * Superpose the outgoing transitions of each stage, read the embedded jump
     * chain off the per-destination rates, and solve the resulting semi-Markov
     * process for its stationary stage probabilities. A semi-Markov stage takes
     * its jump row and mean holding time from `mc::smp_stage_quad` instead, and
     * keeps its arcs so `hold_time_cdf` can evaluate the law.
     *
     * An ABSORBING stage (`is_absorbing_stage`, no outgoing arc) gets an
     * infinite mean holding time, a zero embedded row and rate 0. The
     * environment then has no single steady state, so probEnv and probOrig are
     * left EMPTY; method 'occupation' reads the embedded chain instead.
     */
    void init() {
        const std::size_t E = stages_.size();
        require_declared();
 
        // proc[e][h], the transition process, as an MMAP marked by destination.
        // A semi-Markov stage has no MAP form, so its arcs are the disabled pair.
        sm_stage_.assign(E, 0);
        sm_arcs_.assign(E, std::vector<lang::Distrib<T>>());
        sm_probs_.assign(E, std::vector<double>());
        for (std::size_t e = 0; e < E; ++e) sm_stage_[e] = is_semi_markov_stage(e) ? 1 : 0;
        proc.assign(E, std::vector<mam::Mmap<double>>(E));
        for (std::size_t e = 0; e < E; ++e)
            for (std::size_t h = 0; h < E; ++h)
                proc[e][h] = sm_stage_[e] ? disabled_mmap(E) : arc_mmap(e, h, E);
 
        hold_time.assign(E, mam::Mmap<double>());
        hold_mean.assign(E, 0.0);
        Matrix<double> Pemb(E, E, 0.0);
        std::vector<double> lambda(E, 0.0);
        for (std::size_t e = 0; e < E; ++e) {
            if (is_absorbing_stage(e)) {
                // No exit: the stage is never left, so it has no jump law and an
                // infinite mean holding time (lambda[e] = 0, its Pemb row zero).
                hold_time[e] = disabled_mmap(E);
                hold_mean[e] = std::numeric_limits<double>::infinity();
                continue;
            }
            if (sm_stage_[e]) {
                // The embedded jump probabilities and the mean holding time come
                // from quadrature on the arc CDFs.
                stage_arcs(e, sm_arcs_[e], sm_probs_[e]);
                const mc::SmpStageQuad q = mc::smp_stage_quad(sm_arcs_[e], sm_probs_[e]);
                for (std::size_t h = 0; h < E; ++h) Pemb(e, h) = q.P[h];
                hold_mean[e] = q.m;
                hold_time[e] = disabled_mmap(E);
                lambda[e] = (q.m > 0.0) ? 1.0 / q.m : 0.0;
                continue;
            }
            // The reference seeds the superposition with the SELF transition
            // and folds in every other destination.
            mam::Mmap<double> ht = proc[e][e];
            for (std::size_t h = 0; h < E; ++h) {
                if (h == e) continue;
                ht = superpose_collapsed(ht, proc[e][h]);
            }
            hold_time[e] = ht;
            const std::vector<double> cl = mam::mmap_count_lambda(ht);
            double tot = 0.0;
            for (double v : cl) tot += v;
            if (tot > 0.0)
                for (std::size_t h = 0; h < E; ++h) Pemb(e, h) = cl[h] / tot;
            const double m = mam::map_mean(ht.map());
            hold_mean[e] = m;
            lambda[e] = (m > 0.0) ? 1.0 / m : 0.0;
        }
 
        pemb = Pemb;
        rate = lambda;
        bool all_positive = true;
        for (double v : lambda)
            if (!(v > 0.0)) all_positive = false;
        if (!all_positive) {
            // An absorbing environment has no single steady state: its limit
            // depends on the initial stage, so it is left to the solver
            // (SolverENV method 'occupation').
            prob_env.clear();
            prob_orig = Matrix<double>();
            return;
        }
 
        // A = -lambda_e (I - Pemb): the generator of the semi-Markov process
        // observed at its jumps, whose stationary law is the stage probability.
        Matrix<double> A(E, E, 0.0);
        for (std::size_t e = 0; e < E; ++e)
            for (std::size_t h = 0; h < E; ++h)
                A(e, h) = -lambda[e] * ((e == h ? 1.0 : 0.0) - Pemb(e, h));
        prob_env = mc::ctmc_solve_reducible(A).pi;
 
        prob_orig = Matrix<double>(E, E, 0.0);
        for (std::size_t e = 0; e < E; ++e) {
            for (std::size_t h = 0; h < E; ++h)
                prob_orig(h, e) = prob_env[h] * lambda[h] * Pemb(h, e);
            if (prob_env[e] > 0.0) {
                double s = 0.0;
                for (std::size_t h = 0; h < E; ++h) s += prob_orig(h, e);
                if (s > 0.0)
                    for (std::size_t h = 0; h < E; ++h) prob_orig(h, e) /= s;
            }
        }
    }
 
    /**
     * `getReliabilityTable`: MTTF, MTTR, MTBF and availability of an
     * environment built out of node breakdowns.
     *
     * It reads the UP and DOWN_* stages BY NAME, as the reference does, so it
     * applies to a breakdown environment and to nothing else; a stage set that
     * carries no such names is refused rather than answered about.
     *
     * The four are not independent readings of the same thing. MTTF is the
     * COMPETING-RISK mean of the arcs out of UP, i.e. one over the SUM of the
     * breakdown rates, so a second failing node shortens it. MTTR is the mean
     * repair time weighted by the CONDITIONAL probability of being in each DOWN
     * stage given that the system is down. Availability is read off probEnv and
     * not from MTTF/(MTTF+MTTR): the two agree for a Markovian environment and
     * the stationary law is the one that stays right when the arcs are not.
     */
    struct Reliability {
        double mttf = 0.0;          ///< mean time to failure, UP -> any DOWN
        double mttr = 0.0;          ///< mean time to repair, DOWN -> UP
        double mtbf = 0.0;          ///< MTTF + MTTR
        double availability = 0.0;  ///< stationary probability of being UP
    };
 
    Reliability reliability() const {
        if (prob_env.empty())
            throw InputError(
                "Environment: reliability metrics are read from probEnv; call init() first");
        const std::size_t E = stages_.size();
        const std::size_t up = find_stage("UP");
        if (up == E)
            throw InputError(
                "Environment: no UP stage, so there is nothing for a breakdown to leave; "
                "reliability metrics apply to an environment built with add_node_breakdown");
        std::vector<std::size_t> down;
        for (std::size_t e = 0; e < E; ++e)
            if (stages_[e].name.compare(0, 5, "DOWN_") == 0) down.push_back(e);
        if (down.empty())
            throw InputError(
                "Environment: no DOWN_<node> stage, so no node breaks down; reliability metrics "
                "apply to an environment built with add_node_breakdown");
 
        double lambda_total = 0.0;
        for (std::size_t h : down)
            if (arcs_[up][h].enabled) lambda_total += arc_rate(up, h);
        if (!(lambda_total > 0.0))
            throw InputError(
                "Environment: no UP -> DOWN transition, so the system never fails and its mean "
                "time to failure is not a number");
 
        std::vector<double> mu, p;
        for (std::size_t e : down)
            if (arcs_[e][up].enabled) {
                mu.push_back(is_semi_markov_stage(e) ? 1.0 / hold_mean[e] : arc_rate(e, up));
                p.push_back(prob_env[e]);
            }
        if (mu.empty())
            throw InputError(
                "Environment: no DOWN -> UP transition, so the system is never repaired and its "
                "mean time to repair is not a number");
 
        Reliability r;
        r.mttf = 1.0 / lambda_total;
        // A semi-Markov UP stage has no rates to add: the time to failure is its holding time.
        if (is_semi_markov_stage(up)) r.mttf = hold_mean[up];
        double ptot = 0.0;
        for (double v : p) ptot += v;
        r.mttr = 0.0;
        for (std::size_t i = 0; i < mu.size(); ++i)
            r.mttr += (ptot > 0.0 ? p[i] / ptot : 1.0 / static_cast<double>(mu.size())) / mu[i];
        r.mtbf = r.mttf + r.mttr;
        double pup = prob_env[up], pdown = 0.0;
        for (std::size_t e : down) pdown += prob_env[e];
        r.availability = (pup + pdown > 0.0) ? pup / (pup + pdown) : 0.0;
        return r;
    }
 
    /**
     * `getStageTable`: one row per stage. The reference's HoldT column holds the
     * MAP or the semi-Markov arc struct; here it is the mean holding time and
     * whether the stage is semi-Markov, which is what a reader of the table uses.
     */
    struct StageRow {
        std::size_t stage = 0;  ///< 1-based, as the reference prints it
        std::string name;
        std::string type;
        double prob = 0.0;       ///< probEnv
        double hold_mean = 0.0;  ///< holdMean
        bool semi_markov = false;
    };
 
    std::vector<StageRow> stage_table() const {
        if (prob_env.empty())
            throw InputError("Environment: the stage table reads probEnv; call init() first");
        std::vector<StageRow> rows(stages_.size());
        for (std::size_t e = 0; e < stages_.size(); ++e) {
            rows[e].stage = e + 1;
            rows[e].name = stages_[e].name;
            rows[e].type = stages_[e].type;
            rows[e].prob = prob_env[e];
            rows[e].hold_mean = hold_mean[e];
            rows[e].semi_markov = hold_time_is_semi_markov(e);
        }
        return rows;
    }
 
    /**
     * `printStageTable`, as text: the stages with their networks, then every
     * enabled arc in (from, to) order, a competing one by its rate and a kernel
     * one by its probability, sojourn family and mean. A layered stage reports
     * its task and entry counts, having no node or class count of its own.
     */
    std::string stage_table_text() const {
        std::string s = "Stage Table:\n============\n";
        char buf[512];
        for (std::size_t e = 0; e < stages_.size(); ++e) {
            const EnvStage<T>& st = stages_[e];
            std::snprintf(buf, sizeof buf, "Stage %zu: %s (Type: %s)\n", e + 1, st.name.c_str(), st.type.c_str());
            s += buf;
            if (st.has_model) {
                std::snprintf(buf, sizeof buf, "  - Network: %s\n  - Nodes: %zu\n  - Classes: %zu\n",
                              st.model.name.c_str(), st.model.nodes.size(), st.model.nclasses);
                s += buf;
            } else if (st.has_lqn) {
                std::snprintf(buf, sizeof buf, "  - Layered network: %zu tasks, %zu entries\n",
                              st.lqn_model.ntasks, st.lqn_model.nentries);
                s += buf;
            }
        }
        std::string tr;
        for (std::size_t e = 0; e < stages_.size(); ++e)
            for (std::size_t h = 0; h < stages_.size(); ++h) {
                const EnvArc<T>& a = arcs_[e][h];
                if (!a.enabled || a.dist.disabled) continue;
                const double m = num_traits<T>::to_double(a.dist.mean);
                if (is_kernel_stage(e))
                    std::snprintf(buf, sizeof buf, "  %s -> %s: prob = %.4f, sojourn %s, mean = %.4f\n",
                                  stages_[e].name.c_str(), stages_[h].name.c_str(), kernel_prob_(e, h),
                                  lang::process_to_text(a.dist.type), m);
                else
                    std::snprintf(buf, sizeof buf, "  %s -> %s: rate = %.4f\n", stages_[e].name.c_str(),
                                  stages_[h].name.c_str(), 1.0 / m);
                tr += buf;
            }
        if (!tr.empty()) s += "\nTransitions:\n" + tr;
        return s;
    }
 
    /** `printStageTable`: writes `stage_table_text()` to `out`. */
    void print_stage_table(std::FILE* out = stdout) const {
        std::fputs(stage_table_text().c_str(), out);
    }
 
    /**
     * `Environment.holdTimeCdf(holdTime{e}, t)`: the holding-time CDF of stage e
     * at the instants `t`, from the MAP for a Markovian stage and from the arc
     * race or kernel mixture (`mc::smp_stage_cif`) for a semi-Markov one.
     */
    std::vector<double> hold_time_cdf(std::size_t e, const std::vector<double>& t) const {
        check(e);
        if (e >= sm_stage_.size())
            throw InputError("Environment: the holding-time law is built by init(); call it first");
        if (sm_stage_[e]) return mc::smp_stage_cif(sm_arcs_[e], sm_probs_[e], t).Fhold;
        return mam::map_cdf(hold_time[e].map(), t);
    }
 
    /** `Environment.holdTimeMean(holdTime{e})`: the mean holding time of stage e. */
    double hold_time_mean(std::size_t e) const {
        check(e);
        if (e >= hold_mean.size())
            throw InputError("Environment: the holding-time law is built by init(); call it first");
        return hold_mean[e];
    }
 
    /** True when init() built stage e as a semi-Markov stage (no MAP holding time). */
    bool hold_time_is_semi_markov(std::size_t e) const {
        check(e);
        return e < sm_stage_.size() && sm_stage_[e];
    }
 
    // ---- what init() produced --------------------------------------------
    std::vector<std::vector<mam::Mmap<double>>> proc;  ///< proc[e][h]; disabled for a semi-Markov stage
    std::vector<mam::Mmap<double>> hold_time;          ///< holdTime[e]; disabled for a semi-Markov stage
    std::vector<double> hold_mean;                     ///< holdMean
    std::vector<double> prob_env;                      ///< probEnv
    Matrix<double> prob_orig;                          ///< probOrig(h, e)
    Matrix<double> pemb;                               ///< probEmb, the embedded jump chain
    std::vector<double> rate;                          ///< 1/E[holding time]
 
private:
    /** `isa(arc, 'Markovian')`: a family with a (D0,D1) form that carries one. */
    static bool arc_is_markovian(const lang::Distrib<T>& d) {
        return lang::process_is_markovian(d.type) && d.has_map();
    }
 
    /** The 1 x 1 zero pair a disabled arc is. */
    static mam::Mmap<double> disabled_mmap(std::size_t E) {
        mam::Mmap<double> m;
        m.D0 = Matrix<double>(1, 1, 0.0);
        m.D1 = Matrix<double>(1, 1, 0.0);
        m.Dc.assign(E, Matrix<double>(1, 1, 0.0));
        return m;
    }
 
    void check(std::size_t e) const {
        if (e >= stages_.size()) throw InputError("Environment: stage index out of range");
    }
 
    /** `1 / env{e,h}.getMean()`: the arc rate the reliability reading uses. */
    double arc_rate(std::size_t e, std::size_t h) const {
        const double m = num_traits<T>::to_double(arcs_[e][h].dist.mean);
        if (!(m > 0.0))
            throw InputError("Environment: the transition " + stages_[e].name + " -> " +
                             stages_[h].name +
                             " has no positive mean, so it carries no rate to report");
        return 1.0 / m;
    }
 
    /** Record a descriptor, replacing the one this node already had. */
    void record_node_failure(const NodeFailure<T>& nf) {
        const std::size_t idx = find_node_failure(nf.node);
        if (idx < node_failures_.size())
            node_failures_[idx] = nf;
        else
            node_failures_.push_back(nf);
    }
 
    /**
     * `emmap{e}{h}`: the arc's process marked by destination -- its D1 in slot
     * h and zero in every other slot. A disabled arc is the 1 x 1 zero pair.
     */
    mam::Mmap<double> arc_mmap(std::size_t e, std::size_t h, std::size_t E) const {
        mam::Mmap<double> m;
        if (!arcs_[e][h].enabled) {
            m.D0 = Matrix<double>(1, 1, 0.0);
            m.D1 = Matrix<double>(1, 1, 0.0);
            m.Dc.assign(E, Matrix<double>(1, 1, 0.0));
            return m;
        }
        const lang::Distrib<T>& d = arcs_[e][h].dist;
        const std::size_t n = d.D0.rows();
        if (n == 0)
            throw InputError("Environment: the transition distribution has no representation");
        m.D0 = Matrix<double>(n, n, 0.0);
        m.D1 = Matrix<double>(n, n, 0.0);
        for (std::size_t a = 0; a < n; ++a)
            for (std::size_t b = 0; b < n; ++b) {
                m.D0(a, b) = num_traits<T>::to_double(d.D0(a, b));
                m.D1(a, b) = num_traits<T>::to_double(d.D1(a, b));
            }
        m.Dc.assign(E, Matrix<double>(n, n, 0.0));
        m.Dc[h] = m.D1;
        return m;
    }
 
    /**
     * One fold of the reference's superposition loop: Kronecker-sum the blocks,
     * then move each block's row sums into its FIRST column.
     *
     * The collapse is what makes the superposed holding time RENEW at every
     * jump: without it the phase surviving from the losing risks would carry
     * into the next stage, which is not the semi-Markov model the solver then
     * solves.
     */
    static mam::Mmap<double> superpose_collapsed(const mam::Mmap<double>& a,
                                                 const mam::Mmap<double>& b) {
        mam::Mmap<double> s;
        s.D0 = mam::krons(a.D0, b.D0);
        s.D1 = collapse(mam::krons(a.D1, b.D1));
        s.Dc.reserve(a.Dc.size());
        for (std::size_t c = 0; c < a.Dc.size(); ++c)
            s.Dc.push_back(collapse(mam::krons(a.Dc[c], b.Dc[c])));
        return mam::mmap_normalize(s);
    }
 
    /** Row sums into column one, the rest zeroed. */
    static Matrix<double> collapse(const Matrix<double>& X) {
        Matrix<double> Y(X.rows(), X.cols(), 0.0);
        for (std::size_t i = 0; i < X.rows(); ++i) {
            double s = 0.0;
            for (std::size_t j = 0; j < X.cols(); ++j) s += X(i, j);
            Y(i, 0) = s;
        }
        return Y;
    }
 
    std::string name_;
    std::vector<EnvStage<T>> stages_;
    std::vector<std::vector<EnvArc<T>>> arcs_;
    std::vector<NodeFailure<T>> node_failures_;
    Matrix<double> kernel_prob_;  ///< kernelProb, NaN for a competing arc or none
    std::vector<char> sm_stage_;  ///< per stage, whether init() built it semi-Markov
    std::vector<std::vector<lang::Distrib<T>>> sm_arcs_;  ///< a semi-Markov stage's arcs
    std::vector<std::vector<double>> sm_probs_;           ///< and its kernel probabilities
};
 
}  // namespace env
}  // namespace line
 
#endif  // LINE_LANG_QN_ENVIRONMENT_H