fluid_petri_terms.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_SOLVERS_FLUID_PETRI_TERMS_H
#define LINE_SOLVERS_FLUID_PETRI_TERMS_H
 
/**
 * @file
 * @ingroup line_solvers
 * Event-based representation of the fluid marking process of a stochastic Petri
 * net. Port of `matlab/src/solvers/FLD/fluid_petri_terms.m`, cross-checked
 * against `jar/src/main/java/jline/solvers/fluid/petri/PetriTerms.java`.
 *
 * A GSPN is already a density-dependent Markov population process, which is the
 * object the moment-closure family of SolverFLD is built on: the marking is the
 * population, a transition mode is a reaction, its incidence column is the jump,
 * and the rate law `lambda*min(enabling degree, servers)` is the same min()
 * non-linearity the min-normal closure exists to smooth. Nothing about the
 * closure changes here; only where the drift comes from.
 *
 *     dx/dt = D * r(x, Sigma, phi, mu)
 *
 * THE STATE, x = [ m ; y ].
 *
 *   m(p,k)  token mass of class k at place p. One coordinate per (place, class)
 *           pair some arc touches, the initial marking loads, or a Source feeds;
 *           a pair nothing reaches is dropped rather than carried as a null
 *           direction of the Newton system.
 *   y(j,h)  the number of mode-j servers running in phase h, for a mode whose
 *           firing time has more than one phase. Their SUM is not free: the
 *           ENABLE synchronization latches it instantaneously to
 *           min(enabling degree, servers), so the latch is an ALGEBRAIC row with
 *           one free-sign unknown mu_j and the phase split evolves
 *           differentially. CARRYING THE DISTRIBUTION INSTEAD OF THE COUNT LOOKS
 *           TIDIER AND IS WRONG: the resulting equation is missing a term and
 *           agrees with the count form only AT a fixed point.
 *   f(c)    the LEVEL of fluid place c (FSPN), one coordinate per fluid place, appended after the phase block.
 *           A fluid place holds one real level and no tokens; it moves only by the flows of the modes that run
 *           and by the impulses of the modes that fire.
 *
 * There is no Source coordinate: an exogenous arrival is a CONSTANT-propensity
 * event depositing one token, as it is in the NRM SPN runner. There is no Sink
 * coordinate either: a firing arc into a sink is mass leaving the net.
 *
 * THE EVENTS, one column of D each: kind 1 a firing of mode j out of phase h
 * into phase h'; kind 2 an internal phase change; kind 3 an exogenous arrival;
 * kind 4 a firing of an IMMEDIATE mode, at the algebraic flow phi_j; kind 5 the
 * server latch of a multi-phase mode, at the free-sign unknown mu_j; kind 6 the FLUID FLOW of mode j into
 * (or out of) fluid place c, a unit column on f(c) signed by the direction of the flow, at |r_jc| * theta_j(x).
 * A flow is DETERMINISTIC: it enters the drift and its Jacobian but not the diffusion, as a server latch and an
 * immediate flow do not. An impulse is not a column of its own -- it is the fluid entry of the mode's firing
 * column, so it is a stochastic jump.
 *
 * Arithmetic: DOUBLE ONLY. The closures evaluate the normal CDF, which is not an
 * element of the field generated by the inputs under any arithmetic, so the
 * templated form this tree uses elsewhere would buy nothing.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <functional>
#include <limits>
#include <map>
#include <string>
#include <vector>
 
#include "line/lang/qn/network_struct.h"
#include "line/util/error.h"
#include "line/util/matrix.h"
 
namespace line {
namespace fluid {
namespace petri {
 
/** GlobalConstants.FineTol, the reference's own "effectively zero". */
inline double petri_fine_tol() { return 1e-8; }
 
/** GlobalConstants.CoarseTol, also the relative width of the seed's smoothed inhibitor gate. */
inline double petri_coarse_tol() { return 1e-3; }
 
/** One firing mode of one transition, with its arcs and its firing process. */
struct PetriMode {
    std::size_t node = 0;   ///< 1-based node index of the transition
    std::size_t mode = 0;   ///< mode index within it, 0-based
    lang::TimingStrategy timing = lang::TimingStrategy::TIMED;
    std::string label;
 
    /** Input arcs, as state coordinates and their multiplicities. */
    std::vector<std::size_t> arc_slot;
    std::vector<double> arc_w;
    /** Inhibitor arcs, as state coordinates and their thresholds. */
    std::vector<std::size_t> inh_slot;
    std::vector<double> inh_thr;
 
    double c = 1.0;      ///< servers of this mode; infinite for none
    std::size_t nph = 1; ///< phases of the firing process; 0 for an immediate mode
    Matrix<double> D0, D1;
    std::vector<double> d1;   ///< row sums of D1, the completion rate out of each phase
    std::vector<double> pie;  ///< entry distribution over the phases
    /** Marking-dependent firing multiplier; empty for the unit one. */
    std::function<double(const Matrix<double>&)> dep;
    int prio = 1;
    double weight = 1.0;
    /** The mode's incidence column, over the whole state. */
    std::vector<double> cvec;
    /** The phase coordinates of a multi-phase mode; empty otherwise. */
    std::vector<std::size_t> zblk;
    /**
     * Whether the enabling degree is worth closing over.
     *
     * A single input arc under an infinite server count makes min() the identity,
     * so the closure would smooth nothing and the covariance entries it reads
     * would be carried for no reason.
     */
    bool closable = false;
 
    /** Fluid flows of the mode: level coordinates and signed rates per running server. */
    std::vector<std::size_t> flow_slot;
    std::vector<double> flow_rate;
    /** Fluid guards: level coordinates, +1 for '>=' and -1 for '<=', and thresholds. */
    std::vector<std::size_t> guard_slot;
    std::vector<int> guard_sgn;
    std::vector<double> guard_thr;
    /** Fluid rate dependences: level coordinates and the multipliers h(level). */
    std::vector<std::size_t> fdep_slot;
    std::vector<std::function<double(double)>> fdep_fun;
    /**
     * Marking-dependent flow multipliers, aligned with `flow_slot` (an empty function is the unit one), each of the
     * (I x K) marking restricted to the 1-based token places of `flow_dep_read` (empty: every coordinate).
     */
    std::vector<std::function<double(const Matrix<double>&)>> flow_dep_fun;
    std::vector<std::vector<std::size_t>> flow_dep_read;
};
 
/** The assembled drift terms of a net. */
struct PetriTerms {
    std::size_t I = 0;  ///< nodes
    std::size_t K = 0;  ///< classes
    std::size_t M = 0;  ///< stations
    std::vector<std::size_t> places, transitions;  ///< 1-based node indices
    std::vector<std::string> names_node;           ///< 0-based, one per node
 
    std::size_t nstate = 0;  ///< marking coordinates plus phase coordinates plus fluid levels
    std::size_t nm = 0;      ///< marking coordinates alone
 
    /** Fluid places (FSPN): 1-based nodes, their level coordinates, capacities and initial levels. */
    std::vector<std::size_t> fluid_places;
    std::size_t nf = 0;
    std::vector<std::ptrdiff_t> fidx;  ///< (I) level coordinate of a 0-based node, or -1
    std::vector<std::size_t> fluid_slot;
    std::vector<double> fluid_cap, fluid_init;
    std::vector<std::ptrdiff_t> fluid_station;  ///< 0-based station of each fluid place
    std::vector<std::string> fluid_name;
 
    /** pidx(p,k) is the state coordinate of (0-based node p, class k), or npos. */
    std::vector<std::vector<std::ptrdiff_t>> pidx;
    std::vector<std::size_t> coord_node, coord_class;
    std::vector<std::ptrdiff_t> coord_station;
 
    std::vector<PetriMode> modes;
    std::vector<std::size_t> timed_idx, imm_idx;
 
    Matrix<double> D;              ///< (nstate x nev) incidence
    std::vector<double> rate_base;
    std::size_t nev = 0;
    std::vector<int> ev_kind, ev_mode, ev_phase, ev_to;
    std::vector<std::ptrdiff_t> ev_station, ev_class;
    /** Kind 6: the flow's position in its mode's `flow_slot`; 0 for every other kind. */
    std::vector<std::size_t> ev_flow;
 
    std::vector<std::size_t> imm_col, latch_col, stoch_col, flow_col;
    std::vector<std::size_t> latch_mode;
 
    /** The Sigma entries the closure reads, and the map back to them. */
    std::vector<std::pair<std::size_t, std::size_t>> cov_pairs;
    std::size_t npair = 0;
    std::vector<std::vector<std::ptrdiff_t>> pair_index;
    std::vector<std::size_t> cov_idx;
 
    /** (station*K + class) -> the events that take tokens out of / into it. */
    std::map<std::size_t, std::vector<std::size_t>> consumers, producers;
    std::map<std::size_t, std::vector<double>> consumer_w;
 
    std::vector<double> x0;
    std::vector<std::vector<double>> m0full;
};
 
namespace terms_detail {
 
/** An (I x K) parameter matrix padded to full size with `fill`. */
template <class T>
inline std::vector<std::vector<double>> pad(const Matrix<T>& A, std::size_t I, std::size_t K,
                                            double fill) {
    std::vector<std::vector<double>> out(I, std::vector<double>(K, fill));
    for (std::size_t i = 0; i < I && i < A.rows(); ++i)
        for (std::size_t k = 0; k < K && k < A.cols(); ++k)
            out[i][k] = num_traits<T>::to_double(A(i, k));
    return out;
}
 
/** Record a covariance pair once, in ascending order. */
inline void add_pair(std::map<std::pair<std::size_t, std::size_t>, std::size_t>& key,
                     std::vector<std::pair<std::size_t, std::size_t>>& pairs, std::size_t a,
                     std::size_t b) {
    const std::pair<std::size_t, std::size_t> k(std::min(a, b), std::max(a, b));
    if (key.find(k) != key.end()) return;
    key[k] = pairs.size();
    pairs.push_back(k);
}
 
/**
 * One mode's arcs, incidence column and firing process.
 *
 * `nm` is the marking width at the time of the call; the column is GROWN once
 * the phase coordinates are appended and the total state width is known.
 */
template <class T>
inline PetriMode build_mode(const qn::NetworkStruct<T>& sn, const qn::TransitionParam<T>& tp,
                            std::size_t ind, std::size_t m, PetriTerms& t, std::size_t nm) {
    PetriMode rec;
    rec.node = ind;
    rec.mode = m;
    rec.timing = (m < tp.timing.size()) ? tp.timing[m] : lang::TimingStrategy::TIMED;
    const std::string mname = (m < tp.modenames.size() && !tp.modenames[m].empty())
                                  ? tp.modenames[m]
                                  : ("Mode" + std::to_string(m + 1));
    rec.label = t.names_node[ind - 1] + "." + mname;
 
    const std::vector<std::vector<double>> en = pad(tp.enabling[m], t.I, t.K, 0.0);
    const std::vector<std::vector<double>> fir = pad(tp.firing[m], t.I, t.K, 0.0);
    const std::vector<std::vector<double>> inh =
        pad(tp.inhibiting[m], t.I, t.K, std::numeric_limits<double>::infinity());
 
    std::vector<double> cvec(nm, 0.0);
    for (std::size_t p = 0; p < t.I; ++p)
        for (std::size_t k = 0; k < t.K; ++k) {
            const double w = en[p][k];
            if (!(w > 0)) continue;
            if (std::isinf(w))
                throw InputError("fluid_petri_terms: mode " + rec.label +
                                 " has a non-finite enabling arc weight at " + t.names_node[p] +
                                 ". An arc that no marking can satisfy disables the mode; declare "
                                 "a finite multiplicity");
            if (t.pidx[p][k] < 0)
                throw InputError("fluid_petri_terms: mode " + rec.label +
                                 " takes an enabling arc from " + t.names_node[p] +
                                 ", which is not a Place");
            const std::size_t s = static_cast<std::size_t>(t.pidx[p][k]);
            rec.arc_slot.push_back(s);
            rec.arc_w.push_back(w);
            cvec[s] -= w;
        }
    for (std::size_t p = 0; p < t.I; ++p)
        for (std::size_t k = 0; k < t.K; ++k) {
            // A NEGATIVE ENTRY MARKS AN INPUT PLACE, whose token the enabling arc
            // already removed, and a firing arc into a Sink has no coordinate:
            // that is mass leaving the net, not mass lost by accident.
            if (fir[p][k] <= 0 || t.pidx[p][k] < 0) continue;
            cvec[static_cast<std::size_t>(t.pidx[p][k])] += fir[p][k];
        }
    for (std::size_t p = 0; p < t.I; ++p)
        for (std::size_t k = 0; k < t.K; ++k) {
            const double th = inh[p][k];
            // JMT writes a MISSING inhibitor arc as 0, as -1 or as Inf, never as
            // a threshold, so none of those three is one.
            if (std::isinf(th) || th <= 0 || t.pidx[p][k] < 0) continue;
            rec.inh_slot.push_back(static_cast<std::size_t>(t.pidx[p][k]));
            rec.inh_thr.push_back(th);
        }
    rec.cvec = cvec;
 
    const double c = (m < tp.nmodeservers.size()) ? tp.nmodeservers[m] : 1.0;
    rec.c = std::isnan(c) ? 1.0 : c;
    rec.prio = (m < tp.firingprio.size()) ? static_cast<int>(tp.firingprio[m]) : 1;
    rec.weight =
        (m < tp.fireweight.size()) ? num_traits<T>::to_double(tp.fireweight[m]) : 1.0;
    if (m < tp.firingdep.size() && tp.firingdep[m]) {
        const std::function<T(const std::vector<T>&)> g = tp.firingdep[m];
        const std::size_t I = t.I, K = t.K;
        rec.dep = [g, I, K](const Matrix<double>& mm) -> double {
            std::vector<T> flat(I * K, num_traits<T>::from_int(0));
            for (std::size_t i = 0; i < I; ++i)
                for (std::size_t k = 0; k < K; ++k)
                    flat[i * K + k] = num_traits<T>::from_double(mm(i, k));
            return num_traits<T>::to_double(g(flat));
        };
    }
 
    if (rec.timing == lang::TimingStrategy::IMMEDIATE) {
        // An immediate mode carries no firing process: its flow is an algebraic
        // unknown of the DAE, not a rate.
        rec.nph = 0;
        rec.d1.assign(1, 0.0);
        rec.pie.assign(1, 1.0);
        rec.closable = false;
        return rec;
    }
 
    if (m >= tp.firingproc.size() || tp.firingproc[m].disabled)
        throw InputError("fluid_petri_terms: mode " + rec.label +
                         " has no Markovian firing process. The renewal families are converted "
                         "to phase type before the solver runs, so this is a distribution the "
                         "fluid Petri route cannot time; use SolverCTMC or SolverLDES");
    const lang::Distrib<T>& fp = tp.firingproc[m];
    Matrix<double> D0, D1;
    if (fp.D0.rows() > 0 && fp.D1.rows() > 0) {
        D0 = Matrix<double>(fp.D0.rows(), fp.D0.cols(), 0.0);
        for (std::size_t i = 0; i < D0.rows(); ++i)
            for (std::size_t j = 0; j < D0.cols(); ++j)
                D0(i, j) = num_traits<T>::to_double(fp.D0(i, j));
        D1 = Matrix<double>(fp.D1.rows(), fp.D1.cols(), 0.0);
        for (std::size_t i = 0; i < D1.rows(); ++i)
            for (std::size_t j = 0; j < D1.cols(); ++j)
                D1(i, j) = num_traits<T>::to_double(fp.D1(i, j));
    } else {
        // An exponential law carried by its mean alone: the one-phase (D0, D1).
        const double mean = num_traits<T>::to_double(fp.mean);
        if (!(mean > 0.0))
            throw InputError("fluid_petri_terms: mode " + rec.label +
                             " has a non-positive mean firing time");
        D0 = Matrix<double>(1, 1, -1.0 / mean);
        D1 = Matrix<double>(1, 1, 1.0 / mean);
    }
    rec.nph = D0.rows();
    rec.D0 = D0;
    rec.D1 = D1;
    rec.d1.assign(rec.nph, 0.0);
    for (std::size_t h = 0; h < rec.nph; ++h) {
        double acc = 0.0;
        for (std::size_t hp = 0; hp < D1.cols(); ++hp) acc += D1(h, hp);
        rec.d1[h] = acc;
    }
    rec.pie.assign(std::max<std::size_t>(rec.nph, 1), 0.0);
    if (rec.nph > 1) {
        double tot = 0.0;
        for (std::size_t h = 0; h < rec.nph && h < fp.params.size(); ++h) {
            rec.pie[h] = num_traits<T>::to_double(fp.params[h]);
            tot += rec.pie[h];
        }
        if (!(tot > 0.0)) {
            std::fill(rec.pie.begin(), rec.pie.end(), 0.0);
            rec.pie[0] = 1.0;
        } else {
            for (std::size_t h = 0; h < rec.nph; ++h) rec.pie[h] /= tot;
        }
    } else {
        rec.pie[0] = 1.0;
    }
    rec.closable = !(rec.arc_slot.size() <= 1 && std::isinf(rec.c));
    return rec;
}
 
}  // namespace terms_detail
 
/**
 * Assemble the drift terms of a net.
 *
 * @param sn a model whose nodes are Places, Transitions, Sources and Sinks only
 */
template <class T>
inline PetriTerms petri_build_terms(const qn::NetworkStruct<T>& sn) {
    PetriTerms t;
    t.I = sn.nodes.size();
    t.K = sn.nclasses;
    t.M = sn.nstations;
    t.names_node.resize(t.I);
    for (std::size_t i = 0; i < t.I; ++i) t.names_node[i] = sn.nodes[i].name;
    for (std::size_t i = 1; i <= t.I; ++i) {
        if (sn.nodes[i - 1].nodetype == lang::NodeType::Place) t.places.push_back(i);
        else if (sn.nodes[i - 1].nodetype == lang::NodeType::Transition) t.transitions.push_back(i);
    }
    // Fluid places carry a level, not a token count per class, so they take no marking coordinate; their
    // level coordinate is appended after the phase block.
    {
        std::vector<std::size_t> tok;
        for (std::size_t pi = 0; pi < t.places.size(); ++pi) {
            if (sn.is_fluid_node(t.places[pi])) t.fluid_places.push_back(t.places[pi]);
            else tok.push_back(t.places[pi]);
        }
        t.places = tok;
    }
 
    // ---- the initial marking ------------------------------------------------
    // A Place's DECLARED marking wins; absent one, each closed class loads its
    // reference station, which is the rule spn_sinvariants and the exact SPN
    // engines use. An infinite marking is refused rather than integrated.
    t.m0full.assign(t.I, std::vector<double>(t.K, 0.0));
    for (std::size_t pi = 0; pi < t.places.size(); ++pi) {
        const std::size_t ind = t.places[pi];
        const typename std::map<std::size_t, std::vector<T>>::const_iterator im =
            sn.initmarking.find(ind);
        if (im == sn.initmarking.end()) continue;
        for (std::size_t k = 0; k < t.K && k < im->second.size(); ++k) {
            const double v = num_traits<T>::to_double(im->second[k]);
            if (std::isinf(v))
                throw InputError("fluid_petri_terms: place " + t.names_node[ind - 1] +
                                 " holds an infinite initial marking of class " +
                                 std::to_string(k + 1));
            t.m0full[ind - 1][k] = v;
        }
    }
    bool any_declared = false;
    for (std::size_t pi = 0; pi < t.places.size() && !any_declared; ++pi)
        any_declared = (sn.initmarking.find(t.places[pi]) != sn.initmarking.end());
    if (!any_declared) {
        for (std::size_t r = 0; r < sn.classes.size(); ++r) {
            const double njobs = num_traits<T>::to_double(sn.classes[r].population);
            if (!(njobs > 0.0) || !std::isfinite(njobs)) continue;
            if (sn.classes[r].refstat < 1 || sn.classes[r].refstat > sn.station_to_node.size())
                continue;
            const std::size_t ref_node = sn.station_to_node[sn.classes[r].refstat - 1];
            if (ref_node >= 1 && ref_node <= t.I &&
                sn.nodes[ref_node - 1].nodetype == lang::NodeType::Place)
                t.m0full[ref_node - 1][r] += njobs;
        }
    }
 
    // ---- which (place, class) pairs carry a coordinate ----------------------
    std::vector<std::vector<bool>> touched(t.I, std::vector<bool>(t.K, false));
    for (std::size_t ti = 0; ti < t.transitions.size(); ++ti) {
        const std::size_t ind = t.transitions[ti];
        const typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
            sn.transparam.find(ind);
        if (it == sn.transparam.end()) continue;
        const qn::TransitionParam<T>& tp = it->second;
        for (std::size_t m = 0; m < tp.nmodes; ++m) {
            const std::vector<std::vector<double>> en =
                terms_detail::pad(tp.enabling[m], t.I, t.K, 0.0);
            const std::vector<std::vector<double>> fir =
                terms_detail::pad(tp.firing[m], t.I, t.K, 0.0);
            const std::vector<std::vector<double>> inh = terms_detail::pad(
                tp.inhibiting[m], t.I, t.K, std::numeric_limits<double>::infinity());
            for (std::size_t i = 0; i < t.I; ++i)
                for (std::size_t k = 0; k < t.K; ++k)
                    if (en[i][k] > 0 || fir[i][k] != 0 ||
                        (!std::isinf(inh[i][k]) && inh[i][k] > 0))
                        touched[i][k] = true;
        }
    }
 
    // An arrival makes its target a coordinate even when no arc mentions it.
    struct SrcArr {
        std::size_t snd, r, qnd, l;
    };
    std::vector<SrcArr> src_arr;
    for (std::size_t ind = 1; ind <= t.I; ++ind) {
        if (sn.nodes[ind - 1].nodetype != lang::NodeType::Source) continue;
        const std::size_t ist = sn.nodes[ind - 1].station;
        if (ist < 1) continue;
        for (std::size_t r = 0; r < t.K; ++r) {
            const double lambda = num_traits<T>::to_double(sn.rates(ist - 1, r));
            if (std::isnan(lambda) || lambda <= 0) continue;
            for (std::size_t pi = 0; pi < t.places.size(); ++pi) {
                const std::size_t jnd = t.places[pi];
                for (std::size_t s = 0; s < t.K; ++s) {
                    const double p = num_traits<T>::to_double(
                        sn.rtnodes((ind - 1) * t.K + r, (jnd - 1) * t.K + s));
                    if (p > 0) {
                        touched[jnd - 1][s] = true;
                        SrcArr a;
                        a.snd = ind;
                        a.r = r;
                        a.qnd = jnd;
                        a.l = s;
                        src_arr.push_back(a);
                    }
                }
            }
        }
    }
 
    t.pidx.assign(t.I, std::vector<std::ptrdiff_t>(t.K, -1));
    std::size_t nm = 0;
    for (std::size_t pi = 0; pi < t.places.size(); ++pi) {
        const std::size_t ind = t.places[pi];
        for (std::size_t k = 0; k < t.K; ++k) {
            if (!(touched[ind - 1][k] || t.m0full[ind - 1][k] > 0)) continue;
            t.pidx[ind - 1][k] = static_cast<std::ptrdiff_t>(nm);
            t.coord_node.push_back(ind - 1);
            t.coord_class.push_back(k);
            t.coord_station.push_back(static_cast<std::ptrdiff_t>(sn.nodes[ind - 1].station) - 1);
            ++nm;
        }
    }
    t.nm = nm;
 
    // ---- the modes, and the phase coordinates of the multi-phase ones -------
    std::size_t nstate = nm;
    for (std::size_t ti = 0; ti < t.transitions.size(); ++ti) {
        const std::size_t ind = t.transitions[ti];
        const typename std::map<std::size_t, qn::TransitionParam<T>>::const_iterator it =
            sn.transparam.find(ind);
        if (it == sn.transparam.end()) continue;
        for (std::size_t m = 0; m < it->second.nmodes; ++m) {
            PetriMode rec = terms_detail::build_mode(sn, it->second, ind, m, t, nm);
            if (rec.nph > 1) {
                rec.zblk.resize(rec.nph);
                for (std::size_t h = 0; h < rec.nph; ++h) rec.zblk[h] = nstate + h;
                nstate += rec.nph;
            }
            t.modes.push_back(rec);
        }
    }
    // ---- the fluid levels, one coordinate per fluid place after the phase block
    t.nf = t.fluid_places.size();
    t.fidx.assign(t.I, -1);
    for (std::size_t q = 0; q < t.nf; ++q) {
        const std::size_t ind = t.fluid_places[q];
        const qn::Station<T>& st = sn.stations[sn.nodes[ind - 1].station - 1];
        t.fidx[ind - 1] = static_cast<std::ptrdiff_t>(nstate + q);
        t.fluid_slot.push_back(nstate + q);
        t.fluid_cap.push_back(st.fluidcap);
        t.fluid_init.push_back(st.fluidinit);
        t.fluid_station.push_back(static_cast<std::ptrdiff_t>(sn.nodes[ind - 1].station) - 1);
        t.fluid_name.push_back(t.names_node[ind - 1]);
    }
    nstate += t.nf;
    t.nstate = nstate;
    // The phase coordinates are appended after every marking coordinate, so a
    // jump column built at nm width has to be GROWN once the total is known.
    for (std::size_t j = 0; j < t.modes.size(); ++j)
        if (t.modes[j].cvec.size() < nstate) t.modes[j].cvec.resize(nstate, 0.0);
    // The fluid arcs of each mode are resolved to level coordinates: an impulse is the fluid entry of the
    // firing column, a flow and a guard and a fluid rate dependence are read by index. `sn_fluid_refusal`
    // has already refused an arc naming a place that is not fluid, so every fidx read here is a coordinate.
    const auto level_of = [&](std::size_t node, const std::string& lbl) -> std::size_t {
        if (node == 0 || node > t.I || t.fidx[node - 1] < 0)
            throw InputError("fluid_petri_terms: mode " + lbl + " has a fluid arc on a node that is not a "
                             "fluid place");
        return static_cast<std::size_t>(t.fidx[node - 1]);
    };
    for (std::size_t j = 0; j < t.modes.size(); ++j) {
        PetriMode& md = t.modes[j];
        const qn::TransitionParam<T>& tp = sn.transparam.find(md.node)->second;
        const std::vector<typename qn::TransitionParam<T>::FluidArc>& imp = tp.fluid_impulses(md.mode);
        for (std::size_t a = 0; a < imp.size(); ++a) md.cvec[level_of(imp[a].node, md.label)] += imp[a].value;
        const std::vector<typename qn::TransitionParam<T>::FluidArc>& fr = tp.fluid_rates(md.mode);
        const std::vector<typename qn::TransitionParam<T>::FluidFlowDep>& fld = tp.fluid_flow_deps(md.mode);
        const std::size_t KK = t.K;
        for (std::size_t a = 0; a < fr.size(); ++a) {
            md.flow_slot.push_back(level_of(fr[a].node, md.label));
            md.flow_rate.push_back(fr[a].value);
            md.flow_dep_fun.push_back(std::function<double(const Matrix<double>&)>());
            md.flow_dep_read.push_back(std::vector<std::size_t>());
            for (std::size_t u = 0; u < fld.size(); ++u) {
                if (fld[u].node != fr[a].node) continue;
                // the setter's domain is the row-major flattening of the (I x K) marking
                const std::function<double(const std::vector<double>&)> phi = fld[u].fun;
                md.flow_dep_fun.back() = [phi, KK](const Matrix<double>& mm) -> double {
                    std::vector<double> flat(mm.rows() * KK, 0.0);
                    for (std::size_t i = 0; i < mm.rows(); ++i)
                        for (std::size_t k = 0; k < KK && k < mm.cols(); ++k) flat[i * KK + k] = mm(i, k);
                    return phi(flat);
                };
                md.flow_dep_read.back() = fld[u].read;
            }
        }
        const std::vector<typename qn::TransitionParam<T>::FluidGuard>& gd = tp.fluid_guards(md.mode);
        for (std::size_t a = 0; a < gd.size(); ++a) {
            md.guard_slot.push_back(level_of(gd[a].node, md.label));
            md.guard_sgn.push_back(gd[a].sgn);
            md.guard_thr.push_back(gd[a].thr);
        }
        const std::vector<typename qn::TransitionParam<T>::FluidDep>& fd = tp.fluid_deps(md.mode);
        for (std::size_t a = 0; a < fd.size(); ++a) {
            md.fdep_slot.push_back(level_of(fd[a].node, md.label));
            md.fdep_fun.push_back(fd[a].fun);
        }
    }
    for (std::size_t j = 0; j < t.modes.size(); ++j) {
        if (t.modes[j].timing == lang::TimingStrategy::TIMED) t.timed_idx.push_back(j);
        else t.imm_idx.push_back(j);
    }
 
    // ---- the event columns --------------------------------------------------
    std::vector<std::vector<double>> cols;
    std::vector<double> rate_base;
    std::vector<int> ev_kind, ev_mode, ev_phase, ev_to;
    std::vector<std::ptrdiff_t> ev_station, ev_class;
    std::vector<std::size_t> ev_flow;
    const auto emit = [&](const std::vector<double>& col, double base, int kind, int mode,
                          int phase, int to, std::ptrdiff_t station, std::ptrdiff_t cls) {
        cols.push_back(col);
        rate_base.push_back(base);
        ev_kind.push_back(kind);
        ev_mode.push_back(mode);
        ev_phase.push_back(phase);
        ev_to.push_back(to);
        ev_station.push_back(station);
        ev_class.push_back(cls);
    };
 
    for (std::size_t q = 0; q < t.timed_idx.size(); ++q) {
        const std::size_t j = t.timed_idx[q];
        const PetriMode& md = t.modes[j];
        if (md.nph == 1) {
            emit(md.cvec, md.d1[0], 1, static_cast<int>(j), 0, 0, -1, -1);
        } else {
            for (std::size_t h = 0; h < md.nph; ++h)
                for (std::size_t hp = 0; hp < md.nph; ++hp) {
                    const double w = md.D1(h, hp);
                    if (w <= 0) continue;
                    std::vector<double> col = md.cvec;
                    col[md.zblk[hp]] += 1.0;
                    col[md.zblk[h]] -= 1.0;
                    emit(col, w, 1, static_cast<int>(j), static_cast<int>(h),
                         static_cast<int>(hp), -1, -1);
                }
            for (std::size_t h = 0; h < md.nph; ++h)
                for (std::size_t hp = 0; hp < md.nph; ++hp) {
                    if (hp == h) continue;
                    const double w = md.D0(h, hp);
                    if (w <= 0) continue;
                    std::vector<double> col(nstate, 0.0);
                    col[md.zblk[hp]] = 1.0;
                    col[md.zblk[h]] = -1.0;
                    emit(col, w, 2, static_cast<int>(j), static_cast<int>(h),
                         static_cast<int>(hp), -1, -1);
                }
        }
    }
    // One latch column per multi-phase mode: mu_j servers per unit time enter at
    // the firing process's own entry distribution. The rate is FREE IN SIGN -- a
    // mode whose enabling degree drops stops servers rather than starting them --
    // and it is zero at any fixed point.
    for (std::size_t q = 0; q < t.timed_idx.size(); ++q) {
        const std::size_t j = t.timed_idx[q];
        const PetriMode& md = t.modes[j];
        if (md.nph <= 1) continue;
        std::vector<double> col(nstate, 0.0);
        for (std::size_t h = 0; h < md.nph; ++h) col[md.zblk[h]] = md.pie[h];
        emit(col, 1.0, 5, static_cast<int>(j), 0, 0, -1, -1);
    }
    for (std::size_t q = 0; q < t.imm_idx.size(); ++q) {
        const std::size_t j = t.imm_idx[q];
        emit(t.modes[j].cvec, 1.0, 4, static_cast<int>(j), 0, 0, -1, -1);
    }
    for (std::size_t a = 0; a < src_arr.size(); ++a) {
        const SrcArr& rec = src_arr[a];
        const std::size_t ist = sn.nodes[rec.snd - 1].station;
        if (ist >= 1 && ist <= sn.stations.size()) {
            const lang::ProcessType pt = sn.procid(ist, rec.r + 1);
            if (pt != lang::ProcessType::EXP)
                throw UnsupportedError(ErrorKind{"NonExponential"},
                    "fluid_petri_terms: source " + t.names_node[rec.snd - 1] +
                    " has a non-exponential arrival for class " + std::to_string(rec.r + 1) +
                    ". The fluid Petri route models an arrival as a constant-propensity event, "
                    "which a renewal stream with memory is not; use SolverCTMC, SolverJMT or "
                    "SolverSSA");
        }
        std::vector<double> col(nstate, 0.0);
        col[static_cast<std::size_t>(t.pidx[rec.qnd - 1][rec.l])] = 1.0;
        const double base =
            num_traits<T>::to_double(sn.rates(ist - 1, rec.r)) *
            num_traits<T>::to_double(
                sn.rtnodes((rec.snd - 1) * t.K + rec.r, (rec.qnd - 1) * t.K + rec.l));
        emit(col, base, 3, -1, 0, 0, static_cast<std::ptrdiff_t>(ist) - 1,
             static_cast<std::ptrdiff_t>(rec.r));
    }
    // One flow column per (timed mode, fluid place): a unit jump on the level, signed by the direction of the
    // flow, at |rate| per running server.
    for (std::size_t q = 0; q < t.timed_idx.size(); ++q) {
        const std::size_t j = t.timed_idx[q];
        const PetriMode& md = t.modes[j];
        for (std::size_t a = 0; a < md.flow_slot.size(); ++a) {
            std::vector<double> col(nstate, 0.0);
            col[md.flow_slot[a]] = md.flow_rate[a] > 0.0 ? 1.0 : (md.flow_rate[a] < 0.0 ? -1.0 : 0.0);
            emit(col, std::fabs(md.flow_rate[a]), 6, static_cast<int>(j), 0, 0, -1, -1);
            ev_flow.resize(rate_base.size(), 0);
            ev_flow.back() = a;
        }
    }
 
    t.nev = cols.size();
    t.D = Matrix<double>(nstate, std::max<std::size_t>(t.nev, 1), 0.0);
    for (std::size_t e = 0; e < t.nev; ++e)
        for (std::size_t s = 0; s < nstate; ++s)
            if (cols[e][s] != 0.0) t.D(s, e) = cols[e][s];
    t.rate_base = rate_base;
    t.ev_kind = ev_kind;
    t.ev_mode = ev_mode;
    t.ev_phase = ev_phase;
    t.ev_to = ev_to;
    t.ev_station = ev_station;
    t.ev_class = ev_class;
    ev_flow.resize(t.nev, 0);
    t.ev_flow = ev_flow;
 
    // ---- which Sigma entries the closure reads ------------------------------
    std::map<std::pair<std::size_t, std::size_t>, std::size_t> pairkey;
    for (std::size_t j = 0; j < t.modes.size(); ++j) {
        const PetriMode& md = t.modes[j];
        if (md.closable)
            for (std::size_t a = 0; a < md.arc_slot.size(); ++a)
                for (std::size_t b = a; b < md.arc_slot.size(); ++b)
                    terms_detail::add_pair(pairkey, t.cov_pairs, md.arc_slot[a], md.arc_slot[b]);
        if (md.timing == lang::TimingStrategy::TIMED) {
            for (std::size_t b = 0; b < md.inh_slot.size(); ++b)
                terms_detail::add_pair(pairkey, t.cov_pairs, md.inh_slot[b], md.inh_slot[b]);
            // a fluid guard is smoothed by the marginal of its level, as an inhibitor gate is by the marginal
            // of its place
            for (std::size_t b = 0; b < md.guard_slot.size(); ++b)
                terms_detail::add_pair(pairkey, t.cov_pairs, md.guard_slot[b], md.guard_slot[b]);
        }
    }
    t.npair = t.cov_pairs.size();
    const std::size_t pn = std::max<std::size_t>(nstate, 1);
    t.pair_index.assign(pn, std::vector<std::ptrdiff_t>(pn, -1));
    for (std::size_t i = 0; i < t.npair; ++i) {
        t.pair_index[t.cov_pairs[i].first][t.cov_pairs[i].second] =
            static_cast<std::ptrdiff_t>(i);
        t.pair_index[t.cov_pairs[i].second][t.cov_pairs[i].first] =
            static_cast<std::ptrdiff_t>(i);
    }
 
    for (std::size_t e = 0; e < t.nev; ++e) {
        if (t.ev_kind[e] == 4) t.imm_col.push_back(e);
        else if (t.ev_kind[e] == 5) t.latch_col.push_back(e);
        // A fluid flow is deterministic given the marking, so it has no noise of its own either: the level
        // inherits its fluctuation from the marking through the drift Jacobian, and from the impulses carried
        // on the firing columns.
        if (t.ev_kind[e] != 4 && t.ev_kind[e] != 5 && t.ev_kind[e] != 6) t.stoch_col.push_back(e);
        if (t.ev_kind[e] == 6) t.flow_col.push_back(e);
    }
    for (std::size_t q = 0; q < t.timed_idx.size(); ++q)
        if (t.modes[t.timed_idx[q]].nph > 1) t.latch_mode.push_back(t.timed_idx[q]);
    t.cov_idx.resize(nstate);
    for (std::size_t s = 0; s < nstate; ++s) t.cov_idx[s] = s;
 
    // ---- per-place consumption and production, for the metric reader --------
    // A Place's throughput is the rate at which TOKENS leave it, so each
    // consuming mode contributes its firing rate times the multiplicity of the
    // arc it takes them through. That is SolverCTMC's convention and the one
    // Little's law needs; SolverSSA's NRM sums the UNWEIGHTED propensity, so the
    // two disagree wherever an input arc has multiplicity above one.
    for (std::size_t e = 0; e < t.nev; ++e) {
        if (t.ev_kind[e] == 3) {
            if (t.ev_station[e] < 0) continue;
            const std::size_t kk =
                static_cast<std::size_t>(t.ev_station[e]) * t.K +
                static_cast<std::size_t>(t.ev_class[e]);
            t.producers[kk].push_back(e);
            continue;
        }
        if (t.ev_kind[e] != 1 && t.ev_kind[e] != 4) continue;
        const PetriMode& md = t.modes[static_cast<std::size_t>(t.ev_mode[e])];
        for (std::size_t a = 0; a < md.arc_slot.size(); ++a) {
            const std::size_t s = md.arc_slot[a];
            if (t.coord_station[s] < 0) continue;
            const std::size_t kk =
                static_cast<std::size_t>(t.coord_station[s]) * t.K + t.coord_class[s];
            t.consumers[kk].push_back(e);
            t.consumer_w[kk].push_back(md.arc_w[a]);
        }
    }
 
    // ---- the initial state --------------------------------------------------
    t.x0.assign(nstate, 0.0);
    for (std::size_t s = 0; s < nm; ++s) t.x0[s] = t.m0full[t.coord_node[s]][t.coord_class[s]];
    for (std::size_t q = 0; q < t.nf; ++q) t.x0[t.fluid_slot[q]] = t.fluid_init[q];
    for (std::size_t j = 0; j < t.modes.size(); ++j) {
        const PetriMode& md = t.modes[j];
        if (md.nph <= 1) continue;
        double e = std::numeric_limits<double>::infinity();
        for (std::size_t a = 0; a < md.arc_slot.size(); ++a)
            e = std::min(e, t.x0[md.arc_slot[a]] / md.arc_w[a]);
        if (md.arc_slot.empty()) e = 1.0;
        for (std::size_t h = 0; h < md.nph; ++h) t.x0[md.zblk[h]] = std::min(e, md.c) * md.pie[h];
    }
    return t;
}
 
}  // namespace petri
}  // namespace fluid
}  // namespace line
 
#endif  // LINE_SOLVERS_FLUID_PETRI_TERMS_H