// ============================================================================
//  SyFox — si_substrate.hpp
//  ---------------------------------------------------------------------------
//  THE CORE. An INDEPENDENT REIMPLEMENTATION of the Synthetic-Intelligence
//  substrate concept (Mr-DS-ML-85/Synthetic-Intelligence) — not a port of its
//  physics. HONEST LINEAGE (v3.9 audit; what is shared vs what differs):
//
//    SHARED with SI: the substrate ARCHITECTURE — concept nodes, weighted
//    lanes, energy injection, dissipative settle, a salience integrator of
//    the physics.hpp cavity lineage (s = tanh(s·decay + gain·|motion|)),
//    Hebbian learning with decay and lane caps, resonance-sweep readout,
//    honest silence. The salience_gating / miller_window selection modes are
//    direct ports of the SI selection layer (TSDA live-cap "7 ± 2").
//
//    DIFFERENT from SI (by design, for scale and bit-determinism): the state
//    variable is a single scalar energy per node — NOT SI's mass-spring
//    oscillator (x, v with x_E/v_E/x_I/v_I coupled fields, physics.hpp:53);
//    edges are undirected weight lanes — NOT SI's tunnels with stiffness /
//    polarity / resonance_omega / Q-gating (physics.hpp:117-131); time is
//    pass-counted settle iterations with an eps early stop — NOT SI's
//    dt=0.02 fixed-tick integration; decay 0.82 and diffusion 0.45 are
//    SyFox-native constants with no SI source (SI damps velocity at
//    friction=0.006 and derives coupling per relation from fan-out). SI's
//    Layer-2 constructive substrate (typed graph, nine reasoning primitives,
//    chaining, proof traces) now lives beside this field in core/reason.hpp —
//    ported faithfully THERE, v3.9.
//
//  Nothing here is a transformer, neural network, or pattern-matching
//  classifier. Decisions emerge from field dynamics only:
//
//    1. INJECTION   — state tokens deposit energy on concept nodes,
//                     weighted by acoustic mass (frequent concepts are
//                     "heavier" and carry less per-token energy).
//    2. SETTLE      — energy diffuses along Hebbian lanes with decay,
//                     iterated K_settle passes; an energy-gated source cap
//                     keeps the working set bounded. The SI salience
//                     integrator (physics.hpp cavity lineage:
//                     s = tanh(s·decay + gain·|motion|)) runs every pass and
//                     is observable on every node. Two flags enable the
//                     faithful SI selection modes: salience_gating ranks
//                     sources by salience instead of raw energy; miller_window
//                     samples the live cap from [source_cap-4, source_cap]
//                     per decision (TSDA "7 ± 2"; SI draws it from a seeded
//                     mt19937, SyFox derives it from the state hash so the
//                     same state still settles bit-for-bit). The field stops
//                     when the total energy delta drops below eps.
//    3. READOUT     — a resonance sweep reads the settled field at probe
//                     nodes (options / anchors / statements). Readout is a
//                     measurement of the field, not a classifier.
//    4. HEBBIAN     — labelled experience binds lanes between co-active
//                     state and option concepts (fire together, wire
//                     together), with weight decay and per-node lane caps.
//    5. HONEST      — if the settled field is below the silence floor the
//         SILENCE    substrate reports "no collapse" instead of guessing.
//
//  Zero dependencies. Deterministic. C++17.
// ============================================================================
#pragma once
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <tuple>
#if defined(_OPENMP)
#include <omp.h>
#endif
#include <fstream>
#include <map>
#include <string>
#include <unordered_map>
#include <unordered_set>
#include <vector>

namespace si {

using NodeId = std::uint32_t;

// ---------------------------------------------------------------------------
// v3.2 SEMANTIC FIELD — the substrate's second field, deterministic, no ML.
// Every concept node carries a 64-dim semantic vector built in two
// deterministic steps (see build_semantics):
//   1. LEXICAL LAYER   — signed character-trigram hashing of the concept
//      string into 64 buckets ("card_arrival" and "card_delivery" share
//      trigrams, so their base vectors align; unrelated words stay
//      near-orthogonal). A fixed function of the string: no training.
//   2. FABRIC GROUNDING — two smoothing passes over the Hebbian lanes
//      (v_i <- normalize(v_i + beta * lane-weighted mean of v_j)). The
//      fabric's own co-occurrence structure shapes the vectors — a
//      deterministic diffusion of meaning along the lanes, not a fit.
// From the vectors, per node:
//   omega_semantic  — scalar projection onto one fixed direction (Stage 1
//                     "concept frequency encoding"); close frequencies
//                     resonate, like coupled oscillators detuned.
//   resonance edges — top-k cosine neighbours (cos >= sem_theta), weight
//                     = cos * frequency-match (Stage 4 "semantic field":
//                     during settle, sources leak a small fraction of
//                     their retained energy to semantically similar nodes
//                     — energy is CONSERVED, the flow is one more
//                     diffusion channel beside the Hebbian lanes).
//   semantic type   — cos >= 0.82 reads as synonym-grade kin, >= 0.66 as
//                     related (Option 2 typed lanes, derived from the
//                     vectors, not from a hand-written table).
// CONTEXT-SENSITIVE LANES (Stage 2): a lane may carry a context signature
// (required/forbidden context words, learned from the co-occurrence
// statistics of the lessons that laid it). At settle, the lane's flow is
// damped when the current decision's injected tokens do not match the
// signature: all required present -> full flow, missing -> ctx_missing,
// any forbidden present -> ctx_forbidden. The same "card" can then route
// toward card_arrival when arrive/receive are in the state and toward
// card_delivery_estimate when estimate/how-long is — contextual
// disambiguation as dampened coupling, not pattern matching.
// Persistence: a v4 tail in substrate.bin. Models saved before v3.2 have no
// tail: has_semantics() is false and every semantic code path is inert, so
// old fabrics replay bit-for-bit. Runtime kill switch: set_semantics(false).
// ---------------------------------------------------------------------------
inline constexpr std::size_t SEM_DIMS = 64;

// v3.4 multi-hop readout walk: max nodes expanded per hop level (smallest
// NodeIds first, deterministic). Bounds the walk at O(levels x 64).
inline constexpr std::size_t kHopWidthCap = 64;

// v3.2 context signature attached to a lane (read-side of Stage 2).
struct LaneCtx {
    std::vector<NodeId> required;    // all of these should be in the state
    std::vector<NodeId> forbidden;   // any of these suppresses the lane
};

// ---------------------------------------------------------------------------
// v3 Milestone 3 — the lane evidence ledger (audit surface, NOT physics).
// A lane's WEIGHT drives the field (settle/readout) exactly as before; the
// ledger records HOW a lane earned its weight, so every decision can show
// its evidence and every derived relation can be audited:
//   support_events  — lessons that strengthened this lane (Hebbian binds)
//   counter_events  — anti-Hebbian weakenings (disconfirming evidence)
//   first_seq/last_seq — the teach-event window that touched the lane
//                        (provenance; full per-lesson ids live in the
//                        model's conflicts/ledger files)
//   context         — the engine context tag at last touch (model/domain)
// The ledger never feeds back into settle/readout: it is bookkeeping over
// the physics, the way a lab notebook accompanies an experiment.
// ---------------------------------------------------------------------------
struct LaneEvidence {
    std::uint32_t support_events = 0;
    std::uint32_t counter_events = 0;
    std::uint64_t first_seq = 0;
    std::uint64_t last_seq  = 0;
    std::string   context;
};

// ---------------------------------------------------------------------------
struct SubstrateConfig {
    float decay          = 0.82f;  // per-pass energy retention (damping)
    float diffusion      = 0.45f;  // fraction of energy that flows along lanes
    int   k_settle       = 8;      // max settle passes
    float eps            = 1e-4f;  // phase-stability threshold (early stop)
    float source_cap     = 24.0f;  // max propagation sources per pass, top-K BY ENERGY
                                   // (SI-inspired bounded working set; fixed by default)
    // -- SI salience lineage (physics.hpp cavity integrator + TSDA SalienceMap) --
    float salience_gain   = 0.05f;  // SI physics.hpp salience_gain
    float salience_decay  = 0.995f; // SI physics.hpp salience_decay
    bool  salience_gating = false;  // true: rank sources by salience (motion history), not energy
    bool  miller_window   = false;  // true: live cap sampled from [source_cap-4, source_cap]
                                    // per decision (TSDA live_cap; deterministic per state)
    float silence_floor  = 0.05f;  // below this total energy => honest silence
    float inject_energy  = 1.0f;   // energy per fresh token
    float learn_eta      = 0.10f;  // Hebbian learning rate
    float lane_decay     = 0.995f; // per-lesson lane weight decay
    std::size_t lane_cap = 256;    // max lanes per node (weakest evicted)
    float hop_coupling   = 0.35f;  // 1-hop lane contribution at readout
    // -- v3.2 semantic field constants ---------------------------------------
    float sem_coupling  = 0.12f;  // resonance outflow fraction of retained energy
    float sem_theta     = 0.50f;  // min cosine for a resonance edge
    float sem_hop       = 0.10f;  // semantic-neighbour contribution at readout
    float sem_beta      = 0.50f;  // fabric-grounding strength (vector smoothing)
    int   sem_neighbors = 6;      // resonance edges per node
    float ctx_missing   = 0.60f;  // lane flow factor when required context is absent
    float ctx_forbidden = 0.20f;  // lane flow factor when forbidden context is present
    // -- M3: COMPLEMENTARY LEARNING SYSTEMS (hippocampal / cortical split) ----
    // One store must currently hold BOTH the taxonomy and every episode, under a
    // single ~138-lesson half-life (lane_decay 0.995). Measured consequence: SST-2
    // one-pass 0.896 vs 4-chunk 0.814 — chunking LOST, because the re-intern
    // penalty roughly cancels the decay benefit when one store serves both roles.
    //
    // McClelland's Complementary Learning Systems says the brain runs two: the
    // hippocampus encodes individual episodes FAST and sparsely, the neocortex
    // extracts slowly-learned regular structure, and replay moves episodes into
    // the slow system where they GENERALISE (only regularities survive).
    //
    // Physical form here: two lane populations.
    //   hippocampal — the existing out_ store, decays per lesson as before.
    //   cortical   — an additional, near-non-decaying conductance population that
    //                only receives a lane once that lane's eligibility trace
    //                crosses cortical_gate. That gate IS the consolidation step:
    //                coincidence across lessons, not a single binding, is what
    //                makes structure cortical. Both feed the same CSR, so
    //                cortical structure propagates like any other conductance.
    // -- M5: NTDA (action-potential thresholding, edge death, modulators) ------
    // Three mechanisms from ntda.md, all expressed as physics on lanes rather
    // than as heuristics in the readout:
    //
    // (a) THRESHOLD GATING. A lane accumulates VOLTAGE on each binding
    //     opportunity and only actually binds when that voltage reaches theta;
    //     then it RESETS. So a lane cannot fire twice without being
    //     re-stimulated. This replaces "one bind per lesson" with a real
    //     refractory period, and it is what lets rare conjunctions accumulate
    //     instead of being averaged into noise on their first appearance.
    // (b) EDGE DEATH. Conductivity decays per lesson and lanes below
    //     edge_death are DELETED, so the ontology physically reshapes around
    //     what the fabric still uses. SyFox had no concept of an edge dying.
    // (c) MODULATORS. theta is lowered by novelty (dopamine: risk a wild
    //     analogy when the input is unfamiliar) and raised by contradiction
    //     (cortisol: hyper-focus on the max-conductivity path). In a substrate
    //     whose accuracy collapses to deferral, an explicit exploration term
    //     keyed to novelty is the missing piece.
    bool  ntda_on        = false;   // M5 master gate (default off == pre-M5)
    float ntda_theta     = 1.0f;    // voltage needed to fire a lane
    float ntda_dopamine  = 0.30f;   // max theta reduction at novelty = 1
    float ntda_cortisol  = 0.50f;   // max theta increase at contradiction = 1
    float edge_decay     = 0.9995f; // per-lesson conductivity decay
    float edge_death     = 0.002f;  // delete a lane below this weight
    bool  cls_split      = false;  // M3 master gate (default off == pre-M3)
    float cortical_decay = 0.9999f;// retained for disclosure; cortical lanes are
                                   // not touched by enforce_lane_decay at all
    float cortical_gate  = 3.0f;   // eligibility (lesson co-occurrences) needed
    float cortical_gain  = 0.5f;   // readout weight of a cortical vs hip lane
};

inline float fnv1a_hash(const std::string& s) {  // deterministic token hash
    std::uint64_t h = 1469598103934665603ull;
    for (unsigned char c : s) { h ^= c; h *= 1099511628211ull; }
    return static_cast<float>(h & 0xFFFFFFull);
}

// ---------------------------------------------------------------------------
// ConceptField — nodes, lanes, and the field dynamics over them.
// ---------------------------------------------------------------------------
class Substrate {
public:
    explicit Substrate(SubstrateConfig cfg = {}) : cfg_(cfg) {}

    // -- vocabulary ---------------------------------------------------------
    NodeId intern(const std::string& concept) {
        auto it = index_.find(concept);
        if (it != index_.end()) { nodes_[it->second].mass += 1.0f; return it->second; }
        // v3.2: the semantic field is a property of the WHOLE fabric — a new
        // node changes every top-k neighbourhood, so the stored field is now
        // stale. Invalidate rather than index out of bounds: the next
        // save_model() rebuilds it over the grown fabric (and in-memory
        // engines run lane-physics only until then — honest degradation).
        if (!semvecs_.empty()) {
            semvecs_.clear(); sem_edges_.clear(); sem_off_.clear();
            sem_edge_count_ = 0;
        }
        NodeId id = static_cast<NodeId>(nodes_.size());
        index_.emplace(concept, id);
        // occurrence mass: a LINEAR count of intern() calls (not log-compressed).
        // The sub-linear behavior lives at the USE sites: injection deposits
        // energy/sqrt(mass) and Hebbian bind strength scales by 1/sqrt(mass),
        // so heavy concepts move less.
        nodes_.push_back(Node{concept, 1.0f, 0.0f, 0.0f});
        csr_dirty_ = true;                       // node count changed (M5 CSR mirror)
        return id;
    }

    bool has(const std::string& c) const { return index_.count(c) != 0; }
    NodeId find(const std::string& c) const { return index_.at(c); }
    std::size_t node_count() const { return nodes_.size(); }
    std::size_t lane_count() const {
        std::size_t n = 0;
        for (const auto& kv : out_) n += kv.second.size();
        return n;
    }
    std::size_t readout_neighbours(NodeId a) const {
        auto it = out_.find(a);
        return it == out_.end() ? 0 : it->second.size();
    }

    // -- v3.6.0 DOCUMENT-FREQUENCY TABLE (CMD §2.2 rarity grounding; MIMS §6
    // fan-out derivation) -----------------------------------------------
    // The engine calls begin_df_lesson() once per lesson and note_df() once
    // per (lesson, node), so df(n) = number of lessons containing n. The
    // rarity factor log1p(L/(1+df))/log1p(L) ∈ (0,1] then scales the
    // STATE-side bind weights in hebbian_lesson: a token that fires in every
    // lesson distinguishes nothing and binds weakly; a rare token binds
    // strongly. This is CMD's "ubiquitous participants are suppressed by
    // construction, not by a rule" — no stopword list, one counted formula.
    // Active only when the table is non-empty (fresh learns note df; fabrics
    // loaded from a pre-v3.6 file have no 'IDF5' tail and every factor is
    // 1.0, so the replay contract holds bit for bit).
    void begin_df_lesson() { ++df_lessons_; }
    void note_df(NodeId id) { ++df_[id]; }
    std::uint32_t df_of(NodeId id) const {
        auto it = df_.find(id);
        return it == df_.end() ? 0u : it->second;
    }
    std::uint64_t df_lesson_count() const { return df_lessons_; }
    std::size_t df_table_size() const { return df_.size(); }
    // M9-LEARNED: the document-frequency table, exposed so the role tagger can
    // learn what a FUNCTION WORD is from the fabric instead of consulting a
    // hardcoded list. A token seen in a large fraction of lessons is a function
    // word; that is a property of the DATA, not of the source file.
    std::uint32_t node_df(NodeId id) const {
        const auto it = df_.find(id);
        return it == df_.end() ? 0u : it->second;
    }
    // Top-N most frequent concepts, as (label, df) sorted by df desc then label.
    // Deterministic: ties broken by label, never by hash order.
    std::vector<std::pair<std::string, std::uint32_t>> df_ranked(
        std::size_t top_n) const {
        std::vector<std::pair<std::string, std::uint32_t>> out;
        out.reserve(df_.size());
        for (const auto& kv : df_)
            if (kv.first < nodes_.size())
                out.emplace_back(nodes_[kv.first].concept, kv.second);
        std::sort(out.begin(), out.end(),
                  [](const std::pair<std::string, std::uint32_t>& x,
                     const std::pair<std::string, std::uint32_t>& y) {
                      if (x.second != y.second) return x.second > y.second;
                      return x.first < y.first;
                  });
        if (out.size() > top_n) out.resize(top_n);
        return out;
    }
    bool idf_active() const { return !df_.empty() && df_lessons_ > 0; }
    float idf_factor(NodeId id) const {
        if (!idf_active()) return 1.0f;
        const double L = static_cast<double>(df_lessons_);
        const double d = static_cast<double>(df_of(id));
        return static_cast<float>(std::log1p(L / (1.0 + d)) / std::log1p(L));
    }

    // -- Hebbian lane binding (fire together -> wire together) ---------------
    void bind(NodeId a, NodeId b, float w) {
        if (a == b || w <= 0.0f) return;
        add_lane(a, b, w);
        add_lane(b, a, w);
    }

    // M3: eligibility + cortical consolidation. A lane reaches the cortical
    // population only after cortical_gate DISTINCT lesson co-occurrences, which
    // is the physical form of "replay consolidates a regularity, not an
    // episode". Saturating, like add_lane, so it is order-insensitive once a
    // lane is established (the CLAMP is why cortical lanes are safe to merge
    // across passes where hippocampal ones are not).
    void note_eligibility(NodeId a, NodeId b) {
        if (!cfg_.cls_split) return;
        const std::uint64_t k = lane_key(a, b);
        const float e = (elig_.count(k) ? elig_[k] : 0.0f) + 1.0f;
        elig_[k] = e;
        if (e < cfg_.cortical_gate) return;
        auto& cl = cortical_[a];
        for (auto& l : cl)
            if (l.first == b) { l.second = std::min(4.0f, std::max(0.0f, l.second + 0.25f)); csr_dirty_ = true; return; }
        cl.emplace_back(b, 0.25f);
        csr_dirty_ = true;
    }
    void set_cls_split(bool on) { cfg_.cls_split = on; if (on) csr_dirty_ = true; }
    bool cls_split() const { return cfg_.cls_split; }
    void set_cortical_gate(float g) { cfg_.cortical_gate = g; }
    void set_cortical_gain(float g) { cfg_.cortical_gain = g; if (g != cfg_.cortical_gain) csr_dirty_ = true; }
    std::size_t cortical_lanes() const {
        std::size_t k = 0;
        for (const auto& kv : cortical_) k += kv.second.size();
        return k;
    }
    std::size_t eligibility_traces() const { return elig_.size(); }

    // -- M5 controls ---------------------------------------------------------
    // M9/encoder: configure the projection. Existing W is PRESERVED so a
    // trained encoder survives a re-learn; pass reset=true to discard.
    void set_encoder(std::size_t dim, std::size_t chans, bool reset = false) {
        if (dim == 0 || chans == 0) return;
        if (reset || enc_dim_ != dim || enc_chans_ != chans) {
            enc_dim_ = dim; enc_chans_ = chans;
            enc_W_.assign(chans * dim, 0.0f);
        }
        enc_y_.assign(chans, 0.0f);
    }
    void set_encoder_eta(float e) { enc_eta_ = e; }
    std::size_t enc_chans() const { return enc_chans_; }
    std::size_t enc_dim() const { return enc_dim_; }
    std::uint64_t enc_updates() const { return enc_updates_; }
    bool has_encoder() const { return !enc_W_.empty(); }

    // Oja update: W += eta * g * (x y^T - y y^T W).
    //   x  pre-synaptic  (encoder embedding, length d)
    //   y  post-synaptic (settled field, length C)
    //   g  ONE global scalar -- reused from the NTDA modulator, so exploration
    //      state already modulates the encoder. No spatial information, so no
    //      update ever needs a backward pass.
    void encoder_update(const std::vector<float>& x, const std::vector<float>& y, float g) {
        if (enc_W_.empty() || x.size() != enc_dim_ || y.size() != enc_chans_) return;
        const float lr = enc_eta_ * g;
        if (lr == 0.0f) return;
        for (std::size_t c = 0; c < enc_chans_; ++c) {
            const float yc = y[c];
            if (yc == 0.0f) continue;
            float* row = &enc_W_[c * enc_dim_];
            for (std::size_t d = 0; d < enc_dim_; ++d)
                row[d] += lr * (y[c] * x[d] - yc * y[c] * row[d]);
        }
        enc_y_ = y;
        ++enc_updates_;
    }

    // Project an encoder embedding into channel space and inject as energy.
    // This is the ONLY path by which an external representation touches the
    // field; the lane fabric is otherwise untouched, so typed/bigram lanes and
    // their order survive verbatim.
    std::vector<float> encoder_project(const std::vector<float>& x) const {
        std::vector<float> y(enc_chans_, 0.0f);
        if (enc_W_.empty() || x.size() != enc_dim_) return y;
        for (std::size_t c = 0; c < enc_chans_; ++c) {
            const float* row = &enc_W_[c * enc_dim_];
            float acc = 0.0f;
            for (std::size_t d = 0; d < enc_dim_; ++d) acc += row[d] * x[d];
            y[c] = acc;
        }
        return y;
    }
    void set_ntda(bool on) { cfg_.ntda_on = on; }
    bool ntda_on() const { return cfg_.ntda_on; }
    void set_ntda_theta(float t) { cfg_.ntda_theta = t; }
    void set_novelty(float n) { novelty_ = std::min(1.0f, std::max(0.0f, n)); }
    void set_conflict(float c) { conflict_ = std::min(1.0f, std::max(0.0f, c)); }
    float current_theta() const {
        if (!cfg_.ntda_on) return cfg_.ntda_theta;
        // dopamine lowers theta on novelty (exploration); cortisol raises it on
        // contradiction (exploitation). Both are multiplicative around 1.
        const float d = 1.0f - cfg_.ntda_dopamine * novelty_;
        const float c = 1.0f + cfg_.ntda_cortisol * conflict_;
        return std::max(0.05f, cfg_.ntda_theta * d * c);
    }
    std::uint64_t ntda_spikes() const { return ntda_spikes_; }
    std::uint64_t ntda_deaths() const { return ntda_deaths_; }
    std::size_t voltage_traces() const { return volt_.size(); }

    // M5 (a): accumulate voltage; fire and reset only at threshold. Returns true
    // when the lane actually bound. With ntda off this is always true, so the
    // pre-M5 path is byte-identical.
    bool ntda_gate(NodeId a, NodeId b) {
        if (!cfg_.ntda_on) return true;
        const std::uint64_t k = lane_key(a, b);
        const float th = current_theta();
        float& v = volt_[k];
        v += 1.0f;
        if (v < th) return false;
        v = 0.0f;              // refractory: cannot re-fire without re-stimulation
        ++ntda_spikes_;
        return true;
    }

    void add_lane(NodeId a, NodeId b, float w) {
        note_eligibility(a, b);
        if (!ntda_gate(a, b)) return;      // M5 (a): below threshold, no bind
        auto& lanes = out_[a];
        for (auto& l : lanes)
            if (l.first == b) { l.second = std::min(4.0f, std::max(0.0f, l.second + w)); csr_dirty_ = true; return; }
        if (w > 0.0f) {
            lanes.emplace_back(b, w);
            csr_dirty_ = true;
            if (lanes.size() > cfg_.lane_cap) {          // evict weakest lane
                auto weakest = std::min_element(lanes.begin(), lanes.end(),
                    [](const auto& x, const auto& y){ return x.second < y.second; });
                drop_lane_ctx(a, weakest->first);        // v3.2: signature dies with its lane
                lanes.erase(weakest);
            }
        }
    }

    // subtract weight from a lane; a lane whose weight hits zero dissolves
    void scale_lane(NodeId a, NodeId b, float dw) {
        auto it = out_.find(a);
        if (it == out_.end()) return;
        auto& lanes = it->second;
        for (std::size_t i = 0; i < lanes.size(); ++i)
            if (lanes[i].first == b) {
                lanes[i].second = std::max(0.0f, lanes[i].second + dw);
                if (lanes[i].second <= 0.0f) {
                    drop_lane_ctx(a, b);                 // v3.2: signature dies with its lane
                    lanes.erase(lanes.begin() + static_cast<std::ptrdiff_t>(i));
                }
                csr_dirty_ = true;
                return;
            }
    }

    // -- INJECTION ------------------------------------------------------------
    // Fresh tokens carry inject_energy each, damped by acoustic mass so that
    // frequent concepts cannot flood the field.
    void inject(const std::vector<std::string>& tokens, float energy = 0.0f) {
        if (energy <= 0.0f) energy = cfg_.inject_energy;
        for (const auto& t : tokens) {
            if (!has(t)) continue;                        // unknown concepts stay dark
            NodeId id = find(t);
            Node& n = nodes_[id];
            float mass_damp = 1.0f / std::sqrt(n.mass);   // heavy => less per-token
            n.energy += energy * mass_damp;
            n.salience = 1.0f;   // TSDA: spike to ceiling on touch
            present_.insert(id);   // v3.2: context-lane visibility (until next reset)
            state_hash_ = state_hash_ * 1099511628211ull
                        + static_cast<std::uint64_t>(fnv1a_hash(t));  // decision fingerprint
        }
    }

    // salience of a concept (SI integrator state; diagnostic + test surface)
    float node_salience(const std::string& c) const {
        auto it = index_.find(c);
        return it == index_.end() ? 0.0f : nodes_[it->second].salience;
    }

    // v3.9.0 accessors: intern-order label lookup (deterministic scans) and
    // the decision fingerprint (CEM variant draws, miller discipline).
    const std::string& label_of(NodeId id) const { return nodes_[id].concept; }
    std::uint64_t state_hash() const { return state_hash_; }

    // the source cap actually used by the last settle() (sampled value when
    // miller_window is on; otherwise always config().source_cap)
    float last_source_cap() const { return last_cap_; }

    // -- v3.8.0 (TASK 4) ENERGY SELF-VERIFICATION: the settle's own trace -----
    // EBT lineage, physics-native: a decision's confidence has a THIRD
    // measurement besides margin and calibration — whether the field had
    // finished MOVING when the readout looked at it. Every pass the energy
    // does two things: it DECAYS (uniform, 0.82 baseline — dissipates but
    // moves nothing between nodes) and it DIFFUSES (routes energy along
    // lanes — the actual computation). The trace records, for the last pass
    // executed:
    //   motion_rate — sum_i |e_new[i] - e_old[i]| / sum_i e_old[i]. The
    //                 decay baseline contributes ~0.18 of this; anything
    //                 measurably above is still-travelling structure.
    //   passes      — passes executed (== k_settle unless the eps break hit).
    // The eps break itself (dissipated < 1e-4) almost never fires under
    // decay-dominated dissipation, so motion_rate is the honest "did it
    // converge" signal; passes and both energies are disclosed for context.
    struct SettleTrace {
        int   passes = 0;
        int   k_requested = 0;
        float initial_energy = 0.0f;
        float final_energy = 0.0f;
        float motion_sum = 0.0f;      // last-pass sum |e_new - e_old|
        float motion_base = 0.0f;     // last-pass sum e_old (rate denominator)
        bool  eps_break = false;      // broke early via the eps relaxation
        float motion_rate() const {
            return motion_base > 0.0f ? motion_sum / motion_base : 0.0f;
        }
    };
    const SettleTrace& last_settle_trace() const { return last_trace_; }

    // -- v3.9.0 (MCPE §1 / CPME §2) MOMENTUM READOUT STATE --------------------
    // Emission should couple to nodes that RECEIVED diffusion, never to the
    // just-injected prompt tokens (the parrot effect). CPME's P_i is the
    // delta of a low-pass envelope; syfox's regime is a FIXED 8-pass burst
    // with decay 0.82, where every node's late trajectory is decay-dominated
    // (both a raw final-pass delta and a beta=0.9 envelope measured EVERY
    // bench candidate echo-suppressed). The faithful adaptation of the same
    // principle to the burst regime: momentum = the node's PEAK DIFFUSION
    // GAIN over passes 1..k,  max_p (E_p - E_{p-1}). An injected token only
    // decays after pass 0 => gain <= 0 => echo-suppressed. An anchor that
    // diffusion reached gained energy at some pass => eligible. Read-only
    // over the physics; off = zero cost, replay bit-identical.
    void set_momentum(bool on) { momentum_on_ = on; }
    bool momentum_on() const { return momentum_on_; }
    float node_momentum(NodeId id) const {
        if (!momentum_on_ || id >= max_gain_.size()) return 0.0f;
        return max_gain_[id];
    }

    // -- v3.9.0 (VGHT §2) COINCIDENCE GATE AT READOUT -------------------------
    // Cross-context bridges conduct only under simultaneous dual excitation:
    // a lane (u,v) traversed by the multi-hop readout walk contributes
    // w * C(u,v) with C = sigmoid(kappa*(E_u*E_v - theta_critical)),
    // kappa = 12, theta_critical = 0.05 — the paper's own constants. Either
    // endpoint dark => C ~ 0: spurious one-sided evidence cannot recruit a
    // distant subgraph. Pure readout arithmetic — the settled field is
    // measured, never rewritten; off = the walk is byte-identical to v3.8.
    void set_vght_gate(bool on) { vght_on_ = on; }
    bool vght_gate_on() const { return vght_on_; }
    void reset_vght_stats() const { vght_gated_sum_ = 0.0f; vght_walked_ = 0; }
    float vght_gated_frac() const {
        return vght_walked_ ? vght_gated_sum_ / static_cast<float>(vght_walked_) : 0.0f;
    }
    static float vght_conductance(float e_u, float e_v) {
        return 1.0f / (1.0f + std::exp(-12.0f * (e_u * e_v - 0.05f)));
    }

    // v3.8.0 (TASK 1) adaptive depth: k_settle is an existing CONFIG knob
    // ("max settle passes"), not a physics equation — adapting it per
    // decision scales THINKING DEPTH with problem difficulty while the
    // settle equations stay byte-identical. The Engine computes the
    // difficulty and sets the depth before each settle.
    void set_k_settle(int k) { cfg_.k_settle = std::max(1, k); }

    // v3.9.0 (P-CMA CEM shell): the source cap is an existing CONFIG knob —
    // planning samples it per variant (deterministically from the state hash,
    // miller_window discipline) and re-settles the SAME injected field.
    void set_source_cap(float c) { cfg_.source_cap = std::max(1.0f, c); }

    // v3.9.1 diagnostics sweep: setters for the remaining CONFIG knobs, so a
    // measured sweep (docs: physics constant sweep) does not require a
    // rebuilt binary per point. Defaults are untouched — cfg_ starts at the
    // frozen SI values and every setter is opt-in per process, so replay and
    // bench bit-identity hold unless a flag explicitly changes a knob
    // (test_reason-era contract: the default configuration IS the shipped
    // physics; knobs only scale/excite it, never rewrite the equations).
    //   decay / diffusion — settle dissipation & lane flow fractions
    //   hop_coupling      — readout-side lane contribution scaling
    //   learn_eta         — Hebbian bind rate (LEARN-side; affects fabrics
    //                        taught in the same process AFTER the set)
    //   lane_cap          — per-node lane count ceiling (LEARN-side eviction)
    void set_decay(float d) {
        cfg_.decay = std::min(0.999f, std::max(0.05f, d));
    }
    void set_diffusion(float d) {
        cfg_.diffusion = std::min(1.0f, std::max(0.0f, d));
    }
    void set_hop_coupling(float h) {
        cfg_.hop_coupling = std::min(2.0f, std::max(0.0f, h));
    }
    void set_learn_eta(float e) {
        cfg_.learn_eta = std::min(1.0f, std::max(0.0f, e));
    }
    void set_lane_cap(std::size_t cap) {
        cfg_.lane_cap = std::max<std::size_t>(8, cap);
    }
    float decay()        const { return cfg_.decay; }
    float diffusion()    const { return cfg_.diffusion; }
    float hop_coupling() const { return cfg_.hop_coupling; }
    float learn_eta()    const { return cfg_.learn_eta; }
    std::size_t lane_cap() const { return cfg_.lane_cap; }

    // v3.9.0 (HTR §2.2 discrete surrogate): participation ratio of the whole
    // energized field, R = (ΣE)² / (N·ΣE²) over nodes with E > 1e-7 — the
    // no-phase coherence proxy the ΔR gate compares the candidate set against.
    float field_participation_ratio() const {
        double s = 0.0, s2 = 0.0;
        std::size_t n = 0;
        for (const auto& nd : nodes_) {
            if (nd.energy <= 1e-7f) continue;
            s += nd.energy; s2 += static_cast<double>(nd.energy) * nd.energy;
            ++n;
        }
        if (n == 0 || s <= 0.0) return 0.0f;
        return static_cast<float>((s * s) / (static_cast<double>(n) * s2));
    }

    // -- v3.8.0 (TASK 2) SPECTRAL CHANNEL GATE --------------------------------
    // Multi-vector state: the SAME fabric settled under K different
    // lane-flow gates, one per leading distvec spectral axis. The gate for
    // lane (a,b) in channel k is exp(-(u_k(a)-u_k(b))^2 / (2 sigma_k^2)) —
    // lanes between nodes with similar k-th spectral coordinate carry full
    // flow, lanes between spectrally distant nodes are damped. Each channel
    // is one VIEW of the same Hebbian graph (same physics, same decay 0.82
    // / diffusion 0.45, energy conserved per pass exactly as always); the
    // readout superposes the K settled fields (mean). Requires a --distvec
    // fabric (the gates are its spectral coordinates); the gate multiplies
    // the context factor inside lane_flow_factor, and when INACTIVE the
    // factor math is bit-identical to v3.7.
    void set_channel_gate(std::vector<float> coords, float sigma) {
        channel_coords_ = std::move(coords);
        channel_sigma_ = sigma > 0.0f ? sigma : 1.0f;
    }
    void clear_channel_gate() { channel_coords_.clear(); channel_sigma_ = 0.0f; }
    bool channel_gate_active() const { return !channel_coords_.empty(); }

    // -- SETTLE ---------------------------------------------------------------
    // One pass = dissipative diffusion along lanes + source gating.
    // Energy is CONSERVED per pass up to decay: each node retains
    // (1 - diffusion) of its (decayed) energy and flows `diffusion` out along
    // its lanes, split proportionally to lane weight. The source cap GATES
    // propagation: only the top-`cap` active nodes send energy each pass —
    // ranked by energy by default, or by salience (motion history) when
    // salience_gating is on. Gated nodes keep their energy (readout-visible)
    // but stay silent as sources — bounded working set, nothing is
    // annihilated. Salience itself integrates every pass exactly as the SI
    // cavity physics does: s = tanh(s·decay + gain·|ΔE|); a node at rest is
    // an exact fixed point (tanh(0) = 0), mirroring SI's sparsity guard.
    void settle() {
        std::vector<float> next(nodes_.size(), 0.0f);
        std::vector<std::size_t> active;
        last_trace_ = SettleTrace{};
        last_trace_.k_requested = cfg_.k_settle;

        // M5: build/refresh the CSR mirror if the fabric moved. The build is
        // order-preserving (per-source vector copied verbatim), so switching
        // between map iteration and CSR iteration cannot change a single float.
        if (csr_dirty_) build_csr();

        // M5: parallel-settle buffers. Allocated once per call, reused across
        // passes. Determinism contract: static schedule partitions sources
        // into contiguous ascending chunks; per-thread scatter buffers are
        // combined in ascending thread order; therefore every next[i]
        // accumulates contributions in ascending SOURCE order — exactly the
        // sequential order. Bit-identical, test-verified.
        int nthreads = 1;
#if defined(_OPENMP)
        if (parallel_settle_) nthreads = omp_get_max_threads();
        if (nthreads > static_cast<int>(nodes_.size())) nthreads = nodes_.size() > 0 ? static_cast<int>(nodes_.size()) : 1;
#endif
        std::vector<std::vector<float>> thread_bufs(
            static_cast<std::size_t>(nthreads > 1 ? nthreads : 0),
            std::vector<float>(nodes_.size(), 0.0f));  // each buffer is n floats, zeroed

        // Live source cap (TSDA live_cap lineage): when miller_window is on,
        // the working set is drawn from [source_cap-4, source_cap] each
        // decision instead of being pinned — narrow focus some decisions,
        // broad others. SI samples a seeded mt19937 per tick; SyFox derives
        // the draw from the decision's state hash, so reproducibility holds:
        // same state -> same cap -> same field, bit for bit.
        float cap = cfg_.source_cap;
        if (cfg_.miller_window) {
            std::uint32_t hi = static_cast<std::uint32_t>(std::max(1.0f, cfg_.source_cap));
            std::uint32_t lo = (hi >= 5) ? hi - 4 : 1;
            cap = static_cast<float>(lo + (state_hash_ % (hi - lo + 1)));
        }
        last_cap_ = cap;
        const std::size_t cap_n = static_cast<std::size_t>(cap);
        const std::size_t n = nodes_.size();
        for (int pass = 0; pass < cfg_.k_settle; ++pass) {
            float total_before = 0.0f;
            for (const auto& nd : nodes_) total_before += nd.energy;
            if (pass == 0) last_trace_.initial_energy = total_before;
            if (total_before <= 0.0f) return;             // nothing to settle

            std::fill(next.begin(), next.end(), 0.0f);
            for (auto& tb : thread_bufs) std::fill(tb.begin(), tb.end(), 0.0f);
            active.clear();
            for (std::size_t i = 0; i < n; ++i)
                if (nodes_[i].energy > 1e-7f) active.push_back(i);

            // source gating: at most `cap` sources per pass (ties break by
            // vocabulary order — deterministic)
            if (active.size() > cap_n) {
                std::nth_element(active.begin(),
                                 active.begin() + static_cast<std::ptrdiff_t>(cap_n),
                                 active.end(), [&](std::size_t a, std::size_t b){
                                     return cfg_.salience_gating
                                         ? nodes_[a].salience > nodes_[b].salience
                                         : nodes_[a].energy   > nodes_[b].energy; });
                active.resize(cap_n);
            }
            std::vector<char> is_source(n, 0);
            for (std::size_t i : active) is_source[i] = 1;

            // -- diffusion pass (CSR slices; optional deterministic OMP) ----
#if defined(_OPENMP)
#pragma omp parallel if (nthreads > 1) num_threads(nthreads)
#endif
            {
                float* mybuf = nullptr;
#if defined(_OPENMP)
                const int tid = omp_get_thread_num();
                if (nthreads > 1) mybuf = thread_bufs[static_cast<std::size_t>(tid)].data();
#pragma omp for schedule(static)
#endif
                for (std::ptrdiff_t ii = 0; ii < static_cast<std::ptrdiff_t>(n); ++ii) {
                    const std::size_t i = static_cast<std::size_t>(ii);
                    const Node& nd = nodes_[i];
                    if (nd.energy <= 0.0f) continue;
                    float* scatter = mybuf ? mybuf : next.data();
                    if (!is_source[i]) { scatter[i] += nd.energy * cfg_.decay; continue; }  // gated: decay only
                    const std::size_t b = csr_off_[i], e = csr_off_[i + 1];
                    float out_w = 0.0f;
                    if (sem_active() && (!lane_ctx_.empty() || channel_gate_active())) {
                        // v3.2 Stage 2: context-sensitive lanes — each lane's
                        // EFFECTIVE weight is w * factor(f). The node's total
                        // outflow scales by (sum w*f / sum w): damped lanes
                        // carry less, the energy stays home. All factors 1.0
                        // reproduces the plain path bit-for-bit (same float
                        // sum order, 1.0f multiplies are exact).
                        // v3.8.0 TASK 2: the branch also runs when ONLY the
                        // spectral channel gate is active (no lane contexts).
                        float raw_w = 0.0f, eff_w = 0.0f;
                        for (std::size_t k = b; k < e; ++k) {
                            const float f = lane_flow_factor(i, csr_dst_[k]);
                            raw_w += csr_w_[k];
                            eff_w += csr_w_[k] * f;
                        }
                        if (raw_w > 0.0f && eff_w > 0.0f) {
                            const float out_frac = cfg_.diffusion * (eff_w / raw_w);
                            scatter[i] += nd.energy * cfg_.decay * (1.0f - out_frac);
                            const float flow = nd.energy * cfg_.decay * out_frac;
                            for (std::size_t k = b; k < e; ++k) {
                                const float f = lane_flow_factor(i, csr_dst_[k]);
                                if (f > 0.0f)
                                    scatter[csr_dst_[k]] += flow * (csr_w_[k] * f / eff_w);
                            }
                        } else {
                            scatter[i] += nd.energy * cfg_.decay;   // nothing flows
                        }
                        continue;
                    }
                    for (std::size_t k = b; k < e; ++k) out_w += csr_w_[k];
                    float retained = (out_w > 0.0f) ? (1.0f - cfg_.diffusion) : 1.0f;
                    scatter[i] += nd.energy * cfg_.decay * retained;   // damped self-retention
                    if (out_w > 0.0f) {
                        const float flow = nd.energy * cfg_.decay * cfg_.diffusion;
                        for (std::size_t k = b; k < e; ++k)
                            scatter[csr_dst_[k]] += flow * (csr_w_[k] / out_w);
                    }
                }
            }
            if (nthreads > 1) {                       // fixed thread-order combine
                for (const auto& tb : thread_bufs)
                    for (std::size_t i = 0; i < n; ++i) next[i] += tb[i];
            }
            // -- v3.2 SEMANTIC RESONANCE (Stage 4) ----------------------------
            // Sequential by design (determinism contract): sources — and only
            // sources, so the working-set cap stays meaningful — leak a small
            // fraction (sem_coupling) of their POST-diffusion energy along
            // their resonance edges, distributed proportionally to
            // cos * frequency-match. Energy is conserved exactly (the pool is
            // subtracted from the source, then distributed); decay applies
            // next pass as always. A node with no resonance edges is an
            // exact fixed point of this pass.
            if (sem_active()) {
                for (std::size_t i = 0; i < n; ++i) {
                    if (!is_source[i]) continue;
                    const std::size_t sb = sem_off_[i], se = sem_off_[i + 1];
                    if (sb == se) continue;
                    float wsum = 0.0f;
                    for (std::size_t k = sb; k < se; ++k)
                        wsum += sem_edges_[k].second * freq_match(i, sem_edges_[k].first);
                    if (wsum <= 0.0f) continue;
                    const float pool = next[i] * cfg_.sem_coupling;
                    next[i] -= pool;
                    for (std::size_t k = sb; k < se; ++k) {
                        const float fm = freq_match(i, sem_edges_[k].first);
                        next[sem_edges_[k].first] += pool * (sem_edges_[k].second * fm / wsum);
                    }
                }
            }
            // SI cavity salience integrator: s = tanh(s·decay + gain·|ΔE|).
            // Motion proxy is how much energy moved through the node this
            // pass. Updated before the energies are committed, so the next
            // pass's source ranking sees it. Element-wise (no accumulation):
            // parallel-safe without touching the determinism contract.
            // v3.8.0 TASK 4: the SAME pre-commit ΔE accumulates (serially,
            // after the parallel salience pass) into the settle trace —
            // relative motion of the field this pass. Must read BEFORE the
            // commit, where next and energy still differ.
            {
                float motion = 0.0f;
                for (std::size_t i = 0; i < n; ++i)
                    motion += std::fabs(next[i] - nodes_[i].energy);
                last_trace_.motion_sum = motion;
            }
            // v3.9.0 MCPE/CPME: track each node's PEAK diffusion gain across
            // passes (next[i] - energy[i] pre-commit; decay makes pure-echo
            // nodes non-positive, diffusion receivers positive).
            if (momentum_on_) {
                if (max_gain_.size() != n) max_gain_.assign(n, 0.0f);
                for (std::size_t i = 0; i < n; ++i) {
                    const float gain = next[i] - nodes_[i].energy;
                    if (gain > max_gain_[i]) max_gain_[i] = gain;
                }
            }
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) if (nthreads > 1)
#endif
            for (std::ptrdiff_t ii = 0; ii < static_cast<std::ptrdiff_t>(n); ++ii) {
                const std::size_t i = static_cast<std::size_t>(ii);
                if (nodes_[i].energy <= 0.0f && next[i] <= 0.0f) continue;
                const float motion = std::fabs(next[i] - nodes_[i].energy);
                nodes_[i].salience = std::tanh(nodes_[i].salience * cfg_.salience_decay
                                             + cfg_.salience_gain * motion);
            }
#if defined(_OPENMP)
#pragma omp parallel for schedule(static) if (nthreads > 1)
#endif
            for (std::ptrdiff_t ii = 0; ii < static_cast<std::ptrdiff_t>(n); ++ii)
                nodes_[ii].energy = next[ii];

            float total_after = 0.0f;
            for (const auto& nd : nodes_) total_after += nd.energy;
            last_trace_.passes = pass + 1;
            last_trace_.motion_base = total_before;
            last_trace_.final_energy = total_after;
            if (total_before - total_after < cfg_.eps) {
                last_trace_.eps_break = true;
                break;  // relaxed
            }
        }
    }

    // -- READOUT (resonance sweep, two lanes) -----------------------------------
    // Both lanes are field measurements — no weights fitted, no classes exist.
    //
    //  SPECIFIC — "mean lane support": rewards precise coupling. A node bound
    //    diffusely to half the vocabulary must not read hotter than a node
    //    bound exactly to this state's concepts. Used for Choice/Score where
    //    the engine must pick THE best-matching option.
    //
    //  SUPPORT — "active lane mass": rewards breadth of corroboration. A
    //    statement (Noul) may be supported by many partial lanes; what matters
    //    is how much lane-connected energy the field offers it. Used for Noul.
    enum class ReadoutMode { Specific, Support };

    float readout(const std::vector<std::string>& probe,
                  ReadoutMode mode = ReadoutMode::Specific) const {
        float r = 0.0f;
        for (const auto& t : probe) {
            auto it = index_.find(t);
            if (it == index_.end()) continue;
            const Node& n = nodes_[it->second];
            r += n.energy;
            const auto lane_it = out_.find(it->second);
            if (lane_it == out_.end() || lane_it->second.empty()) continue;
            if (mode == ReadoutMode::Specific) {
                if (hop_depth_ <= 1) {
                    // legacy single-hop path — bit-identical to pre-v3.4
                    float wsum = 0.0f, flow = 0.0f;
                    for (const auto& lane : lane_it->second) {
                        wsum += lane.second;
                        flow += lane.second * nodes_[lane.first].energy;
                    }
                    if (wsum > 0.0f) r += cfg_.hop_coupling * flow / wsum;
                } else {
                    // v3.4 multi-hop walk: BFS over lanes from this anchor.
                    // Level h contributes hop_coupling/sqrt(h+1) of the
                    // weight-averaged energy there. Deterministic: candidate
                    // nodes collect in ascending NodeId order (std::map),
                    // width-capped, never revisited.
                    std::vector<char> visited(nodes_.size(), 0);
                    visited[it->second] = 1;
                    std::vector<NodeId> frontier{it->second};
                    // v3.9.0 VGHT coincidence gate: when on, each traversed
                    // lane's weight and flow are scaled by C = σ(12(E_u·E_v −
                    // 0.05)); per-pair bookkeeping feeds vght_gated_frac().
                    const bool vght = vght_on_;
                    if (vght) reset_vght_stats();
                    for (int hop = 1; hop <= hop_depth_ && !frontier.empty(); ++hop) {
                        std::map<NodeId, float> next_w;
                        if (vght) {
                            for (NodeId u : frontier) {
                                const auto uit = out_.find(u);
                                if (uit == out_.end()) continue;
                                const float eu = nodes_[u].energy;
                                for (const auto& lane : uit->second) {
                                    if (visited[lane.first]) continue;
                                    const float ev = nodes_[lane.first].energy;
                                    const float c = vght_conductance(eu, ev);
                                    next_w[lane.first] += lane.second * c;
                                    ++vght_walked_;
                                    if (c >= 0.5f) vght_gated_sum_ += 1.0f;
                                }
                            }
                        } else {
                            for (NodeId u : frontier) {
                                const auto uit = out_.find(u);
                                if (uit == out_.end()) continue;
                                for (const auto& lane : uit->second)
                                    if (!visited[lane.first])
                                        next_w[lane.first] += lane.second;
                            }
                        }
                        if (next_w.empty()) break;
                        std::vector<NodeId> level;
                        float wsum = 0.0f, flow = 0.0f;
                        for (const auto& kv : next_w) {
                            if (level.size() >= kHopWidthCap) break;
                            level.push_back(kv.first);
                            wsum += kv.second;
                            flow += kv.second * nodes_[kv.first].energy;
                        }
                        for (NodeId v : level) visited[v] = 1;
                        const float damp = cfg_.hop_coupling
                            / std::sqrt(static_cast<float>(hop + 1));
                        if (wsum > 0.0f) r += damp * flow / wsum;
                        frontier.swap(level);
                    }
                }
            } else {                                      // Support: active lane mass
                for (const auto& lane : lane_it->second)
                    r += cfg_.hop_coupling * lane.second * nodes_[lane.first].energy;
            }
            // v3.2 Stage 4: semantic-neighbour term — energy resting on
            // semantically similar nodes counts for this probe even where no
            // lane connects them. Weighted by the edge cosine; zero when the
            // fabric carries no semantic field (pre-v3.2 models, or --no-semantics).
            if (sem_active()) {
                const std::size_t sb = sem_off_[it->second], se = sem_off_[it->second + 1];
                for (std::size_t k = sb; k < se; ++k)
                    r += cfg_.sem_hop * sem_edges_[k].second * nodes_[sem_edges_[k].first].energy;
            }
        }
        return r;
    }

    float total_energy() const {
        float s = 0.0f;
        for (const auto& n : nodes_) s += n.energy;
        return s;
    }

    // v3.4 — decision-layer field composition (question-context gate): blend
    // the settled state field with an ALREADY-SETTLED context field,
    // e = (1-alpha)*e + alpha*ctx. Settle() physics is untouched — this
    // composes two settled fields at the decision layer, so the question can
    // carry part of the readout field without rewriting the dynamics.
    void blend_field(const std::vector<float>& ctx, float alpha) {
        const std::size_t n = std::min(nodes_.size(), ctx.size());
        for (std::size_t i = 0; i < n; ++i)
            nodes_[i].energy = (1.0f - alpha) * nodes_[i].energy + alpha * ctx[i];
    }

    // v3.6.0 — BED §8 perturbation-contrast support: snapshot the settled
    // field (energies only; salience/present_ belong to the composition, not
    // the measurement) and restore it after a diagnostic re-settle. The
    // perturbation check settles a structure-broken copy of the state in a
    // SCRATCH pass, measures its margins, then restores the real field so the
    // readout that follows measures exactly what the un-checked decide()
    // would have measured. Deterministic; no physics touched.
    std::vector<float> snapshot_field() const {
        std::vector<float> e(nodes_.size());
        for (std::size_t i = 0; i < nodes_.size(); ++i) e[i] = nodes_[i].energy;
        return e;
    }
    void restore_field(const std::vector<float>& e) {
        const std::size_t n = std::min(nodes_.size(), e.size());
        for (std::size_t i = 0; i < n; ++i) nodes_[i].energy = e[i];
    }

    float node_energy(NodeId id) const {
        return id < nodes_.size() ? nodes_[id].energy : 0.0f;
    }

    // acoustic mass (occurrence count); derivation uses the same 1/sqrt(mass)
    // damping law the field applies at its use sites, so hub-mediated
    // composition is suppressed exactly like hub-mediated energy flow
    float node_mass(NodeId id) const {
        return id < nodes_.size() ? nodes_[id].mass : 1.0f;
    }

    bool silent() const { return total_energy() < cfg_.silence_floor; }

    // -- HONEST SILENCE reset ---------------------------------------------------
    void reset_field() {
        for (auto& n : nodes_) { n.energy = 0.0f; n.salience = 0.0f; }
        present_.clear();
        state_hash_ = 0;
        // v3.9.0: the momentum state is per-decision — a fresh field must not
        // inherit the previous decision's trajectory.
        max_gain_.clear();
    }

    // -- ANTI-HEBBIAN WEAKENING --------------------------------------------------
    // Counter-evidence: unbind lanes between active concepts and an outcome
    // (used for Noul "false" lessons so disconfirming routes fade).
    void weaken(const std::vector<std::string>& active_state,
                const std::vector<std::string>& outcome, float eta_scale = 1.0f) {
        for (const auto& ta : active_state) {
            if (!has(ta)) continue;
            for (const auto& tb : outcome) {
                if (!has(tb)) continue;
                NodeId a = find(ta), b = find(tb);
                scale_lane(a, b, -cfg_.learn_eta * eta_scale);
                scale_lane(b, a, -cfg_.learn_eta * eta_scale);
            }
        }
    }

    // -- HEBBIAN LESSON ----------------------------------------------------------
    // Strengthen lanes between currently-active state concepts and the
    // concepts of a rewarded outcome; decay every lane slightly so unused
    // routes fade. This is how SyFox learns: substrate rewiring, not fitting.
    void hebbian_lesson(const std::vector<std::string>& active_state,
                        const std::vector<std::string>& outcome,
                        float eta_scale = 1.0f) {
        std::vector<NodeId> s_ids, o_ids;
        for (const auto& t : active_state)
            if (has(t)) s_ids.push_back(find(t));
        for (const auto& t : outcome) {
            intern(t);                                    // outcome concepts join the field
            o_ids.push_back(find(t));
        }
        for (NodeId a : s_ids)
            for (NodeId b : o_ids) {
                float wa = 1.0f / std::sqrt(nodes_[a].mass);
                float wb = 1.0f / std::sqrt(nodes_[b].mass);
                // v3.6.0 (CMD §2.2 / MIMS §6): rarity weighting on the STATE
                // side only. Outcome anchors keep full weight — their df is
                // per-class, and a class label seen often is not a stopword.
                // idf_factor is 1.0 whenever the df table is absent (pre-v3.6
                // fabrics), so existing lanes are reproduced exactly.
                // v3.9.0 (CPME-II §1, discrete surrogate): phase-consistent
                // consolidation — the order-evidence factor Phi in (−0.2, 1]
                // multiplies the same bind (empty table / flag off => 1.0,
                // replay-identical). Chirality-consistent co-activity
                // consolidates at full strength; reversed co-activity cancels.
                const float ord = order_factor(a, b);
                if (ord <= 0.0f) continue;   // anti-phase: mild depression (no bind)
                bind(a, b, cfg_.learn_eta * eta_scale * (wa + wb) * idf_factor(a) * ord);
            }
        enforce_lane_decay();
    }

    // -- v3.9.0 ORDER EVIDENCE (CPME-II port) ---------------------------------
    void set_order_consolidation(bool on) { order_on_ = on; }
    bool order_on() const { return order_on_; }
    void note_order(NodeId a, NodeId b, bool consistent) {
        auto& s = order_stats_[lane_key(a, b)];
        if (consistent) ++s.first;
        ++s.second;
    }
    float order_factor(NodeId a, NodeId b) const {
        if (!order_on_) return 1.0f;
        const auto it = order_stats_.find(lane_key(a, b));
        if (it == order_stats_.end() || it->second.second == 0) return 1.0f;
        const float frac = static_cast<float>(it->second.first)
                         / static_cast<float>(it->second.second);
        const float phi = 0.4f + 0.6f * (2.0f * frac - 1.0f);   // alpha_Phi = 0.6
        return std::max(0.0f, phi);
    }
    std::size_t order_table_size() const { return order_stats_.size(); }

    // -- LANE EVIDENCE LEDGER (v3 Milestone 3; bookkeeping, not physics) ------
    // record_support / record_counter: called by the ENGINE around its teach
    // events. seq = the engine's monotonic teach counter; context = engine tag.
    void record_support(NodeId a, NodeId b, std::uint64_t seq, const std::string& context) {
        LaneEvidence& e = lane_evidence_[lane_key(a, b)];
        ++e.support_events;
        e.last_seq = seq;
        if (e.first_seq == 0) e.first_seq = seq;
        e.context = context;
    }
    void record_counter(NodeId a, NodeId b, std::uint64_t seq, const std::string& context) {
        LaneEvidence& e = lane_evidence_[lane_key(a, b)];
        ++e.counter_events;
        e.last_seq = seq;
        if (e.first_seq == 0) e.first_seq = seq;
        e.context = context;
    }
    bool evidence_of(NodeId a, NodeId b, LaneEvidence& out) const {
        auto it = lane_evidence_.find(lane_key(a, b));
        if (it == lane_evidence_.end()) return false;
        out = it->second;
        return true;
    }
    std::size_t evidence_count() const { return lane_evidence_.size(); }

    // =====================================================================
    // v3.2 SEMANTIC FIELD API (deterministic; see the header comment block)
    // =====================================================================

    // Stage 2 learn-side accumulation: ctx_word co-occurred with lane (a,b)
    // in the lesson that laid it. Counts live in a scratch table; only
    // finalize_contexts() (called at save) promotes surviving counts into
    // the lane's signature. Bookkeeping — nothing here touches the field.
    void add_ctx_support(NodeId a, NodeId b, NodeId ctx_word) {
        if (a == b || ctx_word == a || ctx_word == b) return;
        lane_ctx_acc_[lane_key(a, b)][ctx_word] += 1.0f;
    }

    // Promote accumulated co-occurrence counts into per-lane signatures.
    // Deterministic pruning over node-id-ascending iteration:
    //   required  — top-4 context words with count >= 2 for this lane
    //   forbidden — up to 2 words that are STRONG required context of the
    //               same source node's lanes into a DIFFERENT target family
    //               (count >= 3 there, absent-or-weak here). This is the
    //               cross-class diff: the words that discriminate WHERE
    //               this source node's energy should go.
    // Idempotent: the scratch table is consumed; calling twice is a no-op.
    void finalize_contexts() {
        if (lane_ctx_acc_.empty()) return;
        // gather per-source: which targets exist and each target's top context
        for (auto& kv : lane_ctx_acc_) {
            const std::uint64_t key = kv.first;
            const NodeId a = static_cast<NodeId>(key >> 32);
            const NodeId b = static_cast<NodeId>(key & 0xffffffffull);
            // required: top-4 by (count, then node id ascending for ties)
            std::vector<std::pair<float, NodeId>> cand;
            for (const auto& cv : kv.second)
                if (cv.second >= 2.0f) cand.emplace_back(cv.second, cv.first);
            if (cand.empty()) continue;
            std::sort(cand.begin(), cand.end(), [](const auto& x, const auto& y) {
                return x.first != y.first ? x.first > y.first : x.second < y.second;
            });
            LaneCtx sig;
            for (std::size_t i = 0; i < cand.size() && i < 4; ++i)
                sig.required.push_back(cand[i].second);
            lane_ctx_[key] = std::move(sig);
            (void)a; (void)b;
        }
        // forbidden pass: for lane (a,b), words that dominate a sibling lane
        // (a,b') with b' != b and are not in this lane's required set.
        std::unordered_map<std::uint64_t, std::vector<std::uint64_t>> by_source;
        for (const auto& kv : lane_ctx_) {
            const NodeId a = static_cast<NodeId>(kv.first >> 32);
            by_source[a].push_back(kv.first);
        }
        for (auto& sv : by_source) {
            auto& keys = sv.second;
            std::sort(keys.begin(), keys.end());
            for (std::size_t i = 0; i < keys.size(); ++i) {
                // collect the strongest context words of the SIBLING lanes
                std::unordered_map<NodeId, float> sibling_ctx;
                for (std::size_t j = 0; j < keys.size(); ++j) {
                    if (i == j) continue;
                    auto it = lane_ctx_acc_.find(keys[j]);
                    if (it == lane_ctx_acc_.end()) continue;   // consumed
                    for (const auto& cv : it->second)
                        if (cv.second >= 3.0f) sibling_ctx[cv.first] += cv.second;
                }
                if (sibling_ctx.empty()) continue;
                // not already required here
                const LaneCtx& mine = lane_ctx_[keys[i]];
                auto is_req = [&](NodeId w) {
                    for (NodeId r : mine.required) if (r == w) return true;
                    return false;
                };
                std::vector<std::pair<float, NodeId>> fc;
                for (const auto& cv : sibling_ctx)
                    if (!is_req(cv.first)) fc.emplace_back(cv.second, cv.first);
                if (fc.empty()) continue;
                std::sort(fc.begin(), fc.end(), [](const auto& x, const auto& y) {
                    return x.first != y.first ? x.first > y.first : x.second < y.second;
                });
                for (std::size_t k = 0; k < fc.size() && k < 2; ++k)
                    lane_ctx_[keys[i]].forbidden.push_back(fc[k].second);
            }
        }
        lane_ctx_acc_.clear();
    }

    bool has_semantics() const { return !semvecs_.empty(); }
    void set_semantics(bool on) { sem_enabled_ = on; }

    // v3.4 — multi-hop readout walk depth. 1 (default) = the legacy single-hop
    // readout, bit-identical to every earlier version (replay contract).
    // Depths 2..8 walk lanes BFS-style from each probe anchor with per-hop
    // damping hop_coupling/sqrt(hop+1): hop 1 x0.707, hop 2 x0.577, hop 3
    // x0.500, hop 4 x0.447. Levels are visited in ascending NodeId order and
    // capped at 64 nodes per level, so the walk is deterministic and
    // bounded. Readout-layer only: settle() physics untouched.
    void set_hop_depth(int d) { hop_depth_ = std::max(1, std::min(8, d)); }
    int hop_depth() const { return hop_depth_; }
    std::size_t resonance_edge_count() const { return sem_edge_count_; }
    std::size_t lane_context_count() const { return lane_ctx_.size(); }

    // ---------------------------------------------------------------------
    // v3.7.0 DISTRIBUTIONAL FIELD (distvec.hpp, Move 1 of the limit analysis):
    // dense PPMI+SVD vectors over the lane co-occurrence structure, plus the
    // public install surface the builder uses. The vectors REPLACE NOTHING:
    // semvecs_ (trigram+grounding) still feeds omega/freq_match; the dense
    // vectors re-select the resonance EDGE SET only. Absent (dims 0) on every
    // pre-v3.7 fabric — bit-identical replay.
    void set_distvecs(std::vector<float>&& v, int dims) {
        distvecs_ = std::move(v);
        dist_dims_ = (distvecs_.size() == static_cast<std::size_t>(dims) * node_count()
                      && dims > 0) ? dims : 0;
        if (dist_dims_ == 0) distvecs_.clear();
    }
    bool has_distvecs() const { return dist_dims_ > 0; }
    int dist_dims() const { return dist_dims_; }
    const std::vector<float>& distvecs() const { return distvecs_; }

    // install a resonance-edge CSR (distvec.hpp edge selection; same storage
    // and settle role as the SEM4 tail edges built by build_semantics)
    void set_semantic_edges(std::vector<std::pair<NodeId, float>>&& edges,
                            std::vector<std::size_t>&& offs) {
        if (offs.size() != nodes_.size() + 1) return;   // shape guard
        sem_edges_ = std::move(edges);
        sem_off_ = std::move(offs);
        sem_edge_count_ = sem_edges_.size();
    }

    // semantic vector of a node (empty vector when semantics absent)
    std::vector<float> sem_vector(NodeId id) const {
        std::vector<float> v;
        if (id >= node_count() || semvecs_.empty()) return v;
        v.assign(semvecs_.begin() + static_cast<std::ptrdiff_t>(id * SEM_DIMS),
                 semvecs_.begin() + static_cast<std::ptrdiff_t>((id + 1) * SEM_DIMS));
        return v;
    }

    // Stage 1: omega_semantic — scalar projection of the semantic vector onto
    // one FIXED deterministic direction, mapped to [0,1]. Fixed point of the
    // frequency axis: same string -> same omega, always.
    float semantic_freq(NodeId id) const {
        if (id >= node_count() || semvecs_.empty()) return 0.0f;
        float d = 0.0f;
        for (std::size_t k = 0; k < SEM_DIMS; ++k)
            d += semvecs_[id * SEM_DIMS + k] * omega_dir()[k];
        return std::max(0.0f, std::min(1.0f, 0.5f + 0.5f * d));
    }

    // Option 2 typed-lane read: cos >= 0.82 synonym-grade kin,
    //             cos >= 0.66 related; 0 = no semantic relation measured.
    static const char* semantic_type(float cos) {
        if (cos >= 0.82f) return "synonym";
        if (cos >= 0.66f) return "related";
        if (cos > 0.0f)   return "resonance";
        return "none";
    }

    // resonance neighbours of a node (id, cosine) — empty when no semantics
    const std::vector<std::pair<NodeId, float>>* resonance_of(NodeId id) const {
        if (semvecs_.empty() || id >= node_count()) return nullptr;
        static const std::vector<std::pair<NodeId, float>> kEmpty;
        const std::size_t b = sem_off_[id], e = sem_off_[id + 1];
        if (b == e) return &kEmpty;
        // slices live in sem_edges_; return pointer into it (stable between builds)
        static thread_local std::vector<std::pair<NodeId, float>> out;
        out.assign(sem_edges_.begin() + static_cast<std::ptrdiff_t>(b),
                   sem_edges_.begin() + static_cast<std::ptrdiff_t>(e));
        return &out;
    }

    const LaneCtx* lane_context(NodeId a, NodeId b) const {
        auto it = lane_ctx_.find(lane_key(a, b));
        return it == lane_ctx_.end() ? nullptr : &it->second;
    }

    // Build the semantic field over the CURRENT fabric. Deterministic.
    // Expensive once (O(n^2 * dims) neighbour search, OMP-parallel with a
    // fixed per-row reduction); stored in the model from then on.
    void build_semantics(int threads = 0);

    // -- persistence (binary, deterministic) -------------------------------------
    // v3.2.1 BUGFIX — save() used to iterate unordered_map members (out_,
    // derived_gen_, lane_evidence_, lane_ctx_) directly, so every load->save
    // cycle permuted the file (measured: three different md5s over two
    // rebuild passes on one fabric; the lane MULTISET is identical, the ORDER
    // is not). Consequences: model dirs were never byte-stable, and the
    // settled field's floating-point sum order drifted across save cycles.
    // Fix: every map-backed section is written in SORTED KEY order. Load
    // rebuilds out_ in file order, so iteration order is a pure function of
    // the file => one sorted save converges the format and every later
    // roundtrip is bit-identical. Same lanes, same weights: no physics change.
    void save(const std::string& path) const {
        std::ofstream f(path, std::ios::binary);
        std::uint32_t n = static_cast<std::uint32_t>(nodes_.size());
        f.write(reinterpret_cast<const char*>(&n), 4);
        for (const auto& node : nodes_) {
            std::uint32_t len = static_cast<std::uint32_t>(node.concept.size());
            f.write(reinterpret_cast<const char*>(&len), 4);
            f.write(node.concept.data(), len);
            f.write(reinterpret_cast<const char*>(&node.mass), 4);
        }
        std::uint32_t lanes = 0;
        for (const auto& kv : out_) lanes += static_cast<std::uint32_t>(kv.second.size());
        f.write(reinterpret_cast<const char*>(&lanes), 4);
        std::vector<const std::pair<const NodeId, std::vector<std::pair<NodeId, float>>>*> sorted_out;
        sorted_out.reserve(out_.size());
        for (const auto& kv : out_) sorted_out.push_back(&kv);
        std::sort(sorted_out.begin(), sorted_out.end(),
                  [](const auto* x, const auto* y) { return x->first < y->first; });
        for (const auto* kv : sorted_out)
            for (const auto& l : kv->second) {
                f.write(reinterpret_cast<const char*>(&kv->first), 4);
                f.write(reinterpret_cast<const char*>(&l.first), 4);
                f.write(reinterpret_cast<const char*>(&l.second), 4);
            }
        // v2 provenance tail — old (v1) files simply end here; load() detects EOF.
        std::vector<std::pair<std::uint64_t, std::uint32_t>> sorted_derived;
        sorted_derived.reserve(derived_gen_.size());
        for (const auto& kv : derived_gen_) sorted_derived.emplace_back(kv.first, kv.second);
        std::sort(sorted_derived.begin(), sorted_derived.end(),
                  [](const auto& x, const auto& y) { return x.first < y.first; });
        std::uint32_t dc = static_cast<std::uint32_t>(sorted_derived.size());
        f.write(reinterpret_cast<const char*>(&dc), 4);
        for (const auto& kv : sorted_derived) {
            std::uint32_t a = static_cast<std::uint32_t>(kv.first >> 32);
            std::uint32_t b = static_cast<std::uint32_t>(kv.first & 0xffffffffull);
            f.write(reinterpret_cast<const char*>(&a), 4);
            f.write(reinterpret_cast<const char*>(&b), 4);
            f.write(reinterpret_cast<const char*>(&kv.second), 4);
        }
        // v3 evidence tail — v2 files end here; load() detects EOF.
        std::vector<std::uint64_t> sorted_ev;
        sorted_ev.reserve(lane_evidence_.size());
        for (const auto& kv : lane_evidence_) sorted_ev.push_back(kv.first);
        std::sort(sorted_ev.begin(), sorted_ev.end());
        std::uint32_t ec = static_cast<std::uint32_t>(sorted_ev.size());
        f.write(reinterpret_cast<const char*>(&ec), 4);
        for (const auto& key : sorted_ev) {
            const auto& kv = *lane_evidence_.find(key);
            std::uint32_t a = static_cast<std::uint32_t>(kv.first >> 32);
            std::uint32_t b = static_cast<std::uint32_t>(kv.first & 0xffffffffull);
            f.write(reinterpret_cast<const char*>(&a), 4);
            f.write(reinterpret_cast<const char*>(&b), 4);
            f.write(reinterpret_cast<const char*>(&kv.second.support_events), 4);
            f.write(reinterpret_cast<const char*>(&kv.second.counter_events), 4);
            f.write(reinterpret_cast<const char*>(&kv.second.first_seq), 8);
            f.write(reinterpret_cast<const char*>(&kv.second.last_seq), 8);
            std::uint32_t clen = static_cast<std::uint32_t>(kv.second.context.size());
            f.write(reinterpret_cast<const char*>(&clen), 4);
            if (clen) f.write(kv.second.context.data(), clen);
        }
        // v3.2 SEMANTIC TAIL — vectors + resonance edges + lane contexts.
        // Magic-guarded: a truncated/garbage read leaves semantics off and the
        // rest of the model fully usable (pre-v3.2 files simply end here).
        const std::uint32_t sem_magic = 0x53454D34u;             // 'SEM4'
        f.write(reinterpret_cast<const char*>(&sem_magic), 4);
        const std::uint32_t has_sem = semvecs_.empty() ? 0u : 1u;
        f.write(reinterpret_cast<const char*>(&has_sem), 4);
        if (has_sem) {
            f.write(reinterpret_cast<const char*>(semvecs_.data()),
                    static_cast<std::streamsize>(semvecs_.size() * sizeof(float)));
            const std::size_t n = nodes_.size();
            std::uint32_t total = static_cast<std::uint32_t>(sem_edges_.size());
            f.write(reinterpret_cast<const char*>(&total), 4);
            for (std::size_t i = 0; i <= n; ++i) {
                std::uint32_t off = static_cast<std::uint32_t>(sem_off_[i]);
                f.write(reinterpret_cast<const char*>(&off), 4);
            }
            for (const auto& e : sem_edges_) {
                f.write(reinterpret_cast<const char*>(&e.first), 4);
                f.write(reinterpret_cast<const char*>(&e.second), 4);
            }
        }
        std::vector<std::uint64_t> sorted_ctx;
        sorted_ctx.reserve(lane_ctx_.size());
        for (const auto& kv : lane_ctx_) sorted_ctx.push_back(kv.first);
        std::sort(sorted_ctx.begin(), sorted_ctx.end());
        const std::uint32_t ctxc = static_cast<std::uint32_t>(sorted_ctx.size());
        f.write(reinterpret_cast<const char*>(&ctxc), 4);
        for (const auto& key : sorted_ctx) {
            const auto& kv = *lane_ctx_.find(key);
            f.write(reinterpret_cast<const char*>(&kv.first), 8);
            std::uint32_t rq = static_cast<std::uint32_t>(kv.second.required.size());
            f.write(reinterpret_cast<const char*>(&rq), 4);
            for (NodeId w : kv.second.required) f.write(reinterpret_cast<const char*>(&w), 4);
            std::uint32_t fb = static_cast<std::uint32_t>(kv.second.forbidden.size());
            f.write(reinterpret_cast<const char*>(&fb), 4);
            for (NodeId w : kv.second.forbidden) f.write(reinterpret_cast<const char*>(&w), 4);
        }
        // v3.6.0 IDF TAIL — document-frequency table + lesson counter (CMD §2.2
        // / MIMS §6 rarity grounding). Magic-guarded 'IDF5' after the context
        // block: pre-v3.6 files end here and load() leaves df_ empty, which
        // makes every idf_factor 1.0 — the replay contract, bit for bit.
        const std::uint32_t idf_magic = 0x49444635u;             // 'IDF5'
        f.write(reinterpret_cast<const char*>(&idf_magic), 4);
        const std::uint64_t lesson_count = df_lessons_;
        f.write(reinterpret_cast<const char*>(&lesson_count), 8);
        std::vector<std::pair<std::uint32_t, std::uint32_t>> sorted_df;
        sorted_df.reserve(df_.size());
        for (const auto& kv : df_) sorted_df.emplace_back(kv.first, kv.second);
        std::sort(sorted_df.begin(), sorted_df.end(),
                  [](const auto& x, const auto& y) { return x.first < y.first; });
        const std::uint32_t dfc = static_cast<std::uint32_t>(sorted_df.size());
        f.write(reinterpret_cast<const char*>(&dfc), 4);
        for (const auto& kv : sorted_df) {
            f.write(reinterpret_cast<const char*>(&kv.first), 4);
            f.write(reinterpret_cast<const char*>(&kv.second), 4);
        }
        // v3.7.0 DSTV TAIL — dense PPMI+SVD vectors (distvec.hpp). Magic-guarded
        // after the IDF block: pre-v3.7 files end here and load() leaves the
        // distributional field absent — the replay contract, bit for bit.
        const std::uint32_t dstv_magic = 0x44535456u;            // 'DSTV'
        f.write(reinterpret_cast<const char*>(&dstv_magic), 4);
        const std::uint32_t ddims = static_cast<std::uint32_t>(dist_dims_);
        f.write(reinterpret_cast<const char*>(&ddims), 4);
        if (ddims > 0) {
            f.write(reinterpret_cast<const char*>(distvecs_.data()),
                    static_cast<std::streamsize>(distvecs_.size() * sizeof(float)));
        }
        // v3.9.0 ORD1 TAIL — per-lane order evidence (CPME-II consolidation).
        // Magic-guarded after DSTV: pre-v3.9 files end here and load() leaves
        // order_stats_ empty, so every order_factor is 1.0 — replay contract.
        const std::uint32_t ord_magic = 0x4F524431u;             // 'ORD1'
        f.write(reinterpret_cast<const char*>(&ord_magic), 4);
        std::vector<std::tuple<std::uint32_t, std::uint32_t, std::uint32_t>> sorted_ord;
        sorted_ord.reserve(order_stats_.size());
        for (const auto& kv : order_stats_)
            sorted_ord.emplace_back(static_cast<std::uint32_t>(kv.first >> 32),
                                    static_cast<std::uint32_t>(kv.first & 0xffffffffull),
                                    0);
        std::sort(sorted_ord.begin(), sorted_ord.end());
        const std::uint32_t orc = static_cast<std::uint32_t>(sorted_ord.size());
        f.write(reinterpret_cast<const char*>(&orc), 4);
        for (const auto& t : sorted_ord) {
            const auto it = order_stats_.find(
                lane_key(static_cast<NodeId>(std::get<0>(t)),
                         static_cast<NodeId>(std::get<1>(t))));
            const std::uint32_t cons = it->second.first, tot = it->second.second;
            f.write(reinterpret_cast<const char*>(&std::get<0>(t)), 4);
            f.write(reinterpret_cast<const char*>(&std::get<1>(t)), 4);
            f.write(reinterpret_cast<const char*>(&cons), 4);
            f.write(reinterpret_cast<const char*>(&tot), 4);
        }
        // -- M3 CLS tail -------------------------------------------------------
        // cortical lanes + eligibility traces. Persisted because the M2 lesson
        // applies directly: consolidated structure that is not written out is
        // structure that does not exist after a reload. Magic-guarded so a
        // pre-M3 file (which ends at the previous tail) still loads.
        // ALWAYS write the CLS header, even when empty. The loader reads tail
        // magics unconditionally and advances the stream, so an OMITTED CLS
        // block (an NTDA-only fabric has no cortical lanes) left the reader 4
        // bytes ahead and every later tail mis-matched -- NTDA silently failed
        // to load. Tails are a fixed sequence; skipping one breaks alignment.
        {
            const std::uint32_t cls_magic = 0x434C5331u;         // 'CLS1'
            f.write(reinterpret_cast<const char*>(&cls_magic), 4);
            using CRow = std::unordered_map<NodeId,
                        std::vector<std::pair<NodeId, float>>>::value_type;
            std::vector<const CRow*> cs;
            cs.reserve(cortical_.size());
            for (const auto& kv : cortical_) cs.push_back(&kv);
            std::sort(cs.begin(), cs.end(),
                      [](const auto* x, const auto* y) { return x->first < y->first; });
            std::uint32_t cn = static_cast<std::uint32_t>(cs.size());
            f.write(reinterpret_cast<const char*>(&cn), 4);
            for (const auto* kv : cs) {
                f.write(reinterpret_cast<const char*>(&kv->first), 4);
                std::uint32_t m = static_cast<std::uint32_t>(kv->second.size());
                f.write(reinterpret_cast<const char*>(&m), 4);
                for (const auto& l : kv->second) {
                    f.write(reinterpret_cast<const char*>(&l.first), 4);
                    f.write(reinterpret_cast<const char*>(&l.second), 4);
                }
            }
            std::vector<std::pair<std::uint64_t, float>> es(elig_.begin(), elig_.end());
            std::sort(es.begin(), es.end(),
                      [](const std::pair<std::uint64_t, float>& x,
                         const std::pair<std::uint64_t, float>& y) { return x.first < y.first; });
            std::uint32_t en = static_cast<std::uint32_t>(es.size());
            f.write(reinterpret_cast<const char*>(&en), 4);
            for (const auto& kv : es) {
                std::uint32_t a = static_cast<std::uint32_t>(kv.first >> 32);
                std::uint32_t b = static_cast<std::uint32_t>(kv.first & 0xffffffffull);
                f.write(reinterpret_cast<const char*>(&a), 4);
                f.write(reinterpret_cast<const char*>(&b), 4);
                f.write(reinterpret_cast<const char*>(&kv.second), 4);
            }
        }

        // -- M5 NTDA tail ------------------------------------------------------
        // volt_ MUST persist: a lane that had accumulated 1.5V toward a
        // threshold of 2 restarts from ZERO after a reload, so learning
        // behaviour after reload differs from behaviour before it. That is the
        // same defect class as M2's unserialised MetaLearner, reintroduced in
        // new code -- which is exactly why it needed catching again.
        // BOTH tails are written UNCONDITIONALLY, in a fixed order. Conditional
        // tails are what caused the CLS1 misalignment bug: the loader reads
        // magics unconditionally and advances the stream, so an omitted block
        // desynchronised every later one. Same lesson, third application.
        {
            const std::uint32_t ntd_magic = 0x4E544441u;        // 'NTDA'
            f.write(reinterpret_cast<const char*>(&ntd_magic), 4);
            const float th = cfg_.ntda_theta, ed = cfg_.edge_decay, dt = cfg_.edge_death;
            f.write(reinterpret_cast<const char*>(&th), 4);
            f.write(reinterpret_cast<const char*>(&ed), 4);
            f.write(reinterpret_cast<const char*>(&dt), 4);
            f.write(reinterpret_cast<const char*>(&ntda_spikes_), 8);
            f.write(reinterpret_cast<const char*>(&ntda_deaths_), 8);
            std::vector<std::pair<std::uint64_t, float>> vs(volt_.begin(), volt_.end());
            std::sort(vs.begin(), vs.end(),
                      [](const std::pair<std::uint64_t, float>& x,
                         const std::pair<std::uint64_t, float>& y) { return x.first < y.first; });
            std::uint32_t vn = static_cast<std::uint32_t>(vs.size());
            f.write(reinterpret_cast<const char*>(&vn), 4);
            for (const auto& kv : vs) {
                std::uint32_t a = static_cast<std::uint32_t>(kv.first >> 32);
                std::uint32_t b = static_cast<std::uint32_t>(kv.first & 0xffffffffull);
                f.write(reinterpret_cast<const char*>(&a), 4);
                f.write(reinterpret_cast<const char*>(&b), 4);
                f.write(reinterpret_cast<const char*>(&kv.second), 4);
            }
        }

        // -- L1 ENCODER tail ----------------------------------------------------
        // THIRD instance of the same defect (after M2's MetaLearner and M5's
        // voltage): runtime-only state that is not written out does not exist
        // after a reload. W learned by thousands of Oja updates is the whole
        // point of the encoder, and losing it silently makes --encoder-dim look
        // like it "did nothing" (stats reported encoder_dim: 0).
        {
            const std::uint32_t enc_magic = 0x454E4331u;        // 'ENC1'
            f.write(reinterpret_cast<const char*>(&enc_magic), 4);
            const bool have_W = !enc_W_.empty();
            const std::uint64_t z = 0;
            const std::uint64_t dd = have_W ? enc_dim_ : z;
            const std::uint64_t cc = have_W ? enc_chans_ : z;
            const std::uint64_t n  = have_W ? enc_W_.size() : z;
            f.write(reinterpret_cast<const char*>(&dd), 8);
            f.write(reinterpret_cast<const char*>(&cc), 8);
            f.write(reinterpret_cast<const char*>(&enc_eta_), 4);
            f.write(reinterpret_cast<const char*>(&enc_updates_), 8);
            f.write(reinterpret_cast<const char*>(&n), 8);
            if (have_W)
                f.write(reinterpret_cast<const char*>(enc_W_.data()),
                        static_cast<std::streamsize>(n * sizeof(float)));
        }

    }

    // v3.2.1: reports failure. A missing substrate used to load SILENTLY as an
    // empty fabric (every decide then deferred with honest_silence, and a
    // --router mis-path surfaced as the misleading "lacks router.json anchors").
    // Honest failure: false <=> the file could not be opened.
    bool load(const std::string& path) {
        std::ifstream f(path, std::ios::binary);
        if (!f) return false;
        nodes_.clear(); index_.clear(); out_.clear();
        csr_dirty_ = true;
        std::uint32_t n = 0; f.read(reinterpret_cast<char*>(&n), 4);
        nodes_.reserve(n);
        for (std::uint32_t i = 0; i < n; ++i) {
            std::uint32_t len = 0; f.read(reinterpret_cast<char*>(&len), 4);
            std::string c(len, '\0'); f.read(c.data(), len);
            float mass = 0; f.read(reinterpret_cast<char*>(&mass), 4);
            index_.emplace(c, static_cast<NodeId>(nodes_.size()));
            nodes_.push_back(Node{c, mass, 0.0f, 0.0f});
        }
        std::uint32_t lanes = 0; f.read(reinterpret_cast<char*>(&lanes), 4);
        for (std::uint32_t i = 0; i < lanes; ++i) {
            NodeId a = 0, b = 0; float w = 0;
            f.read(reinterpret_cast<char*>(&a), 4);
            f.read(reinterpret_cast<char*>(&b), 4);
            f.read(reinterpret_cast<char*>(&w), 4);
            out_[a].emplace_back(b, w);
        }
        // v2 provenance tail (optional): v1 files hit EOF here and stay clean.
        std::uint32_t dc = 0;
        if (f.read(reinterpret_cast<char*>(&dc), 4)) {
            for (std::uint32_t i = 0; i < dc; ++i) {
                std::uint32_t a = 0, b = 0, g = 0;
                if (!f.read(reinterpret_cast<char*>(&a), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&b), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&g), 4)) break;
                derived_gen_.emplace(lane_key(static_cast<NodeId>(a), static_cast<NodeId>(b)), g);
            }
        }
        // v3 evidence tail (optional): v2 files hit EOF here and stay clean.
        std::uint32_t ec = 0;
        if (f.read(reinterpret_cast<char*>(&ec), 4)) {
            for (std::uint32_t i = 0; i < ec; ++i) {
                std::uint32_t a = 0, b = 0, sup = 0, cnt = 0, clen = 0;
                std::uint64_t fs = 0, ls = 0;
                if (!f.read(reinterpret_cast<char*>(&a), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&b), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&sup), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&cnt), 4)) break;
                if (!f.read(reinterpret_cast<char*>(&fs), 8)) break;
                if (!f.read(reinterpret_cast<char*>(&ls), 8)) break;
                if (!f.read(reinterpret_cast<char*>(&clen), 4)) break;
                std::string ctx(clen, '\0');
                if (clen && !f.read(ctx.data(), clen)) break;
                LaneEvidence e;
                e.support_events = sup; e.counter_events = cnt;
                e.first_seq = fs; e.last_seq = ls; e.context = ctx;
                lane_evidence_[lane_key(static_cast<NodeId>(a), static_cast<NodeId>(b))] = e;
            }
        }
        // v3.2 semantic tail (optional, magic-guarded): pre-v3.2 files end at
        // the evidence tail; the failed magic read leaves semantics absent.
        std::uint32_t magic = 0;
        if (f.read(reinterpret_cast<char*>(&magic), 4) && magic == 0x53454D34u) {
            std::uint32_t has_sem = 0;
            if (f.read(reinterpret_cast<char*>(&has_sem), 4) && has_sem == 1) {
                const std::size_t n = nodes_.size();
                semvecs_.assign(n * SEM_DIMS, 0.0f);
                if (f.read(reinterpret_cast<char*>(semvecs_.data()),
                           static_cast<std::streamsize>(semvecs_.size() * sizeof(float)))) {
                    std::uint32_t total = 0;
                    if (f.read(reinterpret_cast<char*>(&total), 4)) {
                        sem_edges_.clear();
                        sem_edges_.reserve(total);
                        sem_off_.assign(n + 1, 0);
                        bool ok = true;
                        for (std::size_t i = 0; i <= n && ok; ++i) {
                            std::uint32_t off = 0;
                            ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&off), 4));
                            sem_off_[i] = off;
                        }
                        for (std::uint32_t k = 0; k < total && ok; ++k) {
                            std::uint32_t id = 0; float w = 0;
                            ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&id), 4));
                            if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&w), 4));
                            if (ok) sem_edges_.emplace_back(id, w);
                        }
                        if (ok) sem_edge_count_ = total;
                        else { sem_edges_.clear(); sem_off_.clear(); }
                    }
                }
                if (sem_edges_.empty()) semvecs_.clear();   // torn tail: stay off
            }
            std::uint32_t ctxc = 0;
            if (f.read(reinterpret_cast<char*>(&ctxc), 4)) {
                bool ok = true;
                for (std::uint32_t k = 0; k < ctxc && ok; ++k) {
                    std::uint64_t key = 0; std::uint32_t rq = 0, fb = 0;
                    ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&key), 8));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&rq), 4));
                    LaneCtx sig;
                    for (std::uint32_t j = 0; j < rq && ok; ++j) {
                        std::uint32_t w = 0;
                        ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&w), 4));
                        if (ok) sig.required.push_back(w);
                    }
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&fb), 4));
                    for (std::uint32_t j = 0; j < fb && ok; ++j) {
                        std::uint32_t w = 0;
                        ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&w), 4));
                        if (ok) sig.forbidden.push_back(w);
                    }
                    if (ok) lane_ctx_[key] = std::move(sig);
                }
            }
        }
        // v3.6.0 IDF tail (optional, magic-guarded 'IDF5'): pre-v3.6 files end
        // after the context block; the failed magic read leaves df_ empty and
        // every idf_factor at 1.0.
        std::uint32_t idf_magic = 0;
        if (f.read(reinterpret_cast<char*>(&idf_magic), 4) && idf_magic == 0x49444635u) {
            std::uint64_t lesson_count = 0;
            std::uint32_t dfc = 0;
            if (f.read(reinterpret_cast<char*>(&lesson_count), 8) &&
                f.read(reinterpret_cast<char*>(&dfc), 4)) {
                df_lessons_ = lesson_count;
                bool ok = true;
                for (std::uint32_t k = 0; k < dfc && ok; ++k) {
                    std::uint32_t id = 0, c = 0;
                    ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&id), 4));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&c), 4));
                    if (ok) df_[static_cast<NodeId>(id)] = c;
                }
                if (!ok) { df_.clear(); df_lessons_ = 0; }   // torn tail: stay inert
            }
        }
        // v3.7.0 DSTV tail (optional, magic-guarded 'DSTV'): pre-v3.7 files end
        // after the IDF block; the failed magic read leaves the distributional
        // field absent and the fabric replays exactly as before.
        std::uint32_t dstv_magic = 0;
        if (f.read(reinterpret_cast<char*>(&dstv_magic), 4) && dstv_magic == 0x44535456u) {
            std::uint32_t ddims = 0;
            if (f.read(reinterpret_cast<char*>(&ddims), 4) && ddims > 0 && ddims <= 4096) {
                const std::size_t need = nodes_.size() * static_cast<std::size_t>(ddims);
                distvecs_.assign(need, 0.0f);
                if (!f.read(reinterpret_cast<char*>(distvecs_.data()),
                            static_cast<std::streamsize>(need * sizeof(float)))) {
                    distvecs_.clear(); dist_dims_ = 0;       // torn tail: stay absent
                } else {
                    dist_dims_ = static_cast<int>(ddims);
                }
            }
        }
        // v3.9.0 ORD1 tail (optional, magic-guarded 'ORD1'): pre-v3.9 files end
        // after DSTV; the failed magic read leaves order_stats_ empty and every
        // order_factor at 1.0 — the replay contract.
        std::uint32_t ord_magic = 0;
        if (f.read(reinterpret_cast<char*>(&ord_magic), 4) && ord_magic == 0x4F524431u) {
            std::uint32_t orc = 0;
            if (f.read(reinterpret_cast<char*>(&orc), 4)) {
                bool ok = true;
                for (std::uint32_t k = 0; k < orc && ok; ++k) {
                    std::uint32_t a = 0, b = 0, cons = 0, tot = 0;
                    ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&a), 4));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&b), 4));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&cons), 4));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&tot), 4));
                    if (ok) order_stats_[lane_key(static_cast<NodeId>(a), static_cast<NodeId>(b))] =
                                {cons, tot};
                }
                if (!ok) order_stats_.clear();   // torn tail: stay inert
            }
        }
        // v3.9.1 M3 CLS tail (optional, magic-guarded 'CLS1'): pre-M3 files end
        // after ORD1, the failed magic read leaves cortical_ and elig_ empty,
        // and the substrate replays exactly as before.
        std::uint32_t cls_magic = 0;
        if (f.read(reinterpret_cast<char*>(&cls_magic), 4) && cls_magic == 0x434C5331u) {
            std::uint32_t cn = 0;
            if (f.read(reinterpret_cast<char*>(&cn), 4)) {
                bool ok = true;
                for (std::uint32_t k = 0; k < cn && ok; ++k) {
                    NodeId a = 0;
                    std::uint32_t m = 0;
                    ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&a), 4));
                    if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&m), 4));
                    if (ok) a = static_cast<NodeId>(a);
                    for (std::uint32_t j = 0; j < m && ok; ++j) {
                        NodeId b = 0;
                        float w = 0.0f;
                        ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&b), 4));
                        if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&w), 4));
                        if (ok) cortical_[a].emplace_back(static_cast<NodeId>(b), w);
                    }
                }
                if (!ok) cortical_.clear();       // torn tail: stay inert
                std::uint32_t en = 0;
                if (ok && f.read(reinterpret_cast<char*>(&en), 4)) {
                    for (std::uint32_t k = 0; k < en && ok; ++k) {
                        std::uint32_t a = 0, b = 0;
                        float e = 0.0f;
                        ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&a), 4));
                        if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&b), 4));
                        if (ok) ok = static_cast<bool>(f.read(reinterpret_cast<char*>(&e), 4));
                        if (ok) elig_[lane_key(static_cast<NodeId>(a), static_cast<NodeId>(b))] = e;
                    }
                    if (!ok) elig_.clear();
                }
                if (!cortical_.empty()) csr_dirty_ = true;
            }
        }
        // v3.9.1 M5 NTDA tail (optional, magic-guarded 'NTDA'): absent => ntda
        // stays off and the fabric replays exactly as pre-M5.
        std::uint32_t ntd_magic = 0;
        if (f.read(reinterpret_cast<char*>(&ntd_magic), 4) && ntd_magic == 0x4E544441u) {
            float th = 0, ed = 0, dt = 0;
            bool ok = f.read(reinterpret_cast<char*>(&th), 4)
                   && f.read(reinterpret_cast<char*>(&ed), 4)
                   && f.read(reinterpret_cast<char*>(&dt), 4);
            if (ok) ok = f.read(reinterpret_cast<char*>(&ntda_spikes_), 8)
                     && f.read(reinterpret_cast<char*>(&ntda_deaths_), 8);
            cfg_.ntda_on = true; cfg_.ntda_theta = th; cfg_.edge_decay = ed; cfg_.edge_death = dt;
            std::uint32_t vn = 0;
            if (ok && f.read(reinterpret_cast<char*>(&vn), 4)) {
                for (std::uint32_t k = 0; k < vn && ok; ++k) {
                    std::uint32_t a = 0, b = 0; float v = 0.0f;
                    ok = f.read(reinterpret_cast<char*>(&a), 4)
                      && f.read(reinterpret_cast<char*>(&b), 4)
                      && f.read(reinterpret_cast<char*>(&v), 4);
                    if (ok) volt_[lane_key(static_cast<NodeId>(a), static_cast<NodeId>(b))] = v;
                }
                if (!ok) volt_.clear();          // torn tail: start clean, not half
            }
        }
        std::uint32_t enc_magic_pre = 0;
        // ENC1 is written immediately after NTDA; both headers are always
        // present so the tail sequence can never desynchronise.
        std::uint32_t enc_magic = 0;
        if (f.read(reinterpret_cast<char*>(&enc_magic), 4) && enc_magic == 0x454E4331u) {
            std::uint64_t dd = 0, cc = 0, n = 0;
            if (f.read(reinterpret_cast<char*>(&dd), 8)
                && f.read(reinterpret_cast<char*>(&cc), 8)
                && f.read(reinterpret_cast<char*>(&enc_eta_), 4)
                && f.read(reinterpret_cast<char*>(&enc_updates_), 8)
                && f.read(reinterpret_cast<char*>(&n), 8)) {
                enc_dim_ = static_cast<std::size_t>(dd);
                enc_chans_ = static_cast<std::size_t>(cc);
                enc_W_.resize(static_cast<std::size_t>(n));
                // explicit bool(): the istream& overload is not contextually
                // convertible here.
                bool ok = static_cast<bool>(
                    f.read(reinterpret_cast<char*>(enc_W_.data()),
                           static_cast<std::streamsize>(n * sizeof(float))));
                if (!ok) enc_W_.clear();
                enc_y_.assign(enc_chans_, 0.0f);
            }
        }
        return true;
    }

    const SubstrateConfig& config() const { return cfg_; }

    // -- LANE PROVENANCE + SURGERY (derivation-layer surface, see derive.hpp) --
    // Derived lanes are laid by rule composition over existing fabric, not by
    // direct experience. They carry a generation counter: observed (Hebbian /
    // human-promoted) lanes are premise-grade (generation 0) and never lose
    // that grade; each composition step adds one generation. The generation
    // cap and use-time re-verification live in derive.hpp.
    void bind_derived(NodeId a, NodeId b, float w, std::uint32_t gen) {
        if (a == b || w <= 0.0f) return;
        set_lane_derived(a, b, w, gen);
        set_lane_derived(b, a, w, gen);
    }

    // Verifier-side weight set: assign exactly (down too), dissolve at zero.
    // Only meaningful for derived lanes; observed lanes are never touched by
    // callers (composition uses bind_derived, which is set-never-weaken).
    void set_derived_weight(NodeId a, NodeId b, float w) {
        auto it = out_.find(a);
        if (it == out_.end()) return;
        for (auto& l : it->second) {
            if (l.first == b) {
                if (w <= 0.0f) { dissolve_lane(a, b); return; }
                l.second = std::min(4.0f, w);
                csr_dirty_ = true;
                // keep symmetry: the mirror lane tracks the same evidence
                auto mit = out_.find(b);
                if (mit != out_.end())
                    for (auto& ml : mit->second)
                        if (ml.first == a) { ml.second = l.second; break; }
                return;
            }
        }
    }

    // Remove a lane (both directions) and its provenance.
    void dissolve_lane(NodeId a, NodeId b) {
        auto erase_one = [](std::vector<std::pair<NodeId, float>>& lanes, NodeId t) {
            for (std::size_t i = 0; i < lanes.size(); ++i)
                if (lanes[i].first == t) { lanes.erase(lanes.begin() + static_cast<std::ptrdiff_t>(i)); return; }
        };
        auto it = out_.find(a);
        if (it != out_.end()) erase_one(it->second, b);
        auto mit = out_.find(b);
        if (mit != out_.end()) erase_one(mit->second, a);
        derived_gen_.erase(lane_key(a, b));
        derived_gen_.erase(lane_key(b, a));
        drop_lane_ctx(a, b);                          // v3.2: signature dies with its lane
        drop_lane_ctx(b, a);
        csr_dirty_ = true;
    }

    std::uint32_t generation_of(NodeId a, NodeId b) const {
        auto it = derived_gen_.find(lane_key(a, b));
        return it == derived_gen_.end() ? 0u : it->second;
    }

    float lane_weight(NodeId a, NodeId b) const {
        auto it = out_.find(a);
        if (it == out_.end()) return 0.0f;
        for (const auto& l : it->second) if (l.first == b) return l.second;
        return 0.0f;
    }

    const std::string& concept_of(NodeId id) const { return nodes_[id].concept; }

    // Visitor over every lane (a -> b, weight). Directional view of the
    // symmetric fabric: each direction is visited once, like the save format.
    template <class F>
    void for_each_lane(F&& f) const {
        for (const auto& kv : out_)
            for (const auto& l : kv.second) f(kv.first, l.first, l.second);
    }

    // -- FABRIC SNAPSHOT (transactional derivation surface, see gate.hpp) ------
    // Exact, order-preserving capture of the lane fabric. settle() sums over
    // each node's lane vector in insertion order, so a faithful restore must
    // reproduce that order bit-for-bit: a (a,b,w) set is NOT enough, the
    // per-node vector sequences are the state. Used by the no-regression gate
    // to guarantee "revert" means "decide() output is byte-identical again".
    struct FabricSnapshot {
        std::vector<std::vector<std::pair<NodeId, float>>> adj;   // per node, exact order
        std::unordered_map<std::uint64_t, std::uint32_t> gen;     // lane_key -> generation
        std::size_t lanes = 0;
    };

    FabricSnapshot snapshot_fabric() const {
        FabricSnapshot snap;
        snap.adj.resize(nodes_.size());
        for (const auto& kv : out_) {
            if (kv.first >= snap.adj.size()) continue;
            snap.adj[kv.first] = kv.second;
            snap.lanes += kv.second.size();
        }
        snap.gen = derived_gen_;
        return snap;
    }

    void restore_fabric(const FabricSnapshot& snap) {
        out_.clear();
        for (std::size_t a = 0; a < snap.adj.size(); ++a)
            if (!snap.adj[a].empty()) out_[static_cast<NodeId>(a)] = snap.adj[a];
        derived_gen_ = snap.gen;
        csr_dirty_ = true;
    }

    void lanes_of(NodeId a, std::vector<std::pair<NodeId, float>>& out_lanes) const {
        auto it = out_.find(a);
        if (it == out_.end()) return;
        out_lanes = it->second;
    }

    // node-energy injection by id (dreamer probing; no decision fingerprint)
    void inject_energy(NodeId id, float energy) {
        if (id >= nodes_.size() || energy <= 0.0f) return;
        Node& n = nodes_[id];
        n.energy += energy / std::sqrt(n.mass);
        n.salience = 1.0f;   // TSDA: spike to ceiling on touch
    }

    // Runtime selection-mode overrides (CLI surface for the SI-faithful modes).
    // Physics constants stay fixed; only the two gating/window modes are
    // settable, and neither is persisted into substrate.bin — the saved model
    // is mode-neutral, so every binary can replay it under any mode.
    void set_source_modes(bool salience_gating_mode, bool miller_window_mode) {
        cfg_.salience_gating = salience_gating_mode;
        cfg_.miller_window   = miller_window_mode;
    }

private:
    struct Node { std::string concept; float mass; float energy; float salience; };

    static std::uint64_t lane_key(NodeId a, NodeId b) {
        return (static_cast<std::uint64_t>(a) << 32) | static_cast<std::uint32_t>(b);
    }

    // Composition-time lane write: raise-only (never weakens existing fabric),
    // and provenance is recorded once — a lane already premise-grade (laid by
    // experience) keeps generation 0; experience dominates derivation.
    void set_lane_derived(NodeId a, NodeId b, float w, std::uint32_t gen) {
        auto& lanes = out_[a];
        for (auto& l : lanes) {
            if (l.first == b) {
                if (w > l.second) { l.second = std::min(4.0f, w); csr_dirty_ = true; }
                derived_gen_.emplace(lane_key(a, b), gen);  // keep first (lowest) grade
                return;
            }
        }
        lanes.emplace_back(b, w);
        csr_dirty_ = true;
        derived_gen_.emplace(lane_key(a, b), gen);
        if (lanes.size() > cfg_.lane_cap) {              // evict weakest lane
            auto weakest = std::min_element(lanes.begin(), lanes.end(),
                [](const auto& x, const auto& y){ return x.second < y.second; });
            derived_gen_.erase(lane_key(a, weakest->first));
            drop_lane_ctx(a, weakest->first);            // v3.2: signature dies with its lane
            lanes.erase(weakest);
        }
    }

    void enforce_lane_decay() {
        for (auto& kv : out_)
            for (auto& l : kv.second) l.second *= cfg_.lane_decay;
        // M5 (b) EDGE DEATH. Conductivity decays faster than the weight does, and
        // a lane that falls below edge_death is DELETED outright rather than
        // lingering as a near-zero edge the readout must still consider. The
        // ontology therefore reshapes around what the fabric still uses, which
        // is a mechanism syfox did not previously have at all.
        if (cfg_.ntda_on) {
            for (auto& kv : out_)
                for (auto& l : kv.second) l.second *= cfg_.edge_decay;
            for (auto it = out_.begin(); it != out_.end();) {
                auto& vec = it->second;
                vec.erase(std::remove_if(vec.begin(), vec.end(),
                             [&](const std::pair<NodeId, float>& l) {
                                 if (l.second < cfg_.edge_death) { ++ntda_deaths_; return true; }
                                 return false;
                             }),
                          vec.end());
                if (vec.empty()) it = out_.erase(it); else ++it;
            }
        }
        csr_dirty_ = true;
    }

    SubstrateConfig cfg_;
    std::unordered_map<std::string, NodeId> index_;
    std::vector<Node> nodes_;
    std::unordered_map<NodeId, std::vector<std::pair<NodeId, float>>> out_;
    // M3: cortical lanes (consolidated structure) + the per-lane eligibility
    // traces that decide when a hippocampal lane is CONSOLIDATED. Eligibility is
    // monotone and never decayed: it counts distinct lesson co-occurrences,
    // which is the coincidence signal consolidation is supposed to key on.
    std::unordered_map<NodeId, std::vector<std::pair<NodeId, float>>> cortical_;
    std::unordered_map<std::uint64_t, float> elig_;
    // M5: per-lane accumulated voltage (the action potential) and the last
    // novelty/contradiction levels fed in by the caller.
    std::unordered_map<std::uint64_t, float> volt_;
    mutable float novelty_ = 0.0f;      // 0..1, set by decide-time OOD signal
    mutable float conflict_ = 0.0f;      // 0..1, set from the tension gate
    std::uint64_t ntda_spikes_ = 0;      // lanes that actually fired
    // -- L1 ENCODER PROJECTION (GHL/Oja) --------------------------------------
    // W : C x d, learned by LOCAL Hebbian plasticity -- Oja's subtractive term
    // plus a single global sign scalar g. No gradient, no backward pass, no
    // dependence on any other row. Rationale and citations in
    // ARCHITECTURE.md §15 (CLAPP, NeurIPS 2021; GHL).
    //
    // Channels are how an encoder ADDRESSES the lane fabric; they never replace
    // a lane with a vector. A 768 x 26,000 projection would be a second network
    // wearing a physics hat; 768 x 32 is not.
    std::size_t enc_chans_ = 0;          // C
    std::size_t enc_dim_ = 0;            // d (encoder embedding width)
    std::vector<float> enc_W_;           // C*d, row-major
    std::vector<float> enc_y_;           // C, last post-synaptic response
    float enc_eta_ = 0.01f;
    std::uint64_t enc_updates_ = 0;
    std::uint64_t ntda_deaths_ = 0;      // lanes deleted by edge death
    // Derivation provenance: lane key -> generation (sparse; observed lanes absent = 0)
    std::unordered_map<std::uint64_t, std::uint32_t> derived_gen_;
    // v3.6.0 document-frequency table (node -> lessons containing it) and the
    // lesson counter; persisted in the 'IDF5' tail. Empty table = IDF inert.
    std::unordered_map<NodeId, std::uint32_t> df_;
    std::uint64_t df_lessons_ = 0;
    // v3 Milestone 3: audit ledger (support/counter events, provenance window)
    std::unordered_map<std::uint64_t, LaneEvidence> lane_evidence_;
    // v3.7.0 distributional field (PPMI+SVD over the lane fabric); empty = absent
    std::vector<float> distvecs_;                      // nodes_ * dist_dims_ (row-major)
    int dist_dims_ = 0;
    // v3.2 semantic field state ------------------------------------------------
    std::vector<float> semvecs_;                       // nodes_ * SEM_DIMS (row-major); empty = absent
    std::vector<std::pair<NodeId, float>> sem_edges_;  // resonance edges, CSR payload
    std::vector<std::size_t> sem_off_;                 // CSR offsets, nodes_.size()+1
    std::size_t sem_edge_count_ = 0;
    std::unordered_map<std::uint64_t, LaneCtx> lane_ctx_;            // final signatures
    std::unordered_map<std::uint64_t, std::unordered_map<NodeId, float>> lane_ctx_acc_; // learn-time scratch
    std::unordered_set<NodeId> present_;               // tokens injected since last reset_field()
    bool sem_enabled_ = true;                          // runtime kill switch (--no-semantics)
    int  hop_depth_ = 1;                               // v3.4 readout walk depth (1 = legacy)
    std::uint64_t state_hash_ = 0;   // decision fingerprint for Miller sampling
    float last_cap_ = 0.0f;          // cap used by the last settle()
    SettleTrace last_trace_;         // v3.8.0 TASK 4: settle self-verification trace
    // v3.8.0 TASK 2: spectral channel gate (empty = inactive)
    std::vector<float> channel_coords_;
    float channel_sigma_ = 0.0f;
    // v3.9.0: MCPE/CPME peak-diffusion-gain + VGHT readout-gate bookkeeping
    bool momentum_on_ = false;
    std::vector<float> max_gain_;    // per-node peak pass-over-pass gain
    bool vght_on_ = false;
    mutable float vght_gated_sum_ = 0.0f;   // lanes with C >= 0.5 (const readout)
    mutable std::size_t vght_walked_ = 0;   // lanes traversed by gated walks
    // v3.9.0 (CPME-II discrete surrogate) per-lane order evidence:
    // key = lane_key(a,b); value = {consistent_exposures, total_exposures}.
    // The engine notes, per (state token, anchor) lane, whether the token
    // PRECEDED the anchor in the lesson text. hebbian_lesson scales the
    // state-side bind by Phi = (1-alpha) + alpha*(2c/t - 1) with alpha=0.6:
    // order-consistent lanes bind at full strength (Phi=1), reversed lanes
    // at 0.4-0.6 = -0.2 -> floored at 0, mixed lanes in between. Persisted
    // in the 'ORD1' tail. Empty map = factor 1.0 everywhere (replay no-op).
    std::unordered_map<std::uint64_t, std::pair<std::uint32_t, std::uint32_t>> order_stats_;
    bool order_on_ = false;

    // -- M5: CSR MIRROR of the out-lane fabric (settle hot path) -------------
    // Flattened contiguous buffers (offsets + targets + weights) built lazily
    // from out_ and invalidated by every lane/node mutation above. The build
    // copies each source's lane vector VERBATIM, so per-source iteration order
    // is exactly the map order settle() used before: every float sum keeps its
    // exact operand order and results stay bit-identical. This is the
    // cache-friendly layout the GPU design generalizes (ARCHITECTURE.md S15).
    // mutable: a lazy cache — logically const rebuild inside const settle paths.
    mutable std::vector<std::size_t> csr_off_;   // nodes_.size()+1
    mutable std::vector<NodeId>      csr_dst_;
    mutable std::vector<float>       csr_w_;
    mutable bool csr_dirty_ = true;
    bool parallel_settle_ = true;        // OMP builds honor this; sequential builds ignore it

private:
    // -- v3.2 semantic helpers --------------------------------------------------
    bool sem_active() const {
        return sem_enabled_ && !semvecs_.empty() && !sem_edges_.empty();
    }

    // fixed deterministic direction for omega_semantic (Stage 1); computed once
    static const float* omega_dir() {
        static const std::vector<float> dir = [] {
            std::vector<float> d(SEM_DIMS, 0.0f);
            for (std::size_t k = 0; k < SEM_DIMS; ++k) {
                std::uint64_t h = 1469598103934665603ull;
                const std::string seed = "syfox-omega-axis-" + std::to_string(k);
                for (unsigned char c : seed) { h ^= c; h *= 1099511628211ull; }
                d[k] = (static_cast<float>(h & 0xFFFFu) / 32767.5f) - 1.0f;  // [-1,1)
            }
            float n2 = 0.0f;
            for (float v : d) n2 += v * v;
            if (n2 > 0.0f) { const float inv = 1.0f / std::sqrt(n2); for (auto& v : d) v *= inv; }
            return d;
        }();
        return dir.data();
    }

    // coupled-oscillator detuning: equal frequencies exchange fully, far ones don't
    float freq_match(NodeId a, NodeId b) const {
        return 1.0f - 0.5f * std::fabs(semantic_freq(a) - semantic_freq(b));
    }

    // Stage 2 settle-side factor: how much of a lane's weight flows this pass
    float lane_flow_factor(NodeId a, NodeId b) const {
        // v3.8.0 TASK 2: the return value is now f_ctx * f_channel. When the
        // channel gate is inactive the second factor is skipped entirely and
        // every branch below returns the exact v3.7 value (bit-identity).
        float f = 1.0f;
        auto it = lane_ctx_.find(lane_key(a, b));
        if (it != lane_ctx_.end()) {
            const LaneCtx& sig = it->second;
            bool forbidden_hit = false;
            for (const NodeId w : sig.forbidden)
                if (present_.count(w)) { forbidden_hit = true; break; }
            if (forbidden_hit) {
                f = cfg_.ctx_forbidden;
            } else if (!sig.required.empty()) {
                std::size_t matched = 0;
                for (const NodeId w : sig.required)
                    if (present_.count(w)) ++matched;
                if (matched < sig.required.size())
                    f = cfg_.ctx_missing
                      + (1.0f - cfg_.ctx_missing)
                      * (static_cast<float>(matched) / static_cast<float>(sig.required.size()));
            }
        }
        if (!channel_coords_.empty()) {
            const float d = channel_coords_[a] - channel_coords_[b];
            f *= std::exp(-(d * d) / (2.0f * channel_sigma_ * channel_sigma_));
        }
        return f;
    }

    void drop_lane_ctx(NodeId a, NodeId b) {
        lane_ctx_.erase(lane_key(a, b));
        lane_ctx_acc_.erase(lane_key(a, b));
    }

    void build_csr() const {
        const std::size_t n = nodes_.size();
        // M3: cortical lanes are an ADDITIONAL conductance population. They go
        // into the same CSR so they propagate exactly like hippocampal lanes —
        // the difference is not how they are read but that enforce_lane_decay
        // never touches them. Order is hippocampal-then-cortical per node, which
        // keeps the float accumulation order deterministic.
        auto row_len = [&](std::size_t i) -> std::size_t {
            std::size_t k = 0;
            auto it = out_.find(static_cast<NodeId>(i));
            if (it != out_.end()) k += it->second.size();
            auto ct = cortical_.find(static_cast<NodeId>(i));
            if (ct != cortical_.end()) k += ct->second.size();
            return k;
        };
        csr_off_.assign(n + 1, 0);
        for (std::size_t i = 0; i < n; ++i) csr_off_[i + 1] = csr_off_[i] + row_len(i);
        csr_dst_.resize(csr_off_[n]);
        csr_w_.resize(csr_off_[n]);
        for (std::size_t i = 0; i < n; ++i) {
            std::size_t k = csr_off_[i];
            auto it = out_.find(static_cast<NodeId>(i));
            if (it != out_.end())
                for (const auto& lane : it->second) {
                    csr_dst_[k] = lane.first;
                    csr_w_[k]   = lane.second;
                    ++k;
                }
            auto ct = cortical_.find(static_cast<NodeId>(i));
            if (ct != cortical_.end())
                for (const auto& lane : ct->second) {
                    // cortical_gain scales the cortical population down so it
                    // supports rather than dominates the hippocampal trace.
                    csr_dst_[k] = lane.first;
                    csr_w_[k]   = lane.second * (cfg_.cls_split ? cfg_.cortical_gain : 1.0f);
                    ++k;
                }
        }
        csr_dirty_ = false;
    }

public:
    // parallel settle toggle (M5): on OMP builds, >1 threads use per-thread
    // scatter buffers combined in fixed thread order — bit-identical to the
    // sequential path (test-verified). Off builds always run sequential.
    void set_parallel_settle(bool on) { parallel_settle_ = on; }
    bool parallel_settle_enabled() const {
#if defined(_OPENMP)
        return parallel_settle_;
#else
        return false;
#endif
    }
    // fabric density for the GPU-gate report (M5): edges / possible directed
    // edges, plus mean out-degree. Pure measurements, no claims.
    double fabric_density() const {
        const double n = static_cast<double>(nodes_.size());
        if (n < 2.0) return 0.0;
        return static_cast<double>(lane_count()) / (n * (n - 1.0));
    }
    double mean_out_degree() const {
        const double n = static_cast<double>(nodes_.size());
        return n > 0.0 ? static_cast<double>(lane_count()) / n : 0.0;
    }
};

// ---------------------------------------------------------------------------
// v3.2 build_semantics — defined outside the class for OMP clarity.
// Deterministic end to end: lexical hashing is a pure function of the string;
// the grounding passes iterate node ids ascending; the neighbour search
// parallelizes over ROWS (each row's top-k depends only on that row), with
// ties broken by (cosine desc, neighbour id asc).
// ---------------------------------------------------------------------------
inline void Substrate::build_semantics(int threads) {
    const std::size_t n = node_count();
    semvecs_.clear(); sem_edges_.clear(); sem_off_.clear(); sem_edge_count_ = 0;
    if (n == 0) return;
    // Cost guard: the neighbour search is O(n^2 * dims). Dedicated fabrics are
    // 10^3..10^4 nodes (seconds); a 43k-node fabric is minutes once. Beyond
    // 65k the build is refused — semantics stays absent, the model stays
    // fully usable (honest scope, not silent failure).
    if (n > 65000) return;
    (void)threads;
#if defined(_OPENMP)
    if (threads > 0) omp_set_num_threads(threads);
#endif

    // Step 1 — LEXICAL LAYER: signed character-trigram hashing, L2-normalized.
    semvecs_.assign(n * SEM_DIMS, 0.0f);
    for (std::size_t i = 0; i < n; ++i) {
        float* v = &semvecs_[i * SEM_DIMS];
        const std::string& s = concept_of(static_cast<NodeId>(i));
        auto bump = [&](const std::string& g, float w) {
            std::uint64_t h = 1469598103934665603ull;
            for (unsigned char c : g) { h ^= c; h *= 1099511628211ull; }
            const std::size_t bucket = static_cast<std::size_t>(h % SEM_DIMS);
            const float sign = ((h >> 40) & 1u) ? 1.0f : -1.0f;
            v[bucket] += sign * w;
        };
        if (s.size() < 3) {
            bump(s, 2.0f);
        } else {
            for (std::size_t p = 0; p + 3 <= s.size(); ++p) bump(s.substr(p, 3), 1.0f);
            bump(s, 1.5f);                       // whole-token identity, too
        }
        float n2 = 0.0f;
        for (std::size_t k = 0; k < SEM_DIMS; ++k) n2 += v[k] * v[k];
        if (n2 > 0.0f) {
            const float inv = 1.0f / std::sqrt(n2);
            for (std::size_t k = 0; k < SEM_DIMS; ++k) v[k] *= inv;
        }
    }

    // Step 2 — FABRIC GROUNDING: two smoothing passes over the Hebbian lanes.
    // v_i <- normalize(v_i + beta * lane-weighted mean of v_j). This is how
    // the fabric's own co-occurrence structure shapes the vectors — meaning
    // diffusing along the lanes, not a fit.
    for (int pass = 0; pass < 2; ++pass) {
        std::vector<float> nx(n * SEM_DIMS, 0.0f);
        for (std::size_t i = 0; i < n; ++i) {
            const float* v = &semvecs_[i * SEM_DIMS];
            float* o = &nx[i * SEM_DIMS];
            for (std::size_t k = 0; k < SEM_DIMS; ++k) o[k] = v[k];
            auto it = out_.find(static_cast<NodeId>(i));
            if (it != out_.end() && !it->second.empty()) {
                float wsum = 0.0f;
                float avg[SEM_DIMS];
                for (std::size_t k = 0; k < SEM_DIMS; ++k) avg[k] = 0.0f;
                for (const auto& l : it->second) {
                    const float* u = &semvecs_[l.first * SEM_DIMS];
                    for (std::size_t k = 0; k < SEM_DIMS; ++k) avg[k] += l.second * u[k];
                    wsum += l.second;
                }
                if (wsum > 0.0f)
                    for (std::size_t k = 0; k < SEM_DIMS; ++k)
                        o[k] = v[k] + cfg_.sem_beta * (avg[k] / wsum);
            }
            float n2 = 0.0f;
            for (std::size_t k = 0; k < SEM_DIMS; ++k) n2 += o[k] * o[k];
            if (n2 > 0.0f) {
                const float inv = 1.0f / std::sqrt(n2);
                for (std::size_t k = 0; k < SEM_DIMS; ++k) o[k] *= inv;
            }
        }
        semvecs_.swap(nx);
    }
    (void)omega_dir();   // materialize the fixed omega axis once

    // Step 3 — RESONANCE EDGES: top-k cosine neighbours with cos >= theta.
    const int kmax = std::max(1, cfg_.sem_neighbors);
    std::vector<std::vector<std::pair<float, NodeId>>> per(n);
#if defined(_OPENMP)
#pragma omp parallel for schedule(static)
#endif
    for (std::ptrdiff_t ii = 0; ii < static_cast<std::ptrdiff_t>(n); ++ii) {
        const std::size_t i = static_cast<std::size_t>(ii);
        const float* vi = &semvecs_[i * SEM_DIMS];
        float vi2 = 0.0f;
        for (std::size_t d = 0; d < SEM_DIMS; ++d) vi2 += vi[d] * vi[d];
        if (vi2 <= 0.0f) continue;                       // dark node: no resonance
        std::vector<std::pair<float, NodeId>> best;      // kept sorted: cos desc, id asc
        best.reserve(static_cast<std::size_t>(kmax) + 1);
        auto better = [](const std::pair<float, NodeId>& x, const std::pair<float, NodeId>& y) {
            return x.first != y.first ? x.first > y.first : x.second < y.second;
        };
        for (std::size_t j = 0; j < n; ++j) {
            if (j == i) continue;
            const float* vj = &semvecs_[j * SEM_DIMS];
            float dot = 0.0f;
            for (std::size_t d = 0; d < SEM_DIMS; ++d) dot += vi[d] * vj[d];
            if (dot < cfg_.sem_theta) continue;
            if (best.size() < static_cast<std::size_t>(kmax)) {
                best.emplace_back(dot, static_cast<NodeId>(j));
                std::sort(best.begin(), best.end(), better);
            } else if (better({dot, static_cast<NodeId>(j)}, best.back())) {
                best.back() = {dot, static_cast<NodeId>(j)};
                std::sort(best.begin(), best.end(), better);
            }
        }
        per[i] = std::move(best);
    }
    sem_off_.assign(n + 1, 0);
    std::size_t running = 0;
    for (std::size_t i = 0; i < n; ++i) {
        sem_off_[i] = running;
        running += per[i].size();
    }
    sem_off_[n] = running;
    sem_edges_.reserve(running);
    for (std::size_t i = 0; i < n; ++i)
        for (const auto& p : per[i]) sem_edges_.emplace_back(p.second, p.first);
    sem_edge_count_ = running;
}

// ---------------------------------------------------------------------------
// Tokenizer — plain lexical split (I/O concern, not pattern matching):
// lowercase, split on non-alphanumeric, keep tokens of length 2..24.
// ---------------------------------------------------------------------------
// Light deterministic lexical folding (string normalization, not ML):
//   plurals:    "refunds"->"refund", "invoices"->"invoice"  (len>=4, not "ss")
//   past/verb:  "charged"->"charg", "charge"->"charg"      (len>=5)
//   gerund:     "billing"->"bill", "running"->"runn"       (len>=6)
// Consistency matters, not linguistic truth: the same fold applies to state,
// options, and instructions, so matched forms always meet.
inline std::string fold(const std::string& w) {
    std::string s = w;
    if (s.size() >= 4 && s.compare(s.size() - 2, 2, "ss") != 0 && s.back() == 's') s.pop_back();
    if (s.size() >= 6 && s.compare(s.size() - 3, 3, "ing") == 0) { s.erase(s.size() - 3); return s; }
    if (s.size() >= 5 && s.compare(s.size() - 2, 2, "ed") == 0) s.erase(s.size() - 2);
    if (s.size() >= 5 && s.back() == 'e') s.pop_back();
    return s;
}

inline std::vector<std::string> tokenize(const std::string& text) {
    std::vector<std::string> out;
    std::string cur;
    auto emit = [&]() {
        if (cur.size() >= 2 && cur.size() <= 24) out.push_back(fold(cur));
        cur.clear();
    };
    for (char raw : text) {
        unsigned char c = static_cast<unsigned char>(raw);
        if (std::isalnum(c)) cur.push_back(static_cast<char>(std::tolower(c)));
        else if (!cur.empty()) emit();
    }
    if (!cur.empty()) emit();
    return out;
}

} // namespace si