fluid_petri.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_H
#define LINE_SOLVERS_FLUID_PETRI_H
 
/**
 * @file
 * @ingroup line_solvers
 * Fluid analysis of a stochastic Petri net: one simultaneous algebraic solve per
 * active set.
 *
 * Port of `matlab/src/solvers/FLD/solver_fluid_petri.m`, together with
 * `fluid_petri_conservation.m`, `fluid_petri_constraints.m`,
 * `fluid_petri_immediate.m` and `fluid_petri_applicable.m`. Cross-checked
 * against `jar/src/main/java/jline/solvers/fluid/petri/PetriSolver.java` and the
 * python `solver_fld/methods/petri.py`.
 *
 * A net's fluid limit is NOT the queueing drift with places in place of
 * stations. Its immediate transitions have no rate at all: their limit is a
 * FLOW, an algebraic unknown pinned by the constraint that the input place they
 * bind holds no mass. The marking, the diffusion covariance, those flows, the
 * multi-server phase latches and the capacity gates therefore solve
 * SIMULTANEOUSLY, as one algebraic system per active set, rather than by
 * integrating an ODE to its fixed point.
 *
 * THE ACTIVE SET IS WHAT ITERATES. A pass solves the system for a fixed choice
 * of which immediate modes fire, which coordinate each pins, and which
 * capacities bind; the answer says whether that choice was right (a negative
 * flow, a negative marking, a violated or a released capacity), and the next
 * pass makes ONE move. The moves are ordered so that a failure of each
 * invalidates the next.
 */
 
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstddef>
#include <cstdio>
#include <limits>
#include <map>
#include <set>
#include <string>
#include <vector>
 
#include "line/api/sn/sn_fluid_refusal.h"
#include "line/solvers/fluid/fluid_moments.h"
#include "line/solvers/fluid/fluid_petri_system.h"
#include "line/solvers/fluid/fluid_petri_terms.h"
#include "line/util/error.h"
#include "line/util/linalg.h"
#include "line/util/matrix.h"
#include "line/util/ode.h"
#include "line/util/svd.h"
#include "rodas.hpp"
 
namespace line {
namespace fluid {
namespace petri {
 
// ============================ conservation ===============================
 
/** The conserved quantities of a net, as equations. */
struct PetriConservation {
    Matrix<double> C;
    std::vector<double> N;
    double leak = 0.0;
    std::vector<std::string> label;
};
 
namespace cons_detail {
 
/**
 * A rational basis of the null space, as MATLAB's `null(A,'r')` returns.
 *
 * D is INTEGRAL -- arc multiplicities and unit phase moves -- so the reduced row
 * echelon form gives exact rational rows, which keeps each conservation row
 * readable as a statement about named places instead of an arbitrary orthogonal
 * mixture. An SVD basis would span the same space and say nothing.
 */
inline Matrix<double> null_rational(const Matrix<double>& A) {
    const std::size_t rows = A.rows(), cols = A.cols();
    if (rows == 0 || cols == 0) return Matrix<double>(cols, 0, 0.0);
    Matrix<double> R = A;
    std::vector<std::size_t> piv;
    std::size_t r = 0;
    const double tol = 1e-12;
    for (std::size_t c = 0; c < cols && r < rows; ++c) {
        std::size_t k = r;
        double best = std::fabs(R(r, c));
        for (std::size_t i = r + 1; i < rows; ++i)
            if (std::fabs(R(i, c)) > best) {
                best = std::fabs(R(i, c));
                k = i;
            }
        if (best <= tol) {
            for (std::size_t i = r; i < rows; ++i) R(i, c) = 0.0;
            continue;
        }
        if (k != r)
            for (std::size_t j = 0; j < cols; ++j) std::swap(R(r, j), R(k, j));
        const double p = R(r, c);
        for (std::size_t j = 0; j < cols; ++j) R(r, j) /= p;
        for (std::size_t i = 0; i < rows; ++i) {
            if (i == r) continue;
            const double f = R(i, c);
            if (f == 0.0) continue;
            for (std::size_t j = 0; j < cols; ++j) R(i, j) -= f * R(r, j);
        }
        piv.push_back(c);
        ++r;
    }
    std::vector<std::size_t> free;
    for (std::size_t c = 0; c < cols; ++c)
        if (std::find(piv.begin(), piv.end(), c) == piv.end()) free.push_back(c);
    Matrix<double> Z(cols, free.size(), 0.0);
    for (std::size_t i = 0; i < free.size(); ++i) {
        Z(free[i], i) = 1.0;
        for (std::size_t rr = 0; rr < piv.size(); ++rr) Z(piv[rr], i) = -R(rr, free[i]);
    }
    return Z;
}
 
/** %g formatting, so an integral weight prints without a decimal tail. */
inline std::string trim(double v) {
    if (v == std::rint(v) && std::fabs(v) < 1e15)
        return std::to_string(static_cast<long long>(v));
    char buf[32];
    std::snprintf(buf, sizeof(buf), "%g", v);
    return std::string(buf);
}
 
/** One conservation row, written out over the coordinates it touches. */
inline std::string row_label(const PetriTerms& t, const Matrix<double>& C, std::size_t r) {
    std::string out;
    for (std::size_t s = 0; s < C.cols(); ++s) {
        const double w = C(r, s);
        if (w == 0.0) continue;
        std::string nm;
        if (s < t.nm) {
            nm = t.names_node[t.coord_node[s]] + "(class " +
                 std::to_string(t.coord_class[s] + 1) + ")";
        } else {
            nm = "phase";
            for (std::size_t q = 0; q < t.nf; ++q)
                if (t.fluid_slot[q] == s) nm = "level of " + t.fluid_name[q];
            for (std::size_t j = 0; j < t.modes.size(); ++j)
                for (std::size_t q = 0; q < t.modes[j].zblk.size(); ++q)
                    if (t.modes[j].zblk[q] == s)
                        nm = t.modes[j].label + " phase " + std::to_string(q + 1);
        }
        if (!out.empty()) out += " + ";
        out += (w == 1.0) ? nm : (trim(w) + "*" + nm);
    }
    return out;
}
 
}  // namespace cons_detail
 
/**
 * The conserved quantities, as equations: `u'D = 0` implies `u'x` is constant.
 *
 * On the marking coordinates those u are the net's P-invariants; on a mode's
 * phase block the all-ones vector is one of them, which is the statement that
 * the phase coordinates are a distribution. Both come out of the SAME null
 * space, so the phase normalisation needs no separate row.
 *
 * AN OPEN NET LOSES THE ROWS ITS ARRIVALS BREAK, automatically: the arrival
 * columns are part of D, so a u an arrival moves is not in the null space.
 */
inline PetriConservation petri_conservation(const PetriTerms& t) {
    PetriConservation cons;
    if (t.D.rows() == 0 || t.D.cols() == 0 || t.nev == 0) {
        cons.C = Matrix<double>(0, t.nstate, 0.0);
        return cons;
    }
    Matrix<double> Dt(t.D.cols(), t.D.rows(), 0.0);
    for (std::size_t i = 0; i < t.D.rows(); ++i)
        for (std::size_t j = 0; j < t.D.cols(); ++j) Dt(j, i) = t.D(i, j);
    const Matrix<double> Z = cons_detail::null_rational(Dt);
    Matrix<double> C(Z.cols(), Z.rows(), 0.0);
    for (std::size_t i = 0; i < Z.rows(); ++i)
        for (std::size_t j = 0; j < Z.cols(); ++j) C(j, i) = Z(i, j);
    // A rational basis has exact zeros; anything below this is a rounding
    // artefact of the elimination, not a coefficient.
    for (std::size_t i = 0; i < C.rows(); ++i)
        for (std::size_t j = 0; j < C.cols(); ++j)
            if (std::fabs(C(i, j)) < 1e-12) C(i, j) = 0.0;
    // Scale each row by its smallest nonzero magnitude, so a row reads as a
    // statement about named places with small integer weights.
    std::vector<std::size_t> keep;
    for (std::size_t i = 0; i < C.rows(); ++i) {
        double smallest = std::numeric_limits<double>::infinity();
        for (std::size_t j = 0; j < C.cols(); ++j)
            if (C(i, j) != 0.0) smallest = std::min(smallest, std::fabs(C(i, j)));
        if (!std::isfinite(smallest)) continue;
        for (std::size_t j = 0; j < C.cols(); ++j) C(i, j) /= smallest;
        keep.push_back(i);
    }
    cons.C = Matrix<double>(keep.size(), t.nstate, 0.0);
    for (std::size_t r = 0; r < keep.size(); ++r)
        for (std::size_t j = 0; j < t.nstate; ++j) cons.C(r, j) = C(keep[r], j);
    cons.N.assign(cons.C.rows(), 0.0);
    for (std::size_t r = 0; r < cons.C.rows(); ++r) {
        double v = 0.0;
        for (std::size_t j = 0; j < t.nstate; ++j) v += cons.C(r, j) * t.x0[j];
        cons.N[r] = v;
    }
    for (std::size_t r = 0; r < cons.C.rows(); ++r)
        for (std::size_t c = 0; c < t.nev; ++c) {
            double v = 0.0;
            for (std::size_t j = 0; j < t.nstate; ++j) v += cons.C(r, j) * t.D(j, c);
            cons.leak = std::max(cons.leak, std::fabs(v));
        }
    for (std::size_t r = 0; r < cons.C.rows(); ++r)
        cons.label.push_back(cons_detail::row_label(t, cons.C, r));
    return cons;
}
 
// ============================ capacities =================================
 
/**
 * Every finite place capacity as a linear row `A x <= b`, and the two bounds of every fluid place.
 *
 * `cover[c]` names the coordinates row c caps, i.e. the deposits (or, for a LOWER row, the withdrawals) its
 * multiplier throttles. `lower[c]` marks the empty-level row of a fluid place, whose multiplier gates the
 * withdrawals, and `fluid[c]` any bound of a fluid place.
 */
struct PetriConstraints {
    Matrix<double> A;
    std::vector<double> b;
    std::vector<std::string> label;
    std::vector<std::vector<bool>> cover;
    std::vector<bool> lower;
    std::vector<bool> fluid;
};
 
/**
 * THE GATE IS A LOSS ON THE DEPOSIT: LINE loses the tokens a firing would push
 * past a place's capacity, so the fluid analogue scales the DEPOSIT leg of every
 * event adding mass to the capped place and leaves the removal leg alone. The
 * rows here name which coordinates are capped; the scaling lives in the tangent
 * clamp.
 *
 * A FLUID PLACE HAS TWO BOUNDS, and the lower one is new. Its capacity B is a row of the same kind, throttling
 * the inflow legs (flows into it and positive impulses) while it binds. Its EMPTY level is a row the token places
 * never needed: a token drain slows down with the marking through min(), so a place empties asymptotically, but
 * a fluid drain runs at its full rate down to zero. The row -f <= 0 throttles the OUTFLOW legs (flows out of it
 * and negative impulses) while it binds: an empty place supplies no more than its inflow, and a firing is never
 * blocked, only its withdrawal truncated.
 */
template <class T>
inline PetriConstraints petri_constraints(const qn::NetworkStruct<T>& sn, const PetriTerms& t) {
    std::vector<std::vector<double>> rows;
    std::vector<double> bs;
    std::vector<std::string> labels;
    for (std::size_t pi = 0; pi < t.places.size(); ++pi) {
        const std::size_t ind = t.places[pi];
        const std::size_t ist = sn.nodes[ind - 1].station;
        std::vector<std::size_t> slots;
        for (std::size_t k = 0; k < t.K; ++k)
            if (t.pidx[ind - 1][k] >= 0)
                slots.push_back(static_cast<std::size_t>(t.pidx[ind - 1][k]));
        if (slots.empty()) continue;
        if (ist >= 1 && ist <= sn.stations.size()) {
            const double cap = sn.stations[ist - 1].cap;
            if (std::isfinite(cap)) {
                std::vector<double> row(t.nstate, 0.0);
                for (std::size_t s : slots) row[s] = 1.0;
                rows.push_back(row);
                bs.push_back(cap);
                labels.push_back("capacity " + cons_detail::trim(cap) + " of place " +
                                 t.names_node[ind - 1]);
            }
            for (std::size_t k = 0; k < t.K; ++k) {
                if (t.pidx[ind - 1][k] < 0) continue;
                if (k >= sn.stations[ist - 1].classcap.size()) continue;
                const double ccap = sn.stations[ist - 1].classcap[k];
                if (!std::isfinite(ccap)) continue;
                std::vector<double> row(t.nstate, 0.0);
                row[static_cast<std::size_t>(t.pidx[ind - 1][k])] = 1.0;
                rows.push_back(row);
                bs.push_back(ccap);
                labels.push_back("class-" + std::to_string(k + 1) + " capacity " +
                                 cons_detail::trim(ccap) + " of place " + t.names_node[ind - 1]);
            }
        }
    }
    // TWO ROWS THAT SAY THE SAME THING ARE A SINGULAR NEWTON SYSTEM, not a
    // redundancy the least squares absorbs, so an exact duplicate is pruned and
    // the tighter bound survives.
    std::vector<bool> keep(rows.size(), true);
    for (std::size_t c = 0; c < rows.size(); ++c) {
        if (!keep[c]) continue;
        for (std::size_t d = c + 1; d < rows.size(); ++d) {
            if (!keep[d] || rows[c] != rows[d]) continue;
            if (bs[d] < bs[c]) {
                keep[c] = false;
                break;
            }
            keep[d] = false;
        }
    }
    std::vector<std::size_t> idx;
    for (std::size_t i = 0; i < keep.size(); ++i)
        if (keep[i]) idx.push_back(i);
    std::vector<bool> lowers(idx.size(), false), fluids(idx.size(), false);
    for (std::size_t q = 0; q < t.nf; ++q) {
        const std::size_t slot = t.fluid_slot[q];
        if (std::isfinite(t.fluid_cap[q])) {
            std::vector<double> row(t.nstate, 0.0);
            row[slot] = 1.0;
            rows.push_back(row);
            bs.push_back(t.fluid_cap[q]);
            char buf[64];
            std::snprintf(buf, sizeof(buf), "%g", t.fluid_cap[q]);
            labels.push_back("fluid capacity " + std::string(buf) + " of place " + t.fluid_name[q]);
            idx.push_back(rows.size() - 1);
            lowers.push_back(false);
            fluids.push_back(true);
        }
        std::vector<double> row(t.nstate, 0.0);
        row[slot] = -1.0;
        rows.push_back(row);
        bs.push_back(0.0);
        labels.push_back("empty level of fluid place " + t.fluid_name[q]);
        idx.push_back(rows.size() - 1);
        lowers.push_back(true);
        fluids.push_back(true);
    }
    PetriConstraints con;
    con.A = Matrix<double>(idx.size(), t.nstate, 0.0);
    con.b.assign(idx.size(), 0.0);
    con.cover.assign(idx.size(), std::vector<bool>(t.nstate, false));
    for (std::size_t r = 0; r < idx.size(); ++r) {
        for (std::size_t s = 0; s < t.nstate; ++s) {
            con.A(r, s) = rows[idx[r]][s];
            con.cover[r][s] = rows[idx[r]][s] != 0.0;
        }
        con.b[r] = bs[idx[r]];
        con.label.push_back(labels[idx[r]]);
    }
    con.lower = lowers;
    con.fluid = fluids;
    return con;
}
 
// ============================ the immediates =============================
 
/** The active set of the immediate modes and the equations pinning their flows. */
struct PetriImmediate {
    enum Kind { PIN = 0, RATIO = 1, ZERO = 2 };
    struct Row {
        Kind kind = PIN;
        std::size_t a = 0, b = 0;
        double wa = 0.0, wb = 0.0;
    };
    std::size_t n = 0;
    std::vector<bool> active;
    std::vector<std::ptrdiff_t> bind;
    std::vector<std::size_t> pins;
    std::vector<Row> rows;
    bool initialized = false;
};
 
namespace imm_detail {
 
/**
 * True when an inhibitor arc of this mode has reached its threshold.
 *
 * A HARD TEST ON THE MEAN, AND A KNOWN WRONG ANSWER WHEN THE MEAN SITS ON THE
 * THRESHOLD -- a timed mode closes the same indicator smoothly, as
 * `Phi((thr-m)/sd)` in `petri_theta`; this path does not. See
 * `_kb/06-solver-catalog.md` for the measurement and for what a real fix costs.
 */
inline bool inhibited(const PetriMode& md, const std::vector<double>& x) {
    for (std::size_t b = 0; b < md.inh_slot.size(); ++b)
        if (x[md.inh_slot[b]] >= md.inh_thr[b]) return true;
    return false;
}
 
}  // namespace imm_detail
 
/**
 * An immediate transition has no rate: its fluid limit is a FLOW, an algebraic
 * unknown pinned by the constraint that its binding input place holds no mass,
 *
 *     phi_j >= 0,   x_b = 0 for the coordinate b that binds mode j
 *
 * with the GSPN conflict rule supplying the extra equation when two modes drain
 * one place: `phi_j*weight_l = phi_l*weight_j` among the enabled modes of highest
 * firing priority, and `phi = 0` below it.
 *
 * THE COUNT IS SQUARE BY CONSTRUCTION: V pins plus (F-V) ratio rows is F
 * equations for F flows.
 */
inline void petri_immediate(const PetriTerms& t, const std::vector<double>& x,
                            PetriImmediate& imm) {
    const std::size_t n = t.imm_idx.size();
    if (!imm.initialized) {
        imm.n = n;
        imm.active.assign(n, true);
        imm.bind.assign(n, -1);
        imm.pins.clear();
        // An inhibited mode never fires, so it neither carries a flow nor empties
        // a place. AN EMPTY INPUT PLACE IS NOT A REASON TO DEACTIVATE -- that is
        // the normal state of an enabled immediate mode, and the whole content of
        // its pin.
        for (std::size_t k = 0; k < n; ++k) {
            const PetriMode& md = t.modes[t.imm_idx[k]];
            if (md.arc_slot.empty())
                throw InputError("fluid_petri_immediate: immediate mode " + md.label +
                                 " has no enabling arc, so nothing bounds its firing flow and "
                                 "the net has no fluid limit. Give it an input place, or make "
                                 "it timed");
            if (imm_detail::inhibited(md, x)) imm.active[k] = false;
        }
        imm.initialized = true;
    }
 
    // ---- the assignment: each active mode binds the input arc it is shortest of
    for (std::size_t k = 0; k < n; ++k) {
        if (!imm.active[k]) {
            imm.bind[k] = -1;
            continue;
        }
        const PetriMode& md = t.modes[t.imm_idx[k]];
        if (imm.bind[k] >= 0 &&
            std::find(md.arc_slot.begin(), md.arc_slot.end(),
                      static_cast<std::size_t>(imm.bind[k])) != md.arc_slot.end())
            continue;  // a binding the caller set explicitly is kept
        std::size_t best = md.arc_slot[0];
        double best_lev = x[md.arc_slot[0]] / md.arc_w[0];
        for (std::size_t a = 1; a < md.arc_slot.size(); ++a) {
            const double lev = x[md.arc_slot[a]] / md.arc_w[a];
            if (lev < best_lev) {
                best_lev = lev;
                best = md.arc_slot[a];
            }
        }
        imm.bind[k] = static_cast<std::ptrdiff_t>(best);
    }
 
    // ---- the equations
    std::set<std::size_t> pinset;
    for (std::size_t k = 0; k < n; ++k)
        if (imm.active[k] && imm.bind[k] >= 0)
            pinset.insert(static_cast<std::size_t>(imm.bind[k]));
    imm.pins.assign(pinset.begin(), pinset.end());
    imm.rows.clear();
    for (std::size_t p : imm.pins) {
        PetriImmediate::Row row;
        row.kind = PetriImmediate::PIN;
        row.a = p;
        imm.rows.push_back(row);
        std::vector<std::size_t> grp;
        for (std::size_t k = 0; k < n; ++k)
            if (imm.active[k] && imm.bind[k] == static_cast<std::ptrdiff_t>(p)) grp.push_back(k);
        if (grp.size() <= 1) continue;
        int topprio = std::numeric_limits<int>::min();
        for (std::size_t k : grp) topprio = std::max(topprio, t.modes[t.imm_idx[k]].prio);
        std::vector<std::size_t> top, low;
        for (std::size_t k : grp)
            (t.modes[t.imm_idx[k]].prio == topprio ? top : low).push_back(k);
        const double w0 = t.modes[t.imm_idx[top[0]]].weight;
        for (std::size_t i = 1; i < top.size(); ++i) {
            PetriImmediate::Row rr;
            rr.kind = PetriImmediate::RATIO;
            rr.a = top[0];
            rr.b = top[i];
            rr.wa = w0;
            rr.wb = t.modes[t.imm_idx[top[i]]].weight;
            imm.rows.push_back(rr);
        }
        for (std::size_t k : low) {
            PetriImmediate::Row rr;
            rr.kind = PetriImmediate::ZERO;
            rr.a = k;
            imm.rows.push_back(rr);
        }
    }
    for (std::size_t k = 0; k < n; ++k)
        if (!imm.active[k]) {
            PetriImmediate::Row rr;
            rr.kind = PetriImmediate::ZERO;
            rr.a = k;
            imm.rows.push_back(rr);
        }
}
 
// ============================ the refusals ===============================
 
/** Whether the fluid Petri route can answer this model, and why not. */
struct PetriVerdict {
    bool ok = true;
    std::string reason;
};
 
/**
 * A QUEUEING STATION IS THE ONE STRUCTURAL EXCLUSION: a net whose tokens also
 * visit a Queue or a Delay is two formalisms at once, and LINE has no reference
 * semantics for the hand-off.
 */
template <class T>
inline PetriVerdict petri_applicable(const qn::NetworkStruct<T>& sn) {
    PetriVerdict v;
    const std::size_t I = sn.nodes.size();
    bool has_transition = false;
    for (std::size_t i = 0; i < I; ++i)
        if (sn.nodes[i].nodetype == lang::NodeType::Transition) has_transition = true;
    if (!has_transition) {
        v.ok = false;
        v.reason = "the model has no Transition node, so it is not a Petri net";
        return v;
    }
    for (std::size_t i = 0; i < I; ++i) {
        const lang::NodeType nt = sn.nodes[i].nodetype;
        if (nt != lang::NodeType::Place && nt != lang::NodeType::Transition &&
            nt != lang::NodeType::Source && nt != lang::NodeType::Sink) {
            v.ok = false;
            v.reason = "node " + sn.nodes[i].name + " is a " +
                       std::string(lang::node_type_to_text(nt)) +
                       ". The fluid Petri route solves the marking of a Petri net, and a model "
                       "that also holds queueing stations is two formalisms at once with no "
                       "reference semantics for the hand-off; use SolverCTMC, SolverJMT, "
                       "SolverSSA or SolverLDES";
            return v;
        }
    }
    // A queueing place declares a service process, which is what turns it into a
    // station with an embedded queue and a depository.
    for (std::size_t i = 0; i < I; ++i) {
        if (sn.nodes[i].nodetype != lang::NodeType::Place) continue;
        const std::size_t ist = sn.nodes[i].station;
        if (ist < 1 || ist > sn.nstations) continue;
        for (std::size_t k = 0; k < sn.nclasses; ++k) {
            const double rate = num_traits<T>::to_double(sn.rates(ist - 1, k));
            if (!std::isnan(rate) && rate > 0) {
                v.ok = false;
                v.reason = "place " + sn.nodes[i].name +
                           " is a QUEUEING place (it declares a service process), whose embedded "
                           "queue this drift does not carry; use SolverLDES";
                return v;
            }
        }
    }
    // -- fluid places (FSPN): a fluid place holds a level, not tokens, so every way of putting a token on it is
    // refused, and so is a fluid guard an immediate mode carries: a guard on a mode that fires in zero time is
    // a sliding mode on the guard surface, which a mean-field drift has no value for. The PDMP simulator of
    // SolverSSA does answer it.
    const api::SnFluidVerdict fv = api::sn_fluid_refusal(sn);
    if (!fv.ok) {
        v.ok = false;
        v.reason = fv.reason;
        return v;
    }
    for (const auto& kv : sn.transparam) {
        const qn::TransitionParam<T>& tp = kv.second;
        for (std::size_t m = 0; m < tp.nmodes; ++m) {
            if (!(m < tp.timing.size() && tp.timing[m] == lang::TimingStrategy::IMMEDIATE)) continue;
            if (tp.fluid_guards(m).empty()) continue;
            const std::string mname = (m < tp.modenames.size() && !tp.modenames[m].empty())
                                          ? tp.modenames[m]
                                          : ("Mode" + std::to_string(m + 1));
            v.ok = false;
            v.reason = "immediate mode " + mname + " of transition " + sn.nodes[kv.first - 1].name +
                       " carries a fluid guard, which fires it on the guard surface itself: a sliding mode the "
                       "mean-field drift has no value for; use SolverSSA";
            return v;
        }
    }
    return v;
}
 
// ============================ the driver =================================
 
/** What the Petri route computes that the station table has no column for. */
struct PetriReport {
    Matrix<double> marking, marking_var;
    std::vector<std::string> mode_label;
    std::vector<double> mode_flow, immediate_flow;
    std::vector<std::string> invariant_label;
    std::vector<double> invariant_value, invariant_error;
    std::vector<std::string> capacity_label;
    std::vector<std::size_t> capacity_active;
    std::vector<double> capacity_fraction;
    std::vector<std::size_t> pinned;
    Matrix<double> Sigma;
    /**
     * MATLAB `moments.petri.fluid`: per fluid place its name, reported level and level variance, its
     * unthrottled mean drift, the variance rate of its net inflow (filled only by the diffusion refinement),
     * the bound it is held at (+1 capacity, -1 empty, 0 none), and the `fluid_boundary` mode in force.
     */
    struct Fluid {
        std::vector<std::string> names;
        std::vector<double> level, level_var, drift, sigma2;
        std::vector<int> bound;
        std::string boundary = "reflect";
    } fluid;
};
 
/** Everything `solver_fluid_petri` returns. */
struct PetriSolution {
    Matrix<double> QN, UN, RN, TN;
    std::vector<double> t;
    std::vector<std::vector<double>> xvec_t;
    /** The trajectory tables, indexed [station][class][time point]; filled along `t`. */
    std::vector<std::vector<std::vector<double>>> QNt, UNt, TNt;
    /**
     * The transient marking covariance at each `tvar` instant, over the `cov_idx` coordinates (row-major nc*nc),
     * filled only when a finite horizon was asked for and the covariance fits `dae_maxcov`.
     */
    std::vector<Matrix<double>> Sigmat;
    std::vector<double> tvar;
    /** The active-set switches of the transient: (time, event index), MATLAB's `switches`. */
    std::vector<std::pair<double, int>> switches;
    std::vector<double> x;
    std::size_t iters = 0;
    double resnorm = std::numeric_limits<double>::infinity();
    bool converged = false;
    double runtime = 0.0;
    Matrix<double> Sigma, QVar, QStd;
    PetriReport petri;
    std::vector<std::string> warnings;
};
 
/** Tuning of the outer solve. */
struct PetriOptions {
    std::size_t dae_maxstate = 100;
    double tol = 1e-8;
    std::size_t newton_max = 50;
    /**
     * `options.config.fluid_boundary`: how a fluid level held at a bound is reported, `reflect` (the bound
     * itself, the fixed point) or `diffusion` (the mean of a reflected Brownian motion, see `fluid_levels`).
     */
    std::string fluid_boundary = "reflect";
    /**
     * `options.timespan`: a finite end integrates the transient DAE over [start, end] (MATLAB `local_transient`);
     * an infinite one reports the seed path, as before.
     */
    double timespan_start = 0.0;
    double timespan_end = std::numeric_limits<double>::infinity();
    /** `options.tranpoints`: instants added to the integrator's own steps, read off its dense output. */
    std::vector<double> tranpoints;
    /** `options.config.dae_maxcov`: the largest covariance block integrated along the trajectory. */
    std::size_t dae_maxcov = 25;
};
 
namespace solver_detail {
 
/** Everything the residual needs that does not change within one active set. */
struct Ctx {
    const PetriTerms* terms = nullptr;
    const PetriImmediate* imm = nullptr;
    std::vector<std::size_t> active;
    const PetriConstraints* con = nullptr;
    Matrix<double> C;
    std::vector<double> N;
    Matrix<double> Dp, Dn;
};
 
/**
 * THE CONSERVATION ROWS A BINDING CAP BREAKS ARE DROPPED: a capped place loses
 * the tokens that do not fit, so a conserved quantity supported on it is not
 * conserved while the cap binds, and keeping its row would state an equation the
 * drift contradicts -- a singular Newton system rather than an inaccuracy.
 */
inline Ctx context(const PetriTerms& terms, const PetriConservation& cons,
                   const PetriConstraints& con, const PetriImmediate& imm,
                   const std::vector<std::size_t>& active) {
    Ctx ctx;
    ctx.terms = &terms;
    ctx.imm = &imm;
    ctx.active = active;
    ctx.con = &con;
    Matrix<double> C = cons.C;
    std::vector<double> N = cons.N;
    if (!active.empty() && C.rows() > 0) {
        std::vector<bool> hit(terms.nstate, false);
        for (std::size_t c : active)
            for (std::size_t s = 0; s < terms.nstate; ++s)
                if (con.cover[c][s]) hit[s] = true;
        std::vector<std::size_t> keep;
        for (std::size_t r = 0; r < C.rows(); ++r) {
            bool drop = false;
            for (std::size_t s = 0; s < terms.nstate && !drop; ++s)
                if (hit[s] && C(r, s) != 0.0) drop = true;
            if (!drop) keep.push_back(r);
        }
        Matrix<double> Ck(keep.size(), terms.nstate, 0.0);
        std::vector<double> Nk(keep.size(), 0.0);
        for (std::size_t i = 0; i < keep.size(); ++i) {
            for (std::size_t s = 0; s < terms.nstate; ++s) Ck(i, s) = C(keep[i], s);
            Nk[i] = N[keep[i]];
        }
        C = Ck;
        N = Nk;
    }
    ctx.C = C;
    ctx.N = N;
    ctx.Dp = Matrix<double>(terms.D.rows(), terms.D.cols(), 0.0);
    ctx.Dn = Matrix<double>(terms.D.rows(), terms.D.cols(), 0.0);
    for (std::size_t i = 0; i < terms.D.rows(); ++i)
        for (std::size_t j = 0; j < terms.D.cols(); ++j) {
            ctx.Dp(i, j) = std::max(terms.D(i, j), 0.0);
            ctx.Dn(i, j) = std::min(terms.D(i, j), 0.0);
        }
    return ctx;
}
 
/** The unknown vector, split into its blocks. */
struct Unpacked {
    std::vector<double> x, s2, phi, mu, zeta;
};
 
inline Unpacked unpack(const std::vector<double>& u, const PetriTerms& t,
                       const PetriImmediate& imm, std::size_t na) {
    const std::size_t n = t.nstate, npair = t.npair, ni = imm.n, nl = t.latch_mode.size();
    Unpacked up;
    up.x.assign(u.begin(), u.begin() + n);
    up.s2.assign(u.begin() + n, u.begin() + n + npair);
    up.phi.assign(u.begin() + n + npair, u.begin() + n + npair + ni);
    for (std::size_t i = 0; i < up.phi.size(); ++i) up.phi[i] = std::max(0.0, up.phi[i]);
    up.mu.assign(u.begin() + n + npair + ni, u.begin() + n + npair + ni + nl);
    up.zeta.assign(u.begin() + n + npair + ni + nl, u.begin() + n + npair + ni + nl + na);
    for (std::size_t i = 0; i < up.zeta.size(); ++i) up.zeta[i] = std::max(0.0, up.zeta[i]);
    return up;
}
 
inline std::ptrdiff_t index_of(const std::vector<std::size_t>& a, std::size_t v) {
    for (std::size_t i = 0; i < a.size(); ++i)
        if (a[i] == v) return static_cast<std::ptrdiff_t>(i);
    return -1;
}
 
/**
 * The reduction of the fluctuation onto the manifold the fast and clamped
 * directions leave free.
 *
 * AN IMMEDIATE PIN REDUCES OBLIQUELY, ALONG THE FAST REACTION ITSELF, and this
 * is the one place where the orthogonal projector the queueing twin uses is
 * WRONG rather than merely different: a slow event depositing into a pinned
 * place is answered instantly by the immediate transition, so its effective jump
 * is its own plus the immediate flow it triggers -- the token is forwarded, not
 * lost. An orthogonal projection deletes the deposit and destroys mass in the
 * diffusion.
 *
 *     P = I - Cf * G * E_B,   G = G0 * (E_B Cf G0)^-1
 *
 * A CAPACITY CAP REDUCES ORTHOGONALLY -- mass that does not fit is genuinely
 * lost, so there is nothing to forward it to. A SERVER LATCH REDUCES
 * ORTHOGONALLY TOO, on the LINEARISED row `[-dtheta_j/dm, 1 over the phase
 * block]`, which is why this projector depends on the iterate.
 *
 * THE PSEUDO-INVERSE IS THE RANK-REVEALING ONE, not a regularised normal-equation
 * solve: the oblique reduction must ANNIHILATE the pinned coordinate exactly, and
 * a Tikhonov term of 1e-12 leaves it at 1e-12, after which the reduced generator
 * keeps a marginal eigenvalue and `fluid_lyapunov` refuses the fixed point as
 * non-hyperbolic -- a failure to reduce reported as a property of the model.
 */
inline bool clamp_tangent(const PetriTerms& terms, const PetriImmediate& imm,
                          const PetriConstraints& con, const std::vector<std::size_t>& active,
                          const PetriTheta* th, Matrix<double>& T) {
    const std::vector<std::size_t>& idx = terms.cov_idx;
    const std::size_t nc = idx.size();
    T = Matrix<double>(nc, nc, 0.0);
    for (std::size_t i = 0; i < nc; ++i) T(i, i) = 1.0;
 
    std::vector<std::size_t> actk;
    for (std::size_t k = 0; k < imm.n; ++k)
        if (imm.active[k] && imm.bind[k] >= 0) actk.push_back(k);
    const std::vector<std::size_t>& B = imm.pins;
    if (!actk.empty() && !B.empty()) {
        Matrix<double> Cf(nc, actk.size(), 0.0);
        for (std::size_t a = 0; a < actk.size(); ++a) {
            const std::vector<double>& cv = terms.modes[terms.imm_idx[actk[a]]].cvec;
            for (std::size_t i = 0; i < nc; ++i) Cf(i, a) = cv[idx[i]];
        }
        Matrix<double> G0(actk.size(), B.size(), 0.0);
        for (std::size_t jb = 0; jb < B.size(); ++jb) {
            std::vector<std::size_t> grp;
            double wsum = 0.0;
            for (std::size_t q = 0; q < actk.size(); ++q)
                if (imm.bind[actk[q]] == static_cast<std::ptrdiff_t>(B[jb])) {
                    grp.push_back(q);
                    wsum += terms.modes[terms.imm_idx[actk[q]]].weight;
                }
            if (grp.empty()) continue;
            const bool uniform = !(wsum > 0);
            for (std::size_t q : grp) {
                const double w =
                    uniform ? 1.0 : terms.modes[terms.imm_idx[actk[q]]].weight;
                G0(q, jb) = w / (uniform ? static_cast<double>(grp.size()) : wsum);
            }
        }
        Matrix<double> EB(B.size(), nc, 0.0);
        for (std::size_t jb = 0; jb < B.size(); ++jb) {
            const std::ptrdiff_t at = index_of(idx, B[jb]);
            if (at >= 0) EB(jb, static_cast<std::size_t>(at)) = 1.0;
        }
        const Matrix<double> Mb = matmul(matmul(EB, Cf), G0);
        const Matrix<double> corr = matmul(matmul(Cf, matmul(G0, pinv(Mb))), EB);
        for (std::size_t i = 0; i < nc; ++i)
            for (std::size_t j = 0; j < nc; ++j) T(i, j) -= corr(i, j);
    }
 
    std::vector<std::vector<double>> R;
    for (std::size_t c : active) {
        std::vector<double> row(nc, 0.0);
        for (std::size_t i = 0; i < nc; ++i) row[i] = con.A(c, idx[i]);
        R.push_back(row);
    }
    if (th != nullptr) {
        for (std::size_t j : terms.latch_mode) {
            std::vector<double> row(nc, 0.0);
            for (std::size_t z : terms.modes[j].zblk) {
                const std::ptrdiff_t at = index_of(idx, z);
                if (at >= 0) row[static_cast<std::size_t>(at)] = 1.0;
            }
            for (std::size_t q = 0; q < th->dslot[j].size(); ++q) {
                const std::ptrdiff_t at = index_of(idx, th->dslot[j][q]);
                if (at >= 0) row[static_cast<std::size_t>(at)] -= th->dval[j][q];
            }
            R.push_back(row);
        }
    }
    if (!R.empty()) {
        bool any = false;
        for (const std::vector<double>& row : R)
            for (double v : row)
                if (std::fabs(v) > 1e-14) any = true;
        if (any) {
            Matrix<double> Rm(R.size(), nc, 0.0);
            for (std::size_t i = 0; i < R.size(); ++i)
                for (std::size_t j = 0; j < nc; ++j) Rm(i, j) = R[i][j];
            Matrix<double> Rt(nc, R.size(), 0.0);
            for (std::size_t i = 0; i < R.size(); ++i)
                for (std::size_t j = 0; j < nc; ++j) Rt(j, i) = Rm(i, j);
            const Matrix<double> RRt = matmul(Rm, Rt);
            const Matrix<double> corr = matmul(matmul(Rt, pinv(RRt)), Rm);
            Matrix<double> P(nc, nc, 0.0);
            for (std::size_t i = 0; i < nc; ++i)
                for (std::size_t j = 0; j < nc; ++j) P(i, j) = (i == j ? 1.0 : 0.0) - corr(i, j);
            T = matmul(P, T);
        }
    }
    double worst = 0.0;
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < nc; ++j)
            worst = std::max(worst, std::fabs(T(i, j) - (i == j ? 1.0 : 0.0)));
    return worst > 1e-14;
}
 
/**
 * One Lyapunov solve, over the marking coordinates.
 *
 * THE DIFFUSION COUNTS THE STOCHASTIC EVENTS ONLY: an immediate flow is not a
 * Poisson stream with an intensity but the limit of an infinitely fast one whose
 * fluctuation is slaved, and its pinned coordinate is projected out.
 */
inline Matrix<double> sigma_of(const PetriTerms& terms, const Matrix<double>& A,
                               const std::vector<double>& r, const Matrix<double>* clampT) {
    const std::vector<std::size_t>& idx = terms.cov_idx;
    const std::size_t nc = idx.size(), ns = terms.stoch_col.size(), nfl = terms.flow_col.size();
    Matrix<double> Dc(nc, std::max<std::size_t>(ns, 1), 0.0);
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < ns; ++j) Dc(i, j) = terms.D(idx[i], terms.stoch_col[j]);
    Matrix<double> Df(nc, nfl, 0.0);
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < nfl; ++j) Df(i, j) = terms.D(idx[i], terms.flow_col[j]);
    Matrix<double> Am(nc, nc, 0.0);
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < nc; ++j) Am(i, j) = A(idx[i], idx[j]);
    if (clampT != nullptr) {
        // BOTH the jump directions and the generator are reduced: reducing Dc
        // alone would fix the subspace but leave the generator's orthogonal
        // component on it, which is not the reduced dynamics when the reduction
        // is oblique.
        Dc = matmul(*clampT, Dc);
        if (nfl > 0) Df = matmul(*clampT, Df);
        Am = matmul(*clampT, Am);
    }
    Matrix<double> Q(nc, nc, 0.0);
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < nc; ++j) {
            double v = 0.0;
            for (std::size_t e = 0; e < ns; ++e)
                v += Dc(i, e) * r[terms.stoch_col[e]] * Dc(j, e);
            Q(i, j) = v;
        }
    // A FLUID LEVEL FLUCTUATES THROUGH THE DRIFT, NOT THROUGH A JUMP OF ITS OWN: its flows carry no noise, so
    // range(Dc) alone would leave the level out of the reachable subspace and report it as deterministic. The
    // flow directions are added to the basis, which is where the marking's fluctuation is carried to.
    Matrix<double> basis(nc, Dc.cols() + nfl, 0.0);
    for (std::size_t i = 0; i < nc; ++i) {
        for (std::size_t j = 0; j < Dc.cols(); ++j) basis(i, j) = Dc(i, j);
        for (std::size_t j = 0; j < nfl; ++j) basis(i, Dc.cols() + j) = Df(i, j);
    }
    FluidLyapunovInfo info;
    const Matrix<double> Sc = fluid_lyapunov(Am, Q, basis, info);
    Matrix<double> Sigma(terms.nstate, terms.nstate, 0.0);
    for (std::size_t i = 0; i < nc; ++i)
        for (std::size_t j = 0; j < nc; ++j) Sigma(idx[i], idx[j]) = Sc(i, j);
    return Sigma;
}
 
/**
 * The multipliers the binding rows put on the two legs of the drift. A capacity row throttles the DEPOSIT leg of
 * the coordinates it covers, an empty-level row the WITHDRAWAL leg, each by the fraction zeta that keeps its
 * place on the bound. Port of `local_gains`.
 */
inline void gains(const PetriConstraints& con, const std::vector<std::size_t>& active,
                  const std::vector<double>& zeta, std::size_t n, std::vector<double>& gP,
                  std::vector<double>& gN) {
    gP.assign(n, 1.0);
    gN.assign(n, 1.0);
    for (std::size_t k = 0; k < active.size(); ++k) {
        const std::size_t c = active[k];
        const double z = k < zeta.size() ? zeta[k] : 1.0;
        const bool low = c < con.lower.size() && con.lower[c];
        for (std::size_t s = 0; s < n; ++s)
            if (con.cover[c][s]) (low ? gN[s] : gP[s]) *= z;
    }
}
 
/** The stacked residual, the rates it was evaluated at, and the covariance. */
struct Residual {
    bool ok = false;
    std::vector<double> G, r;
    PetriTheta th;
    Matrix<double> Sigma;
};
 
/**
 * The coupled algebraic system, stacked.
 *
 * `ok` is false when the closure cannot be evaluated at this iterate, so the line
 * search can back off; the first evaluation of a pass runs with `quiet=false`,
 * where a genuine failure surfaces.
 */
inline Residual residual(const std::vector<double>& u, const Ctx& ctx, bool quiet) {
    const PetriTerms& terms = *ctx.terms;
    const PetriImmediate& imm = *ctx.imm;
    const std::size_t n = terms.nstate;
    const Unpacked up = unpack(u, terms, imm, ctx.active.size());
 
    Residual res;
    res.th = petri_theta(terms, up.x, up.s2);
    res.r = petri_rates(terms, up.x, up.phi, up.mu, res.th);
 
    // the deposit gate of every binding capacity, and the withdrawal gate of every binding empty level, as
    // products of fractions
    std::vector<double> gP, gN;
    gains(*ctx.con, ctx.active, up.zeta, n, gP, gN);
    std::vector<double> drift(n, 0.0);
    for (std::size_t s = 0; s < n; ++s) {
        double neg = 0.0, pos = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) {
            neg += ctx.Dn(s, e) * res.r[e];
            pos += ctx.Dp(s, e) * res.r[e];
        }
        drift[s] = gN[s] * neg + gP[s] * pos;
    }
 
    try {
        const Matrix<double> A = petri_jacobian(terms, res.th);
        Matrix<double> T;
        const bool reduced =
            clamp_tangent(terms, imm, *ctx.con, ctx.active, &res.th, T);
        res.Sigma = sigma_of(terms, A, res.r, reduced ? &T : nullptr);
    } catch (const Error&) {
        if (!quiet) throw;
        res.ok = false;
        return res;
    }
 
    std::vector<double> G;
    G.reserve(n + ctx.C.rows() + terms.npair + imm.rows.size() + terms.latch_mode.size() +
              ctx.active.size());
    for (double d : drift) G.push_back(d);
    for (std::size_t r = 0; r < ctx.C.rows(); ++r) {
        double v = 0.0;
        for (std::size_t s = 0; s < n; ++s) v += ctx.C(r, s) * up.x[s];
        G.push_back(v - ctx.N[r]);
    }
    for (std::size_t p = 0; p < terms.npair; ++p)
        G.push_back(up.s2[p] - res.Sigma(terms.cov_pairs[p].first, terms.cov_pairs[p].second));
    for (const PetriImmediate::Row& row : imm.rows) {
        if (row.kind == PetriImmediate::PIN) G.push_back(up.x[row.a]);
        else if (row.kind == PetriImmediate::RATIO)
            G.push_back(up.phi[row.a] * row.wb - up.phi[row.b] * row.wa);
        else G.push_back(up.phi[row.a]);
    }
    for (std::size_t j : terms.latch_mode) {
        double s = 0.0;
        for (std::size_t z : terms.modes[j].zblk) s += up.x[z];
        G.push_back(s - res.th.theta[j]);
    }
    for (std::size_t k = 0; k < ctx.active.size(); ++k) {
        double v = 0.0;
        for (std::size_t s = 0; s < n; ++s) v += ctx.con->A(ctx.active[k], s) * up.x[s];
        G.push_back(v - ctx.con->b[ctx.active[k]]);
    }
    res.G = G;
    res.ok = true;
    return res;
}
 
inline double inf_norm(const std::vector<double>& v) {
    double m = 0.0;
    for (double d : v) m = std::max(m, std::fabs(d));
    return m;
}
 
/**
 * An iterate projected onto its feasible box, the lower bound only.
 *
 * The bound is a VECTOR, one entry per unknown, never a count: a scalar "nfree"
 * is indistinguishable from a one-unknown bound vector, which is how the MATLAB
 * twin crashed on the simplest net in the tree.
 */
inline std::vector<double> project(const std::vector<double>& u, const std::vector<double>& lb) {
    if (lb.empty()) return u;
    if (lb.size() != u.size())
        throw InputError("fluid_petri: the bound vector has " + std::to_string(lb.size()) +
                         " entries for " + std::to_string(u.size()) + " unknowns");
    std::vector<double> w = u;
    for (std::size_t i = 0; i < w.size(); ++i)
        if (std::isfinite(lb[i])) w[i] = std::max(lb[i], w[i]);
    return w;
}
 
/** Forward-difference Jacobian of the residual. */
inline Matrix<double> fdjac(const Ctx& ctx, const std::vector<double>& u,
                            const std::vector<double>& G) {
    const std::size_t n = u.size(), m = G.size();
    Matrix<double> J(std::max<std::size_t>(m, 1), std::max<std::size_t>(n, 1), 0.0);
    for (std::size_t k = 0; k < n; ++k) {
        const double h = 1e-7 * std::max(1.0, std::fabs(u[k]));
        std::vector<double> up = u;
        up[k] += h;
        Residual rp = residual(up, ctx, true);
        if (rp.ok) {
            for (std::size_t i = 0; i < m; ++i) J(i, k) = (rp.G[i] - G[i]) / h;
            continue;
        }
        up[k] = u[k] - h;
        Residual rm = residual(up, ctx, true);
        if (!rm.ok) continue;
        for (std::size_t i = 0; i < m; ++i) J(i, k) = (G[i] - rm.G[i]) / h;
    }
    return J;
}
 
struct NewtonResult {
    std::vector<double> u;
    std::size_t iterations = 0;
    bool converged = false;
    double resnorm = std::numeric_limits<double>::infinity();
};
 
/** Damped projected Newton with a finite-difference Jacobian. */
inline NewtonResult newton(const Ctx& ctx, const std::vector<double>& u0, double tol,
                           std::size_t maxit, const std::vector<double>& lb) {
    NewtonResult out;
    std::vector<double> u = project(u0, lb);
    Residual res = residual(u, ctx, false);
    if (!res.ok) {
        out.u = u;
        return out;
    }
    std::vector<double> G = res.G;
    double resnorm = inf_norm(G);
    std::size_t it = 0;
    for (it = 1; it <= maxit; ++it) {
        if (resnorm <= tol) {
            out.u = u;
            out.iterations = it - 1;
            out.converged = true;
            out.resnorm = resnorm;
            return out;
        }
        const Matrix<double> J = fdjac(ctx, u, G);
        Matrix<double> rhs(G.size(), 1, 0.0);
        for (std::size_t i = 0; i < G.size(); ++i) rhs(i, 0) = -G[i];
        const Matrix<double> du = matmul(pinv(J), rhs);
        double lam = 1.0;
        bool improved = false;
        for (int b = 0; b < 30; ++b) {
            std::vector<double> un(u.size(), 0.0);
            for (std::size_t i = 0; i < u.size(); ++i) un[i] = u[i] + lam * du(i, 0);
            un = project(un, lb);
            Residual rn = residual(un, ctx, true);
            if (rn.ok) {
                const double v = inf_norm(rn.G);
                if (v < resnorm) {
                    u = un;
                    G = rn.G;
                    resnorm = v;
                    improved = true;
                    break;
                }
            }
            lam *= 0.5;
        }
        if (!improved) break;
    }
    out.u = u;
    out.iterations = it;
    out.converged = resnorm <= tol;
    out.resnorm = resnorm;
    return out;
}
 
/**
 * The first-order drift: the same rates at zero variance, with an immediate mode
 * firing at LAM times its enabling degree and its firing weight, and the server
 * latch relaxed at LAM towards the enabling degree instead of solved. Both are
 * approximations of an algebraic constraint by a fast reaction, and both are
 * confined to the seed.
 */
inline std::vector<double> seed_drift(const PetriTerms& terms, const std::vector<double>& xin,
                                      double lam) {
    std::vector<double> x = xin;
    for (std::size_t i = 0; i < x.size(); ++i) x[i] = std::max(0.0, x[i]);
    // The inhibitor gate is smoothed here (see petri_theta): a hard step chatters.
    const PetriTheta th = petri_theta(
        terms, x, std::vector<double>(std::max<std::size_t>(terms.npair, 1), 0.0), petri_coarse_tol());
    std::vector<double> phi(terms.imm_idx.size(), 0.0);
    for (std::size_t k = 0; k < phi.size(); ++k) {
        const std::size_t j = terms.imm_idx[k];
        phi[k] = lam * terms.modes[j].weight * th.theta[j];
    }
    std::vector<double> mu(terms.latch_mode.size(), 0.0);
    for (std::size_t q = 0; q < mu.size(); ++q) {
        const std::size_t j = terms.latch_mode[q];
        double s = 0.0;
        for (std::size_t z : terms.modes[j].zblk) s += x[z];
        mu[q] = lam * (th.theta[j] - s);
    }
    const std::vector<double> r = petri_rates(terms, x, phi, mu, th);
    std::vector<double> d(terms.nstate, 0.0);
    for (std::size_t s = 0; s < terms.nstate; ++s) {
        double v = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) v += terms.D(s, e) * r[e];
        d[s] = v;
    }
    // A fluid level stops at its bounds: the flow that would carry it past the empty level or the capacity is
    // the part the binding multiplier would cut. The bound test reads the UNCLIPPED state, as the reference does.
    for (std::size_t q = 0; q < terms.nf; ++q) {
        const std::size_t slot = terms.fluid_slot[q];
        if ((xin[slot] <= 0.0 && d[slot] < 0.0) || (xin[slot] >= terms.fluid_cap[q] && d[slot] > 0.0))
            d[slot] = 0.0;
    }
    return d;
}
 
/**
 * The bound rows of the fluid places that the seed ended on, and the refusals of the two passive levels no
 * fixed point describes. Port of `local_fluid_seed_active`.
 *
 * A level that NOTHING READS (no guard, no fluid rate dependence) has a drift column of zeros, so the algebraic
 * solve cannot place it: it sits on a bound when its drift points out of the interval, and otherwise it has no
 * stationary value at all -- it grows without limit under a positive drift and an infinite capacity, and it
 * wanders with no restoring force under a zero drift in the interior. Both are refused by name rather than
 * answered with the seed's end.
 */
inline std::vector<std::size_t> fluid_seed_active(const PetriTerms& terms, const PetriConstraints& con,
                                                  const std::vector<double>& x) {
    std::vector<std::size_t> active;
    if (terms.nf == 0) return active;
    std::vector<double> xp = x;
    for (std::size_t i = 0; i < xp.size(); ++i) xp[i] = std::max(0.0, xp[i]);
    const PetriTheta th = petri_theta(
        terms, xp, std::vector<double>(std::max<std::size_t>(terms.npair, 1), 0.0), petri_coarse_tol());
    const std::vector<double> r = petri_rates(terms, xp, std::vector<double>(terms.imm_idx.size(), 0.0),
                                              std::vector<double>(terms.latch_mode.size(), 0.0), th);
    // The reference evaluates the seed drift first and then overwrites every fluid entry of it with this free
    // drift, so only the free drift is read.
    double maxbase = 0.0;
    for (double b : terms.rate_base) maxbase = std::max(maxbase, std::fabs(b));
    for (std::size_t q = 0; q < terms.nf; ++q) {
        const std::size_t slot = terms.fluid_slot[q];
        double dfree = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) dfree += terms.D(slot, e) * r[e];
        bool read = false;
        for (const PetriMode& md : terms.modes) {
            if (std::find(md.guard_slot.begin(), md.guard_slot.end(), slot) != md.guard_slot.end()) read = true;
            if (std::find(md.fdep_slot.begin(), md.fdep_slot.end(), slot) != md.fdep_slot.end()) read = true;
        }
        const double cap = terms.fluid_cap[q];
        const double scale = std::max(1.0, std::max(std::fabs(dfree), maxbase));
        const double tolx = petri_coarse_tol() * std::max(1.0, std::min(cap, std::max(std::fabs(x[slot]), 1.0)));
        const double told = petri_coarse_tol() * scale;
        std::ptrdiff_t low_row = -1, up_row = -1;
        for (std::size_t c = 0; c < con.b.size(); ++c) {
            if (!(con.fluid[c] && con.cover[c][slot])) continue;
            if (con.lower[c]) low_row = static_cast<std::ptrdiff_t>(c);
            else up_row = static_cast<std::ptrdiff_t>(c);
        }
        if (x[slot] <= tolx && dfree <= told) {
            if (low_row >= 0) active.push_back(static_cast<std::size_t>(low_row));
        } else if (std::isfinite(cap) && x[slot] >= cap - tolx && dfree >= -told) {
            if (up_row >= 0) active.push_back(static_cast<std::size_t>(up_row));
        } else if (!read) {
            char buf[64];
            std::snprintf(buf, sizeof(buf), "%g", dfree);
            if (!std::isfinite(cap) && dfree > told)
                throw NumericError(ErrorKind{"Divergent"},
                    "solver_fluid_petri: Fluid place " + terms.fluid_name[q] + " has a positive mean drift (" +
                    std::string(buf) + ") and no capacity, so its level grows without bound and has no steady "
                    "state. Give it a capacity with setFluidCapacity.");
            if (std::fabs(dfree) <= told)
                throw NumericError(ErrorKind{"Degenerate"},
                    "solver_fluid_petri: Fluid place " + terms.fluid_name[q] + " has a zero mean drift and "
                    "nothing reads its level, so its stationary level is undetermined (the level is null "
                    "recurrent).");
        }
    }
    return active;
}
 
/**
 * The reported level of every fluid place. Port of `local_fluid_levels`.
 *
 * Under `reflect` (the default) it is the fixed point: a level on a bound is the bound. Under `diffusion` a
 * level on a bound is replaced by the stationary mean of a Brownian motion reflected on [0, B], with the drift
 * mu the level would have without the bound and the variance rate sigma2 of its net inflow. The inflow
 * fluctuates through its impulses (the jumps delta) and through the marking that drives its flows (gradient g),
 * and the long-run variance rate of that integrated linear-noise process is
 *
 *     sigma2 = (delta - g W) diag(r) (delta - g W)',   W = V (V' A V)^-1 V' Dc
 *
 * on the span V of the stochastic jumps Dc of the marking. With kappa = 2 mu / sigma2 the reflected law is
 * exponential on [0, B]:
 *
 *     E = B / (1 - exp(-kappa B)) - 1/kappa,   Var = 1/kappa^2 - B^2 / (4 sinh^2(kappa B / 2))
 *
 * and E = -1/kappa, Var = 1/kappa^2 when B = Inf. The refinement is reported, not fed back into the fixed point.
 */
inline PetriReport::Fluid fluid_levels(const PetriTerms& terms, const PetriConstraints& con,
                                       const PetriImmediate& imm, const std::vector<std::size_t>& active,
                                       const std::vector<double>& x, const std::vector<double>& s2,
                                       const std::vector<double>& phi, const std::vector<double>& r,
                                       const Matrix<double>& Sigma, const std::string& boundary) {
    std::string mode = boundary.empty() ? std::string("reflect") : boundary;
    for (char& ch : mode) ch = static_cast<char>(std::tolower(static_cast<unsigned char>(ch)));
    if (mode != "reflect" && mode != "diffusion")
        throw InputError(ErrorKind{"ConfigInvalid"},
                         "solver_fluid_petri: options.config.fluid_boundary must be 'reflect' or 'diffusion', "
                         "got '" + mode + "'.");
    const std::size_t nf = terms.nf;
    PetriReport::Fluid fl;
    fl.names = terms.fluid_name;
    fl.level.assign(nf, 0.0);
    fl.level_var.assign(nf, 0.0);
    fl.drift.assign(nf, 0.0);
    fl.sigma2.assign(nf, 0.0);
    fl.bound.assign(nf, 0);
    fl.boundary = mode;
    if (nf == 0) return fl;
    bool any_bound = false;
    for (std::size_t q = 0; q < nf; ++q) {
        const std::size_t slot = terms.fluid_slot[q];
        fl.level[q] = x[slot];
        fl.level_var[q] = std::max(0.0, Sigma(slot, slot));
        double d = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) d += terms.D(slot, e) * r[e];
        fl.drift[q] = d;
        for (std::size_t c : active)
            if (con.fluid[c] && con.cover[c][slot]) fl.bound[q] = con.lower[c] ? -1 : 1;
        if (fl.bound[q] != 0) any_bound = true;
    }
    if (mode != "diffusion" || !any_bound) return fl;
 
    // the marking that drives the flows, reduced as the covariance is, but with the fluid bounds left out: they
    // are what the refinement replaces
    const PetriTheta th = petri_theta(terms, x, s2);
    const Matrix<double> A = petri_jacobian(terms, th);
    (void)phi;
    std::vector<std::size_t> tokact;
    for (std::size_t c : active)
        if (!con.fluid[c]) tokact.push_back(c);
    Matrix<double> T;
    clamp_tangent(terms, imm, con, tokact, &th, T);
    std::vector<bool> isfl(terms.nstate, false);
    for (std::size_t q = 0; q < nf; ++q) isfl[terms.fluid_slot[q]] = true;
    // tok: positions in cov_idx of the coordinates that are not fluid levels (cov_idx is the identity here)
    std::vector<std::size_t> tokpos, tokslot;
    for (std::size_t i = 0; i < terms.cov_idx.size(); ++i)
        if (!isfl[terms.cov_idx[i]]) {
            tokpos.push_back(i);
            tokslot.push_back(terms.cov_idx[i]);
        }
    const std::size_t nt = tokslot.size(), ns = terms.stoch_col.size();
    Matrix<double> Ttok(nt, nt, 0.0), Dtok(nt, ns, 0.0), Atok(nt, nt, 0.0);
    for (std::size_t i = 0; i < nt; ++i) {
        for (std::size_t j = 0; j < nt; ++j) {
            Ttok(i, j) = T(tokpos[i], tokpos[j]);
            Atok(i, j) = A(tokslot[i], tokslot[j]);
        }
        for (std::size_t e = 0; e < ns; ++e) Dtok(i, e) = terms.D(tokslot[i], terms.stoch_col[e]);
    }
    const Matrix<double> Dc = (nt > 0 && ns > 0) ? matmul(Ttok, Dtok) : Matrix<double>(nt, ns, 0.0);
    const Matrix<double> Ar = nt > 0 ? matmul(Ttok, Atok) : Matrix<double>(0, 0, 0.0);
    const Matrix<double> V = (nt > 0 && ns > 0) ? fluid::detail::fluid_orth(Dc) : Matrix<double>(nt, 0, 0.0);
    Matrix<double> W(nt, ns, 0.0);
    if (V.cols() > 0) {
        const Matrix<double> Vt = V.transpose();
        const Matrix<double> VAV = matmul(matmul(Vt, Ar), V);
        W = matmul(V, matmul(inverse(VAV), matmul(Vt, Dc)));
    }
    for (std::size_t q = 0; q < nf; ++q) {
        if (fl.bound[q] == 0) continue;
        const std::size_t slot = terms.fluid_slot[q];
        std::vector<double> h(ns, 0.0);
        for (std::size_t e = 0; e < ns; ++e) h[e] = terms.D(slot, terms.stoch_col[e]);
        if (V.cols() > 0)
            for (std::size_t e = 0; e < ns; ++e) {
                double gw = 0.0;
                for (std::size_t i = 0; i < nt; ++i) gw += A(slot, tokslot[i]) * W(i, e);
                h[e] -= gw;
            }
        double sig2 = 0.0;
        for (std::size_t e = 0; e < ns; ++e) sig2 += h[e] * h[e] * r[terms.stoch_col[e]];
        fl.sigma2[q] = sig2;
        const double mu = fl.drift[q];
        const double B = terms.fluid_cap[q];
        if (!(sig2 > 1e-14)) continue;  // GlobalConstants.Zero
        const double kappa = 2.0 * mu / sig2;
        if (!std::isfinite(B)) {
            if (kappa >= 0) continue;
            fl.level[q] = -1.0 / kappa;
            fl.level_var[q] = 1.0 / (kappa * kappa);
        } else if (std::fabs(kappa * B) < 1e-8) {
            fl.level[q] = B / 2.0;
            fl.level_var[q] = B * B / 12.0;
        } else {
            fl.level[q] = B / (-std::expm1(-kappa * B)) - 1.0 / kappa;
            const double sh = std::sinh(kappa * B / 2.0);
            fl.level_var[q] = std::max(0.0, 1.0 / (kappa * kappa) - B * B / (4.0 * sh * sh));
        }
    }
    return fl;
}
 
// ============================ the transient ==============================
 
/**
 * The immediate flows the SEED trajectory carried, point by point: the approximation the seed integrated, and
 * therefore the one its throughput table has to be read at.
 */
inline std::vector<std::vector<double>> seed_flows(const PetriTerms& terms,
                                                   const std::vector<std::vector<double>>& xt, double lam) {
    const std::size_t ni = terms.imm_idx.size();
    std::vector<std::vector<double>> phit(xt.size(), std::vector<double>(ni, 0.0));
    if (ni == 0) return phit;
    for (std::size_t a = 0; a < xt.size(); ++a) {
        const PetriTheta th = petri_theta(
            terms, xt[a], std::vector<double>(std::max<std::size_t>(terms.npair, 1), 0.0), petri_coarse_tol());
        for (std::size_t k = 0; k < ni; ++k) {
            const std::size_t j = terms.imm_idx[k];
            phit[a][k] = lam * terms.modes[j].weight * th.theta[j];
        }
    }
    return phit;
}
 
/**
 * MATLAB's `A\b`, which the collapse is defined by. A rectangular system takes the BASIC solution of a
 * column-pivoted Householder QR (LAPACK xGEQP3, first column of largest norm on a tie): with two immediate modes
 * draining one place it routes the whole mass down the first, where a minimum-norm solve would split it. A
 * singular square system returns false, the non-finite answer the caller already falls back on.
 */
inline bool petri_mldivide(std::vector<std::vector<double>> A, std::vector<double> b, std::vector<double>& xs) {
    const std::size_t m = A.size();
    const std::size_t nn = m == 0 ? 0 : A[0].size();
    std::vector<std::size_t> perm(nn);
    for (std::size_t j = 0; j < nn; ++j) perm[j] = j;
    const std::size_t kmax = std::min(m, nn);
    std::vector<double> diagR(kmax, 0.0);
    for (std::size_t k = 0; k < kmax; ++k) {
        std::size_t p = k;
        double best = -1.0;
        for (std::size_t j = k; j < nn; ++j) {
            double nrm = 0.0;
            for (std::size_t i = k; i < m; ++i) nrm += A[i][j] * A[i][j];
            if (nrm > best) {
                best = nrm;
                p = j;
            }
        }
        if (p != k) {
            for (std::size_t i = 0; i < m; ++i) std::swap(A[i][k], A[i][p]);
            std::swap(perm[k], perm[p]);
        }
        double alpha = std::sqrt(best);
        if (alpha == 0.0) {
            diagR[k] = 0.0;
            continue;
        }
        if (A[k][k] > 0) alpha = -alpha;
        std::vector<double> v(m, 0.0);
        for (std::size_t i = k; i < m; ++i) v[i] = A[i][k];
        v[k] -= alpha;
        double vv = 0.0;
        for (std::size_t i = k; i < m; ++i) vv += v[i] * v[i];
        if (vv > 0) {
            for (std::size_t j = k; j < nn; ++j) {
                double d = 0.0;
                for (std::size_t i = k; i < m; ++i) d += v[i] * A[i][j];
                d = 2.0 * d / vv;
                for (std::size_t i = k; i < m; ++i) A[i][j] -= d * v[i];
            }
            double d = 0.0;
            for (std::size_t i = k; i < m; ++i) d += v[i] * b[i];
            d = 2.0 * d / vv;
            for (std::size_t i = k; i < m; ++i) b[i] -= d * v[i];
        }
        diagR[k] = A[k][k];
    }
    const double tol = kmax > 0 ? static_cast<double>(std::max(m, nn)) * std::numeric_limits<double>::epsilon() *
                                      std::fabs(diagR[0])
                                : 0.0;
    std::size_t rank = 0;
    while (rank < kmax && std::fabs(diagR[rank]) > tol) ++rank;
    if (m == nn && rank < nn) return false;
    std::vector<double> z(nn, 0.0);
    for (std::size_t ii = rank; ii-- > 0;) {
        double acc = b[ii];
        for (std::size_t j = ii + 1; j < rank; ++j) acc -= A[ii][j] * z[j];
        z[ii] = acc / A[ii][ii];
    }
    xs.assign(nn, 0.0);
    for (std::size_t j = 0; j < nn; ++j) xs[perm[j]] = z[j];
    return true;
}
 
/**
 * The fluid vanishing-marking collapse: move the initial marking along the immediate incidence columns until
 * every pinned place is empty. Every conserved quantity survives it, since a conserved direction annihilates
 * those columns. Port of the collapse in `local_transient`.
 */
inline std::vector<double> collapse(const PetriTerms& terms, const PetriImmediate& imm,
                                    const std::vector<double>& x0) {
    std::vector<double> x = x0;
    std::vector<std::size_t> act;
    for (std::size_t k = 0; k < imm.n; ++k)
        if (imm.active[k]) act.push_back(k);
    const std::vector<std::size_t>& pins = imm.pins;
    if (act.empty() || pins.empty()) return x;
    const std::size_t n = terms.nstate;
    std::vector<std::vector<double>> Aimm(pins.size(), std::vector<double>(act.size(), 0.0));
    std::vector<double> b(pins.size(), 0.0);
    for (std::size_t i = 0; i < pins.size(); ++i) {
        for (std::size_t a = 0; a < act.size(); ++a)
            Aimm[i][a] = terms.modes[terms.imm_idx[act[a]]].cvec[pins[i]];
        b[i] = -x0[pins[i]];
    }
    std::vector<double> s;
    if (!petri_mldivide(Aimm, b, s)) return x;
    for (double& v : s) {
        if (!std::isfinite(v)) return x;
        v = std::max(0.0, v);
    }
    for (std::size_t r = 0; r < n; ++r) {
        double v = 0.0;
        for (std::size_t a = 0; a < s.size(); ++a) v += terms.modes[terms.imm_idx[act[a]]].cvec[r] * s[a];
        x[r] = x0[r] + v;
    }
    for (std::size_t p : pins) x[p] = 0.0;
    for (std::size_t r = 0; r < terms.nm; ++r) x[r] = std::max(0.0, x[r]);
    // the collapse moved the marking, so the servers a multi-phase mode has latched moved with it
    const PetriTheta th =
        petri_theta(terms, x, std::vector<double>(std::max<std::size_t>(terms.npair, 1), 0.0));
    for (std::size_t j : terms.latch_mode) {
        const PetriMode& md = terms.modes[j];
        for (std::size_t h = 0; h < md.zblk.size(); ++h) x[md.zblk[h]] = th.theta[j] * md.pie[h];
    }
    return x;
}
 
/** The right-hand side of one transient segment, in the layout `petri_transient` documents. */
inline std::vector<double> dae_rhs(const double* z, const Ctx& ctx, bool withcov, std::size_t nc) {
    const PetriTerms& terms = *ctx.terms;
    const PetriImmediate& imm = *ctx.imm;
    const PetriConstraints& con = *ctx.con;
    const std::size_t n = terms.nstate, nimm = imm.n, nact = ctx.active.size(), nlatch = terms.latch_mode.size();
    const std::size_t nvar = n + nimm + nlatch + nact;
    const std::vector<double> x(z, z + n), phi(z + n, z + n + nimm), mu(z + n + nimm, z + n + nimm + nlatch),
        zeta(z + n + nimm + nlatch, z + nvar);
    const PetriTheta th = petri_theta(
        terms, x, std::vector<double>(std::max<std::size_t>(terms.npair, 1), 0.0), 0.0, petri_coarse_tol());
    const std::vector<double> r = petri_rates(terms, x, phi, mu, th);
    std::vector<double> gP, gN;
    gains(con, ctx.active, zeta, n, gP, gN);
    std::vector<double> drift(n, 0.0);
    for (std::size_t s = 0; s < n; ++s) {
        double neg = 0.0, pos = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) {
            neg += ctx.Dn(s, e) * r[e];
            pos += ctx.Dp(s, e) * r[e];
        }
        drift[s] = gN[s] * neg + gP[s] * pos;
    }
    std::vector<double> dz(nvar + (withcov ? nc * nc : 0), 0.0);
    for (std::size_t s = 0; s < n; ++s) dz[s] = drift[s];
    for (std::size_t p : imm.pins) dz[p] = x[p];  // the algebraic level pin
    std::size_t at = n;
    for (const PetriImmediate::Row& row : imm.rows) {
        if (row.kind == PetriImmediate::PIN) dz[at++] = drift[row.a];
        else if (row.kind == PetriImmediate::RATIO) dz[at++] = phi[row.a] * row.wb - phi[row.b] * row.wa;
        else dz[at++] = phi[row.a];
    }
    // the DIFFERENTIATED latch: the running-server count tracks the enabling degree, so their rates agree
    for (std::size_t j : terms.latch_mode) {
        double acc = 0.0;
        for (std::size_t zc : terms.modes[j].zblk) acc += drift[zc];
        for (std::size_t q = 0; q < th.dslot[j].size(); ++q) acc -= th.dval[j][q] * drift[th.dslot[j][q]];
        dz[at++] = acc;
    }
    for (std::size_t k = 0; k < nact; ++k) {
        double v = 0.0;
        for (std::size_t s = 0; s < n; ++s) v += con.A(ctx.active[k], s) * drift[s];
        dz[at++] = v;
    }
    if (withcov) {
        const std::vector<std::size_t>& idx = terms.cov_idx;
        const std::vector<std::size_t>& sc = terms.stoch_col;
        const Matrix<double> A = petri_jacobian(terms, th);
        Matrix<double> Dc(nc, std::max<std::size_t>(sc.size(), 1), 0.0), Am(nc, nc, 0.0);
        for (std::size_t i = 0; i < nc; ++i) {
            for (std::size_t e = 0; e < sc.size(); ++e) Dc(i, e) = terms.D(idx[i], sc[e]);
            for (std::size_t j = 0; j < nc; ++j) Am(i, j) = A(idx[i], idx[j]);
        }
        Matrix<double> T;
        if (clamp_tangent(terms, imm, con, ctx.active, &th, T)) {
            Dc = matmul(T, Dc);
            Am = matmul(T, Am);
        }
        Matrix<double> S(nc, nc, 0.0);
        for (std::size_t i = 0; i < nc; ++i)
            for (std::size_t j = 0; j < nc; ++j) S(i, j) = z[nvar + i * nc + j];
        const Matrix<double> AS = matmul(Am, S), SAt = matmul(S, Am.transpose());
        for (std::size_t i = 0; i < nc; ++i)
            for (std::size_t j = 0; j < nc; ++j) {
                double q = 0.0;
                for (std::size_t e = 0; e < sc.size(); ++e) q += Dc(i, e) * r[sc[e]] * Dc(j, e);
                dz[nvar + i * nc + j] = AS(i, j) + SAt(i, j) + q;
            }
    }
    return dz;
}
 
/**
 * Where a segment ends: a flow that would go negative, a capacity that starts or stops binding, or a marking
 * coordinate that would go negative. A crossing is a value going from positive to non-positive.
 */
inline std::vector<double> event_values(const double* z, const Ctx& ctx) {
    const PetriTerms& terms = *ctx.terms;
    const PetriImmediate& imm = *ctx.imm;
    const PetriConstraints& con = *ctx.con;
    const std::size_t n = terms.nstate, nimm = imm.n, nlatch = terms.latch_mode.size(), ncon = con.b.size();
    std::vector<double> val(nimm + ncon + terms.nm, 0.0);
    for (std::size_t k = 0; k < nimm; ++k) val[k] = imm.active[k] ? z[n + k] : 1.0;
    for (std::size_t c = 0; c < ncon; ++c) {
        const std::ptrdiff_t k = index_of(ctx.active, c);
        if (k < 0) {
            double ax = 0.0;
            for (std::size_t s = 0; s < n; ++s) ax += con.A(c, s) * z[s];
            val[nimm + c] = con.b[c] - ax;
        } else {
            val[nimm + c] = 1.0 - z[n + nimm + nlatch + static_cast<std::size_t>(k)];
        }
    }
    for (std::size_t s = 0; s < terms.nm; ++s)
        val[nimm + ncon + s] = index_of(imm.pins, s) >= 0 ? 1.0 : z[s] + petri_fine_tol();
    return val;
}
 
/** One segment's integration state, reached by the RODAS callbacks through `petri_tran_active`. */
struct PetriTranSeg {
    const Ctx* ctx = nullptr;
    bool withcov = false;
    std::size_t nc = 0, nz = 0;
    std::vector<double> diag, gprev;
    const std::vector<double>* req = nullptr;
    std::size_t jreq = 0;
    std::vector<double> segT;
    std::vector<std::vector<double>> segZ;
    int hit = -1;
};
inline PetriTranSeg*& petri_tran_active() {
    static PetriTranSeg* p = nullptr;
    return p;
}
 
inline int petri_tran_fcn(rodas_impl::integer*, rodas_impl::doublereal*, rodas_impl::doublereal* y,
                          rodas_impl::doublereal* f, rodas_impl::doublereal*, rodas_impl::integer*) {
    PetriTranSeg& S = *petri_tran_active();
    const std::vector<double> d = dae_rhs(y, *S.ctx, S.withcov, S.nc);
    for (std::size_t i = 0; i < S.nz; ++i) f[i] = d[i];
    return 0;
}
inline int petri_tran_jac(rodas_impl::integer*, rodas_impl::doublereal*, rodas_impl::doublereal*,
                          rodas_impl::doublereal*, rodas_impl::integer*, rodas_impl::doublereal*,
                          rodas_impl::integer*) { return 0; }
inline int petri_tran_dfx(rodas_impl::integer*, rodas_impl::doublereal*, rodas_impl::doublereal*,
                          rodas_impl::doublereal*, rodas_impl::doublereal*, rodas_impl::integer*) { return 0; }
/** The singular mass matrix: banded with MLMAS=MUMAS=0, i.e. its diagonal. */
inline int petri_tran_mas(rodas_impl::integer* n, rodas_impl::doublereal* am, rodas_impl::integer* lmas,
                          rodas_impl::doublereal*, rodas_impl::integer*) {
    PetriTranSeg& S = *petri_tran_active();
    for (int j = 0; j < *n; ++j) am[0 + j * (*lmas)] = S.diag[static_cast<std::size_t>(j)];
    return 0;
}
inline std::vector<double> petri_tran_dense(double tq, rodas_impl::doublereal* cont, rodas_impl::integer* lrc,
                                            std::size_t nz) {
    std::vector<double> y(nz, 0.0);
    for (std::size_t i = 0; i < nz; ++i) {
        rodas_impl::integer ii = static_cast<rodas_impl::integer>(i + 1);
        double t = tq;
        y[i] = rodas_impl::contro_(&ii, &t, cont, lrc);
    }
    return y;
}
/** The requested instants (`tranpoints`) the step (xold, xend) passed, read off its interpolant. */
inline void petri_tran_emit(PetriTranSeg& S, double xold, double xend, rodas_impl::doublereal* cont,
                            rodas_impl::integer* lrc) {
    const std::vector<double>& req = *S.req;
    while (S.jreq < req.size() && req[S.jreq] < xend) {
        const double tq = req[S.jreq];
        if (tq > xold && (S.segT.empty() || tq > S.segT.back())) {
            S.segT.push_back(tq);
            S.segZ.push_back(petri_tran_dense(tq, cont, lrc, S.nz));
        }
        ++S.jreq;
    }
}
/**
 * Every accepted step is recorded, the requested instants are added off the dense output, and a crossing is
 * located on that same interpolant. THE OVERSHOOTING STEP IS NOT REPORTED. The first call precedes any step
 * (NR = 1, XOLD = X = t0), where `cont` holds no coefficients yet.
 */
inline int petri_tran_solout(rodas_impl::integer* nr, rodas_impl::doublereal* xold, rodas_impl::doublereal* x,
                             rodas_impl::doublereal* y, rodas_impl::doublereal* cont, rodas_impl::integer* lrc,
                             rodas_impl::integer*, rodas_impl::doublereal*, rodas_impl::integer*,
                             rodas_impl::integer* irtrn) {
    PetriTranSeg& S = *petri_tran_active();
    const std::vector<double> gcur = event_values(y, *S.ctx);
    int cross = -1;
    for (std::size_t c = 0; c < gcur.size(); ++c)
        if (S.gprev[c] > 0 && gcur[c] <= 0) {
            cross = static_cast<int>(c);
            break;
        }
    S.gprev = gcur;
    if (cross < 0) {
        if (*nr > 1) petri_tran_emit(S, *xold, *x, cont, lrc);
        S.segT.push_back(*x);
        S.segZ.push_back(std::vector<double>(y, y + S.nz));
        return 0;
    }
    double lo = *xold, hi = *x;
    for (int bit = 0; bit < 60; ++bit) {
        const double mid = 0.5 * (lo + hi);
        const std::vector<double> ymid = petri_tran_dense(mid, cont, lrc, S.nz);
        if (event_values(ymid.data(), *S.ctx)[static_cast<std::size_t>(cross)] > 0) lo = mid;
        else hi = mid;
        if (hi - lo <= 1e-12 * std::max(1.0, std::fabs(hi))) break;
    }
    S.hit = cross;
    petri_tran_emit(S, *xold, hi, cont, lrc);
    S.segT.push_back(hi);
    S.segZ.push_back(petri_tran_dense(hi, cont, lrc, S.nz));
    // RODAS reads the stop through IRTRN, not through the return value
    if (irtrn) *irtrn = -1;
    return 0;
}
 
/** What `petri_transient` returns. */
struct PetriTransient {
    std::vector<double> t;
    std::vector<std::vector<double>> xt, phit;
    std::vector<Matrix<double>> Sigmat;
    std::vector<double> tvar;
    std::vector<std::pair<double, int>> switches;
};
 
/**
 * The transient closure as an index-1 DAE with a SINGULAR mass matrix, integrated by RODAS segment by segment
 * between active-set switches. Port of `local_transient` in `matlab/src/solvers/FLD/solver_fluid_petri.m`.
 *
 * THE INITIAL MARKING JUMPS: an immediate transition fires in zero time, so the trajectory starts at the
 * collapsed marking. Per segment the unknown is [x; phi; mu; zeta; vec(S)]: a free coordinate integrates its
 * drift; a pinned one carries both `x_b = 0` and the differentiated pin `drift_b = 0`, which determines the flow;
 * the conflict rule, the differentiated server latch and the differentiated binding rows are algebraic; and the
 * covariance `dS = A S + S A' + Q` rides along when it fits `dae_maxcov`.
 *
 * @param imm the steady-state immediate set, copied and moved here
 */
inline PetriTransient petri_transient(const PetriTerms& terms, const PetriConservation& cons,
                                      const PetriConstraints& con, PetriImmediate imm,
                                      const std::vector<std::size_t>& active_ss, const PetriOptions& opt,
                                      const std::vector<double>& phiss, const std::vector<double>& zetass,
                                      std::vector<std::string>& warnings) {
    using namespace rodas_impl;
    const std::size_t n = terms.nstate, nimm = imm.n, nlatch = terms.latch_mode.size(), ncon = con.b.size();
    const double tEnd = opt.timespan_end;
    const double tStart = std::isfinite(opt.timespan_start) ? opt.timespan_start : 0.0;
    const std::size_t maxcov = opt.dae_maxcov > 0 ? opt.dae_maxcov : 25;
    const std::size_t nc = terms.cov_idx.size();
    const bool withcov = nc > 0 && nc <= maxcov;
    if (!withcov && nc > maxcov)
        warnings.push_back("The transient covariance would add " + std::to_string(nc * nc) +
                           " differential states, above the limit of " + std::to_string(maxcov * maxcov) +
                           " set by options.config.dae_maxcov. Integrating the mean as a DAE with the variance "
                           "held at its stationary value.");
 
    const std::vector<double> x0 = collapse(terms, imm, terms.x0);
 
    // THE FLUID BOUNDS ARE RE-READ AT THE INITIAL LEVEL, not taken from the steady state: a bound the initial
    // level sits on starts binding unthrottled, and the events release it when it should.
    std::vector<std::size_t> active;
    std::vector<double> zeta;
    for (std::size_t k = 0; k < active_ss.size(); ++k) {
        const std::size_t c = active_ss[k];
        if (!con.fluid[c]) {
            active.push_back(c);
            zeta.push_back(k < zetass.size() ? zetass[k] : 1.0);
        }
    }
    for (std::size_t q = 0; q < terms.nf; ++q) {
        const std::size_t slot = terms.fluid_slot[q];
        for (std::size_t c = 0; c < ncon; ++c) {
            if (!con.fluid[c] || !con.cover[c][slot]) continue;
            if ((x0[slot] <= 0 && con.lower[c]) ||
                (x0[slot] > 0 && x0[slot] >= terms.fluid_cap[q] && !con.lower[c])) {
                active.push_back(c);
                zeta.push_back(1.0);
                break;
            }
        }
    }
 
    std::vector<double> req = opt.tranpoints;
    std::sort(req.begin(), req.end());
    PetriTransient out;
    double tcur = tStart;
    std::vector<double> xcur = x0, phi = phiss, mu(nlatch, 0.0), Scur(withcov ? nc * nc : 0, 0.0);
    const std::size_t segMax = std::max<std::size_t>(8, 4 * (nimm + ncon + 1));
    for (std::size_t seg = 0; seg < segMax; ++seg) {
        // THE ACTIVE SET IS RE-READ EVERY SEGMENT: a located crossing may have added or released a row.
        const Ctx ctx = context(terms, cons, con, imm, active);
        const std::size_t nact = ctx.active.size();
        const std::size_t nvar = n + nimm + nlatch + nact;
        const std::size_t nz = nvar + (withcov ? nc * nc : 0);
        PetriTranSeg S;
        S.ctx = &ctx;
        S.withcov = withcov;
        S.nc = nc;
        S.nz = nz;
        S.diag.assign(nz, 0.0);
        for (std::size_t s = 0; s < n; ++s) S.diag[s] = 1.0;
        for (std::size_t p : imm.pins) S.diag[p] = 0.0;
        for (std::size_t i = nvar; i < nz; ++i) S.diag[i] = 1.0;
        std::vector<doublereal> y(nz, 0.0);
        std::copy(xcur.begin(), xcur.end(), y.begin());
        std::copy(phi.begin(), phi.end(), y.begin() + n);
        std::copy(mu.begin(), mu.end(), y.begin() + n + nimm);
        std::copy(zeta.begin(), zeta.end(), y.begin() + n + nimm + nlatch);
        if (withcov) std::copy(Scur.begin(), Scur.end(), y.begin() + nvar);
        S.gprev = event_values(y.data(), ctx);
        S.req = &req;
        while (S.jreq < req.size() && req[S.jreq] <= tcur) ++S.jreq;
 
        // LWORK per rodas.f, N*(LJAC+LMAS+LE1+14)+20 with a full Jacobian and a diagonal mass matrix
        const integer N = static_cast<integer>(nz);
        const std::size_t lwork = nz * (nz + 1 + nz + 14) + 20 + 32, liwork = nz + 20 + 32;
        std::vector<doublereal> work(lwork, 0.0), rpar(1, 0.0);
        std::vector<integer> iwork(liwork, 0), ipar(1, 0);
        doublereal rtol = petri_coarse_tol(), atol = petri_fine_tol(), x = tcur, xend = tEnd, h = 1e-6;
        integer itol = 0, ifcn = 1, ijac = 0, mljac = N, mujac = N, idfx = 0, imas = 1, mlmas = 0, mumas = 0;
        integer iout = 1, idid = 0, lw = static_cast<integer>(lwork), liw = static_cast<integer>(liwork);
        bool failed = false;
        petri_tran_active() = &S;
        try {
            rodas_(const_cast<integer*>(&N), (U_fp)petri_tran_fcn, &ifcn, &x, y.data(), &xend, &h, &rtol, &atol,
                   &itol, (U_fp)petri_tran_jac, &ijac, &mljac, &mujac, (U_fp)petri_tran_dfx, &idfx,
                   (U_fp)petri_tran_mas, &imas, &mlmas, &mumas, (U_fp)petri_tran_solout, &iout, work.data(), &lw,
                   iwork.data(), &liw, rpar.data(), ipar.data(), &idid);
            if (idid < 0 && S.hit < 0) {
                warnings.push_back("The transient DAE failed (RODAS returned idid=" + std::to_string(idid) +
                                   "); the trajectory stops at the last accepted point.");
                failed = true;
            }
        } catch (const Error& ex) {
            warnings.push_back(std::string("The transient DAE failed (") + ex.what() +
                               "); the trajectory stops at the last accepted point.");
            failed = true;
        }
        petri_tran_active() = nullptr;
        if (!S.segT.empty()) {
            for (std::size_t a = 0; a < S.segT.size(); ++a) {
                const std::vector<double>& za = S.segZ[a];
                out.t.push_back(S.segT[a]);
                out.xt.push_back(std::vector<double>(za.begin(), za.begin() + n));
                out.phit.push_back(std::vector<double>(za.begin() + n, za.begin() + n + nimm));
                if (withcov) {
                    out.tvar.push_back(S.segT[a]);
                    Matrix<double> Sg(nc, nc, 0.0);
                    for (std::size_t i = 0; i < nc; ++i)
                        for (std::size_t j = 0; j < nc; ++j) Sg(i, j) = za[nvar + i * nc + j];
                    out.Sigmat.push_back(Sg);
                }
            }
            tcur = S.segT.back();
            const std::vector<double>& zend = S.segZ.back();
            xcur.assign(zend.begin(), zend.begin() + n);
            phi.assign(zend.begin() + n, zend.begin() + n + nimm);
            mu.assign(zend.begin() + n + nimm, zend.begin() + n + nimm + nlatch);
            zeta.assign(zend.begin() + n + nimm + nlatch, zend.begin() + nvar);
            if (withcov) Scur.assign(zend.begin() + nvar, zend.begin() + nz);
        }
        if (failed || S.hit < 0 || tcur >= tEnd - petri_fine_tol()) break;
        // the active-set move the located crossing asks for
        bool changed = false;
        const std::size_t ev = static_cast<std::size_t>(S.hit);
        if (ev < nimm) {
            if (imm.active[ev]) {
                imm.active[ev] = false;
                petri_immediate(terms, xcur, imm);
                changed = true;
            }
        } else if (ev - nimm < ncon) {
            const std::size_t c = ev - nimm;
            const std::ptrdiff_t k = index_of(active, c);
            if (k < 0) {
                active.push_back(c);
                zeta.push_back(1.0);
            } else {
                active.erase(active.begin() + k);
                zeta.erase(zeta.begin() + k);
            }
            changed = true;
        } else {
            const std::size_t s = ev - nimm - ncon;
            for (std::size_t k = 0; k < nimm; ++k) {
                const PetriMode& md = terms.modes[terms.imm_idx[k]];
                if (imm.active[k] && std::find(md.arc_slot.begin(), md.arc_slot.end(), s) != md.arc_slot.end() &&
                    imm.bind[k] != static_cast<std::ptrdiff_t>(s)) {
                    imm.bind[k] = static_cast<std::ptrdiff_t>(s);
                    petri_immediate(terms, xcur, imm);
                    changed = true;
                    break;
                }
            }
        }
        out.switches.push_back(std::make_pair(tcur, S.hit));
        if (!changed) break;
    }
    if (out.t.empty()) {
        out.t.assign(1, tStart);
        out.xt.assign(1, x0);
        out.phit.assign(1, phiss);
    }
    return out;
}
 
/**
 * The station table along a trajectory, in the conventions of the steady-state table. The rates are read at the
 * STATIONARY variance and at the trajectory's own immediate flows; the latch flow moves no marking, so it cannot
 * change a place throughput.
 */
inline void metrics_t(const PetriTerms& terms, const std::vector<std::vector<double>>& xt,
                      const std::vector<double>& s2, const std::vector<std::vector<double>>& phit, std::size_t M,
                      std::size_t K, PetriSolution& out) {
    const std::size_t nt = xt.size();
    out.QNt.assign(M, std::vector<std::vector<double>>(K, std::vector<double>(nt, 0.0)));
    out.UNt = out.QNt;
    out.TNt = out.QNt;
    std::vector<std::vector<double>> Rt(nt);
    for (std::size_t a = 0; a < nt; ++a) {
        const PetriTheta th = petri_theta(terms, xt[a], s2);
        Rt[a] = petri_rates(terms, xt[a], a < phit.size() ? phit[a] : std::vector<double>(), std::vector<double>(),
                            th);
    }
    for (std::size_t s = 0; s < terms.nm; ++s) {
        if (terms.coord_station[s] < 0 || static_cast<std::size_t>(terms.coord_station[s]) >= M) continue;
        const std::size_t i = static_cast<std::size_t>(terms.coord_station[s]), k = terms.coord_class[s];
        for (std::size_t a = 0; a < nt; ++a) {
            out.QNt[i][k][a] += xt[a][s];
            out.UNt[i][k][a] = out.QNt[i][k][a];
        }
    }
    for (std::map<std::size_t, std::vector<std::size_t>>::const_iterator it = terms.consumers.begin();
         it != terms.consumers.end(); ++it) {
        const std::size_t i = it->first / K, k = it->first % K;
        const std::vector<double>& w = terms.consumer_w.at(it->first);
        for (std::size_t a = 0; a < nt; ++a)
            for (std::size_t q = 0; q < it->second.size(); ++q) out.TNt[i][k][a] += w[q] * Rt[a][it->second[q]];
    }
    for (std::map<std::size_t, std::vector<std::size_t>>::const_iterator it = terms.producers.begin();
         it != terms.producers.end(); ++it) {
        const std::size_t i = it->first / K, k = it->first % K;
        for (std::size_t a = 0; a < nt; ++a)
            for (std::size_t e : it->second) out.TNt[i][k][a] += Rt[a][e];
    }
    // Along a trajectory the empty-level gate is read off the level itself: at an empty place the outflow is cut
    // to the inflow, which is the value the binding multiplier takes at that instant.
    for (std::size_t q = 0; q < terms.nf; ++q) {
        if (terms.fluid_station[q] < 0 || static_cast<std::size_t>(terms.fluid_station[q]) >= M || K == 0) continue;
        const std::size_t i = static_cast<std::size_t>(terms.fluid_station[q]), slot = terms.fluid_slot[q];
        double levmax = 0.0;
        for (std::size_t a = 0; a < nt; ++a) levmax = std::max(levmax, std::fabs(xt[a][slot]));
        const double tol_empty = petri_coarse_tol() * std::max(1.0, levmax);
        for (std::size_t a = 0; a < nt; ++a) {
            const double lev = xt[a][slot];
            out.QNt[i][0][a] = lev;
            out.UNt[i][0][a] = std::isfinite(terms.fluid_cap[q]) ? lev / terms.fluid_cap[q] : lev;
            double outf = 0.0, inf = 0.0;
            for (std::size_t e = 0; e < terms.nev; ++e) {
                const double d = terms.D(slot, e);
                if (d < 0) outf -= d * Rt[a][e];
                else if (d > 0) inf += d * Rt[a][e];
            }
            if (lev <= tol_empty) outf = std::min(outf, inf);
            out.TNt[i][0][a] = outf;
        }
    }
}
 
}  // namespace solver_detail
 
/**
 * Fluid analysis of a stochastic Petri net.
 *
 * @param sn  a model whose nodes are Places, Transitions, Sources and Sinks only
 * @param opt the state-size cap, the Newton tolerance and its iteration cap
 */
template <class T>
inline PetriSolution solver_fluid_petri(const qn::NetworkStruct<T>& sn,
                                        const PetriOptions& opt = PetriOptions()) {
    const std::size_t M = sn.nstations, K = sn.nclasses;
 
    const PetriVerdict v = petri_applicable(sn);
    if (!v.ok)
        throw UnsupportedError(ErrorKind{"Model"}, "solver_fluid_petri: the fluid Petri route cannot solve this "
                               "model: " + v.reason + ".");
 
    const PetriTerms terms = petri_build_terms(sn);
    const std::size_t n = terms.nstate, npair = terms.npair;
    if (n > opt.dae_maxstate)
        throw UnsupportedError(ErrorKind{"StateSpace"},
            "solver_fluid_petri: the fluid Petri route solves a " + std::to_string(n) +
            "-unknown algebraic system with a finite-difference Jacobian, above the limit of " +
            std::to_string(opt.dae_maxstate) +
            " set by options.config.dae_maxstate. Raise that limit, or use SolverSSA for a net "
            "of this size");
 
    const PetriConservation cons = petri_conservation(terms);
    if (cons.leak > 1e-7)
        throw NumericError(ErrorKind{"Degenerate"}, "solver_fluid_petri: the conserved directions and the jump matrix "
                           "disagree: the largest leak per unit rate is " +
                           std::to_string(cons.leak) + ", where it must be zero");
    const PetriConstraints con = petri_constraints(sn, terms);
    const std::size_t ncon = con.b.size();
 
    // ---- seed ---------------------------------------------------------------
    // Newton needs a point in the basin, not an answer. The immediate modes get a
    // large FINITE rate here and only here, scaled to the model's own timescale
    // rather than taken from GlobalConstants.Immediate: 1e8 against a rate of
    // order one is a stiffness the seed does not need, and the answer does not
    // depend on the seed's accuracy.
    double rmax = 0.0;
    for (std::size_t j : terms.timed_idx) {
        double s = 0.0;
        for (double d : terms.modes[j].d1) s += d;
        rmax = std::max(rmax, s);
    }
    for (std::size_t e = 0; e < terms.nev; ++e)
        if (terms.ev_kind[e] == 3 || terms.ev_kind[e] == 6) rmax = std::max(rmax, terms.rate_base[e]);
    if (!(rmax > 0)) rmax = 1.0;
    const double lam = std::min(1e8, 1e4 * rmax);
    double mass = 0.0;
    for (std::size_t s = 0; s < terms.nm; ++s) mass += terms.x0[s];
    // A fluid level has to reach its bound or its balance within the horizon too, so its initial level and any
    // finite capacity count as mass of the same units.
    for (std::size_t q = 0; q < terms.nf; ++q) {
        mass += terms.fluid_init[q];
        if (std::isfinite(terms.fluid_cap[q])) mass += terms.fluid_cap[q];
    }
    double Thor = 50.0 * (mass + 1.0) / rmax;
 
    const auto drift_fn = [&terms, lam](const double&, const std::vector<double>& y) {
        return solver_detail::seed_drift(terms, y, lam);
    };
    std::vector<double> tseed;
    std::vector<std::vector<double>> xseed;
    for (int attempt = 0; attempt < 6; ++attempt) {
        OdeOptions<double> oo;
        oo.rtol = 1e-6;
        oo.atol = 1e-10;
        oo.store_trajectory = true;
        oo.max_steps = 200000;
        const OdeSolution<double> sol =
            ode_rosenbrock4(drift_fn, 0.0, Thor, terms.x0, oo);
        tseed = sol.t;
        xseed = sol.y;
        const std::vector<double> d = solver_detail::seed_drift(terms, xseed.back(), lam);
        if (solver_detail::inf_norm(d) <= 1e-6 * std::max(1.0, rmax * (mass + 1.0))) break;
        Thor *= 4.0;
    }
    std::vector<double> x = xseed.back();
 
    PetriImmediate imm;
    petri_immediate(terms, x, imm);
    const std::size_t nimm = imm.n;
 
    // THE VARIANCE IS SEEDED POSITIVE: sigma2 = 0 is where min() has no
    // derivative, and a saturated net's first-order fixed point sits there.
    std::vector<double> s2(npair, 0.0);
    std::vector<bool> ondiag(npair, false);
    for (std::size_t p = 0; p < npair; ++p) {
        ondiag[p] = (terms.cov_pairs[p].first == terms.cov_pairs[p].second);
        if (ondiag[p]) s2[p] = std::max(petri_fine_tol(), x[terms.cov_pairs[p].first]);
    }
 
    // ---- steady state -------------------------------------------------------
    // A FLUID BOUND IS READ OFF THE SEED. A level whose flows do not depend on it has a constant column of drift,
    // so a Newton step cannot move it onto a bound; the seed has already run it there, and a bound it presses
    // against is binding.
    std::vector<std::size_t> active = solver_detail::fluid_seed_active(terms, con, x);
    std::vector<double> phi(nimm, 0.0);
    const std::size_t nlatch = terms.latch_mode.size();
    std::vector<double> mu(nlatch, 0.0), zeta;
    std::size_t iters = 0;
    const std::size_t aset_max = std::max<std::size_t>(6, 2 * (ncon + nimm) + 2);
    bool converged = false;
    double resnorm = std::numeric_limits<double>::infinity();
    bool have_best = false;
    std::vector<double> bx, bs2, bphi, bmu, bzeta;
    std::vector<std::size_t> bactive;
    PetriImmediate bimm;
    double bres = std::numeric_limits<double>::infinity();
    bool bconv = false;
 
    for (std::size_t sweep = 0; sweep < aset_max; ++sweep) {
        const solver_detail::Ctx ctx =
            solver_detail::context(terms, cons, con, imm, active);
        std::vector<double> u0;
        u0.insert(u0.end(), x.begin(), x.end());
        u0.insert(u0.end(), s2.begin(), s2.end());
        u0.insert(u0.end(), phi.begin(), phi.end());
        u0.insert(u0.end(), mu.begin(), mu.end());
        u0.insert(u0.end(), active.size(), 1.0);
        std::vector<double> lb(u0.size(), -std::numeric_limits<double>::infinity());
        for (std::size_t k = 0; k < nimm; ++k) lb[n + npair + k] = 0.0;
        for (std::size_t k = 0; k < active.size(); ++k)
            lb[n + npair + nimm + nlatch + k] = 0.0;
        for (std::size_t p = 0; p < npair; ++p)
            if (ondiag[p]) lb[n + p] = 0.0;
 
        const solver_detail::NewtonResult nr =
            solver_detail::newton(ctx, u0, opt.tol, opt.newton_max, lb);
        iters += nr.iterations;
        resnorm = nr.resnorm;
        converged = nr.converged;
        const solver_detail::Unpacked up =
            solver_detail::unpack(nr.u, terms, imm, active.size());
        x = up.x;
        s2 = up.s2;
        phi = up.phi;
        mu = up.mu;
        zeta = up.zeta;
 
        if (!have_best || resnorm < bres) {
            bx = x; bs2 = s2; bphi = phi; bmu = mu; bzeta = zeta;
            bactive = active; bimm = imm; bres = resnorm; bconv = converged;
            have_best = true;
        }
 
        // The active-set moves, in the order that a failure of each invalidates
        // the next: a negative flow means the mode does not fire at all, a
        // negative marking means the wrong coordinate was pinned, and only then
        // is it worth asking which capacity rows bind.
        bool moved = false;
        if (nimm > 0) {
            double scale = 1.0;
            for (double p : phi) scale = std::max(scale, std::fabs(p));
            const double thr = -std::max(opt.tol, 1e-10) * scale;
            for (std::size_t k = 0; k < nimm && !moved; ++k)
                if (phi[k] < thr) {
                    imm.active[k] = false;
                    petri_immediate(terms, x, imm);
                    moved = true;
                }
        }
        if (!moved && nimm > 0) {
            const double thr = -std::max(opt.tol, 1e-10);
            for (std::size_t s = 0; s < terms.nm && !moved; ++s) {
                if (x[s] >= thr) continue;
                for (std::size_t k = 0; k < nimm; ++k) {
                    const PetriMode& md = terms.modes[terms.imm_idx[k]];
                    if (imm.active[k] &&
                        std::find(md.arc_slot.begin(), md.arc_slot.end(), s) !=
                            md.arc_slot.end() &&
                        imm.bind[k] != static_cast<std::ptrdiff_t>(s)) {
                        imm.bind[k] = static_cast<std::ptrdiff_t>(s);
                        petri_immediate(terms, x, imm);
                        moved = true;
                        break;
                    }
                }
            }
        }
        if (!moved && ncon > 0) {
            std::vector<std::size_t> over;
            for (std::size_t c = 0; c < ncon; ++c) {
                double val = 0.0;
                for (std::size_t s = 0; s < n; ++s) val += con.A(c, s) * x[s];
                if (val > con.b[c] + std::max(1e-9, opt.tol) &&
                    std::find(active.begin(), active.end(), c) == active.end())
                    over.push_back(c);
            }
            if (!over.empty()) {
                active.insert(active.end(), over.begin(), over.end());
                moved = true;
            } else {
                std::vector<std::size_t> keep;
                for (std::size_t i = 0; i < active.size(); ++i)
                    if (!(i < zeta.size() && zeta[i] > 1 + std::max(1e-9, opt.tol)))
                        keep.push_back(active[i]);
                if (keep.size() != active.size()) {
                    active = keep;
                    moved = true;
                }
            }
        }
        if (!moved) break;
    }
 
    PetriSolution out;
    if (have_best && !converged && bconv) {
        x = bx; s2 = bs2; phi = bphi; mu = bmu; zeta = bzeta;
        active = bactive; imm = bimm; converged = true; resnorm = bres;
    }
    if (!converged)
        out.warnings.push_back(
            "The simultaneous closure solve stopped at residual " + std::to_string(resnorm) +
            " after " + std::to_string(iters) + " Newton steps without reaching " +
            std::to_string(opt.tol) + ". The reported point is the last iterate.");
    // A NEGATIVE MARKING IS NOT ROUNDING: the fixed point wanted mass a place
    // cannot supply and no immediate mode could be rebound to pin it, so the
    // answer is outside the model's own state space and is reported as such.
    for (std::size_t s = 0; s < terms.nm; ++s)
        if (x[s] < -std::max(opt.tol, 1e-10)) {
            out.warnings.push_back(
                "The fixed point holds " + std::to_string(x[s]) + " tokens at " +
                terms.names_node[terms.coord_node[s]] +
                ", which is negative: no immediate transition could be rebound to pin that "
                "place at zero. Use SolverCTMC, SolverSSA or SolverLDES for this net.");
            break;
        }
 
    const solver_detail::Ctx ctx = solver_detail::context(terms, cons, con, imm, active);
    std::vector<double> ufin;
    ufin.insert(ufin.end(), x.begin(), x.end());
    ufin.insert(ufin.end(), s2.begin(), s2.end());
    ufin.insert(ufin.end(), phi.begin(), phi.end());
    ufin.insert(ufin.end(), mu.begin(), mu.end());
    ufin.insert(ufin.end(), zeta.begin(), zeta.end());
    const solver_detail::Residual fin = solver_detail::residual(ufin, ctx, false);
 
    // ---- the metrics --------------------------------------------------------
    // A PLACE IS AN INF STATION: its queue length is its mean token count and its
    // utilization is the same number. Its throughput is the rate at which TOKENS
    // leave it, so a consuming mode contributes its firing rate times the arc
    // multiplicity -- SolverCTMC's convention, and the one Little's law needs.
    out.QN = Matrix<double>(M, K, 0.0);
    out.UN = Matrix<double>(M, K, 0.0);
    out.RN = Matrix<double>(M, K, 0.0);
    out.TN = Matrix<double>(M, K, 0.0);
    for (std::size_t s = 0; s < terms.nm; ++s) {
        if (terms.coord_station[s] < 0) continue;
        const std::size_t i = static_cast<std::size_t>(terms.coord_station[s]);
        const std::size_t k = terms.coord_class[s];
        out.QN(i, k) += x[s];
        out.UN(i, k) = out.QN(i, k);
    }
    for (std::map<std::size_t, std::vector<std::size_t>>::const_iterator it =
             terms.consumers.begin();
         it != terms.consumers.end(); ++it) {
        const std::size_t i = it->first / K, k = it->first % K;
        const std::vector<double>& w = terms.consumer_w.at(it->first);
        double val = 0.0;
        for (std::size_t q = 0; q < it->second.size(); ++q) val += w[q] * fin.r[it->second[q]];
        out.TN(i, k) += val;
    }
    for (std::map<std::size_t, std::vector<std::size_t>>::const_iterator it =
             terms.producers.begin();
         it != terms.producers.end(); ++it) {
        const std::size_t i = it->first / K, k = it->first % K;
        double val = 0.0;
        for (std::size_t e : it->second) val += fin.r[e];
        out.TN(i, k) += val;
    }
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t k = 0; k < K; ++k)
            if (out.TN(i, k) > 1e-14) out.RN(i, k) = out.QN(i, k) / out.TN(i, k);
    // A FLUID PLACE reports its level in the first class column: QN is the mean level (the reflected-diffusion
    // one under fluid_boundary='diffusion'), UN the fraction of a finite capacity it fills (the level itself when
    // unbounded), and TN the rate at which fluid LEAVES it, flows out plus negative impulses, gated at the empty
    // level, so that RN = QN/TN is the mean time a unit of fluid spends in the place by Little's law.
    std::vector<double> gP, gN;
    solver_detail::gains(con, active, zeta, n, gP, gN);
    const PetriReport::Fluid fl = solver_detail::fluid_levels(terms, con, imm, active, x, s2, phi, fin.r,
                                                              fin.Sigma, opt.fluid_boundary);
    for (std::size_t q = 0; q < terms.nf; ++q) {
        if (terms.fluid_station[q] < 0 || K == 0) continue;
        const std::size_t i = static_cast<std::size_t>(terms.fluid_station[q]);
        const std::size_t slot = terms.fluid_slot[q];
        // the reported level is clamped to [0, B]: a level pinned at a bound comes out of the solve a rounding
        // error past it (-9e-20 under 'reflect'), which the metrics must not show as a negative QLen
        const double lev = std::min(std::max(fl.level[q], 0.0), terms.fluid_cap[q]);
        out.QN(i, 0) = lev;
        out.UN(i, 0) = std::isfinite(terms.fluid_cap[q]) ? lev / terms.fluid_cap[q] : lev;
        double outflow = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e) outflow += std::min(terms.D(slot, e), 0.0) * fin.r[e];
        out.TN(i, 0) = -gN[slot] * outflow;
        out.RN(i, 0) = out.TN(i, 0) > 0.0 ? out.QN(i, 0) / out.TN(i, 0) : 0.0;
    }
 
    // ---- the trajectory -----------------------------------------------------
    // The immediate flow ALONG the reported path is not the steady-state one: on the seed path it is the
    // large-finite-rate approximation the seed itself integrated, so the throughput table and the trajectory it
    // is read off describe the same run. A finite horizon replaces the seed path by the transient DAE.
    out.t = tseed;
    out.xvec_t = xseed;
    std::vector<std::vector<double>> phit = solver_detail::seed_flows(terms, xseed, lam);
    if (std::isfinite(opt.timespan_end)) {
        solver_detail::PetriTransient tr =
            solver_detail::petri_transient(terms, cons, con, imm, active, opt, phi, zeta, out.warnings);
        out.t = tr.t;
        out.xvec_t = tr.xt;
        phit = tr.phit;
        out.Sigmat = tr.Sigmat;
        out.tvar = tr.tvar;
        out.switches = tr.switches;
    }
    solver_detail::metrics_t(terms, out.xvec_t, s2, phit, M, K, out);
    out.x = x;
    out.iters = iters;
    out.resnorm = resnorm;
    out.converged = converged;
    out.Sigma = fin.Sigma;
    out.QVar = Matrix<double>(M, K, 0.0);
    out.QStd = Matrix<double>(M, K, 0.0);
    for (std::size_t s = 0; s < terms.nm; ++s) {
        if (terms.coord_station[s] < 0) continue;
        const std::size_t i = static_cast<std::size_t>(terms.coord_station[s]);
        out.QVar(i, terms.coord_class[s]) = std::max(0.0, fin.Sigma(s, s));
    }
    for (std::size_t q = 0; q < terms.nf; ++q)
        if (terms.fluid_station[q] >= 0 && K > 0)
            out.QVar(static_cast<std::size_t>(terms.fluid_station[q]), 0) = fl.level_var[q];
    for (std::size_t i = 0; i < M; ++i)
        for (std::size_t k = 0; k < K; ++k) out.QStd(i, k) = std::sqrt(out.QVar(i, k));
 
    // ---- the report ---------------------------------------------------------
    out.petri.marking = Matrix<double>(terms.I, terms.K, 0.0);
    out.petri.marking_var = Matrix<double>(terms.I, terms.K, 0.0);
    for (std::size_t s = 0; s < terms.nm; ++s) {
        out.petri.marking(terms.coord_node[s], terms.coord_class[s]) = x[s];
        out.petri.marking_var(terms.coord_node[s], terms.coord_class[s]) =
            std::max(0.0, fin.Sigma(s, s));
    }
    for (const PetriMode& md : terms.modes) out.petri.mode_label.push_back(md.label);
    out.petri.mode_flow.assign(terms.modes.size(), 0.0);
    for (std::size_t j = 0; j < terms.modes.size(); ++j) {
        double val = 0.0;
        for (std::size_t e = 0; e < terms.nev; ++e)
            if (terms.ev_mode[e] == static_cast<int>(j) &&
                (terms.ev_kind[e] == 1 || terms.ev_kind[e] == 4))
                val += fin.r[e];
        out.petri.mode_flow[j] = val;
    }
    out.petri.immediate_flow = phi;
    out.petri.invariant_label = cons.label;
    out.petri.invariant_value = cons.N;
    out.petri.invariant_error.assign(cons.C.rows(), 0.0);
    for (std::size_t c = 0; c < cons.C.rows(); ++c) {
        double val = 0.0;
        for (std::size_t s = 0; s < terms.nstate; ++s) val += cons.C(c, s) * x[s];
        out.petri.invariant_error[c] = val - cons.N[c];
    }
    out.petri.capacity_label = con.label;
    out.petri.capacity_active = active;
    out.petri.capacity_fraction = zeta;
    out.petri.pinned = imm.pins;
    out.petri.Sigma = fin.Sigma;
    out.petri.fluid = fl;
    return out;
}
 
}  // namespace petri
}  // namespace fluid
}  // namespace line
 
#endif  // LINE_SOLVERS_FLUID_PETRI_H