perf(evolution): ⚡ halve the follower-marking epoch stamp

AGENTS.md

@@ -104,6 +104,24 @@ Key files:
   for the mutating/collecting paths, which run on the partitions' own pinned masters; `sum_partitions_`,
   `fold_partitions_`, `first_partition_` for reads off quiescent partitions) rather than hand-rolling a
   `run_on_all` loop — the declarations record which helper is legal where.
+- **The operator memory ledger**: `detail::MPOperatorMemoryBreakdown` in
+  `cpp/monoprop/detail/operator/MPOperator.h`, surfaced to Python as `operator_memory_breakdown()` by
+  `src/monoprop/bindings/binder.h`. Plain fields are summed by `total_bytes()`; `d_`-prefixed ones are
+  diagnostics *of* those fields and are excluded from the sum so they can never double-count. Two rules
+  that have already cost measurements:
+    - **Adding a field to `total_bytes()` makes `total_bytes` step up, which is not a regression.**
+      `matched_scratch_bytes` — the `detail::MatchedEpochSet` stamp array, one entry per term — was
+      counted by no field at all until it was added, so `total_bytes` rises by
+      `sizeof(MatchedEpochSet::Stamp)` per term across that change while the resident memory it names
+      *fell*. Any A/B or benchmark series keyed on `total_bytes` reads it the wrong way round; subtract
+      the new field when comparing across such a change, and say so in the PR body.
+    - **`init_operator_bytes` tracks the OBSERVABLE, not the operator.** It measures `bucket_count()` of
+      `init_op_map`, which `get_operator()` binds and then *releases* by swapping with an empty map —
+      `erase` and `clear` leave `bucket_count()` untouched, so a drained map otherwise keeps its whole
+      bucket array for the life of the run. With a single-site observable the map never grows and the
+      release is worth a few hundred bytes; with a multi-million-term observable it has been worth
+      gigabytes. Never quote one workload's figure as the other's, and read `d_init_operator_entries`
+      beside it — bytes with no entries means retained buckets, not retained terms.
 
 
 ### Environment Management

cpp/include/monoprop/MonomialPropagator.h

@@ -132,7 +132,11 @@ class MonomialPropagator {
         if (partition_group_) {
             return partitioned_operator_memory_usage_();
         }
-        return detail::estimate_memory_usage(mp_op_);
+        // matched_scratch_ is propagator-owned, so estimate_memory_usage() cannot see it; it is filled in
+        // here so that the facade's per-partition sum picks up one stamp array per partition.
+        auto breakdown = detail::estimate_memory_usage(mp_op_);
+        breakdown.matched_scratch_bytes = matched_scratch_.memory_bytes();
+        return breakdown;
     }
 
     auto graph_layers() const -> size_t { return partition_group_ ? partitioned_graph_layers_() : graph_.layers(); }

cpp/monoprop/detail/evolution/layer_build/Common.h

@@ -29,22 +29,59 @@ namespace monoprop::detail {
 // Marks matched followers without a per-gate O(n) memset: one counter bump clears every mark. Reused
 // across gates.
 struct MatchedEpochSet {
-    std::vector<uint32_t> epoch_;
-    uint32_t cur_ = 0;
-
-    // u32 wrap resets the array — once per 2^32-1 gates.
+    // The stamp is 16-bit, not 32-bit. There is one entry per term, so at 4 bytes plus geometric
+    // overshoot the array measured 6-7.6 B/term — of the same order as the whole inverted index — for a
+    // counter that never leaves this struct. Halving it is worth a measured -2.62 B/term of peak RSS at
+    // 1x16 and -2.53 B/term at 8x16, 6/6 with non-overlapping distributions, at no cost in time. Width
+    // is the only thing that changes: a stamp is never serialised, never exchanged between ranks, and
+    // never compared against anything but cur_.
+    //
+    // The shrink_to_fit() below is NOT part of that saving and is NOT called — see the note in
+    // MonomialPropagator<NumModes>::initialize_operator_caches_() for the measurement that rejected it.
+    // So this array carries geometric-growth overshoot on top of the 2 B/term, deliberately.
+    using Stamp = uint16_t;
+
+    std::vector<Stamp> epoch_;
+    Stamp cur_ = 0;
+
+    // Wrap resets the array — once per 65535 gates, where the u32 stamp took 2^32.
+    //
+    // That makes this branch LIVE rather than effectively dead: any run of more than 65535 gate
+    // applications reaches it, and a long circuit reaches it repeatedly. The std::fill is what keeps it
+    // correct — without it a slot last stamped at epoch e aliases onto the next epoch e and reads as
+    // marked. `matched_epoch_stamp_wrap_reached_by_gate_count` in cpp/tests/evolution_detail_tests.cpp
+    // is the only test that pins the fill (it cycles a whole period rather than assigning to cur_); it
+    // and its companion `matched_epoch_stamp_wrap_resets` are load-bearing — do not delete either in a
+    // refactor.
     auto begin_gate(size_t n) -> void {
-        if (cur_ == std::numeric_limits<uint32_t>::max()) {
-            std::fill(epoch_.begin(), epoch_.end(), 0);
+        if (cur_ == std::numeric_limits<Stamp>::max()) {
+            std::fill(epoch_.begin(), epoch_.end(), Stamp{0});
             cur_ = 0;
         }
         ++cur_;
         if (epoch_.size() < n) {
-            epoch_.resize(n, 0);
+            epoch_.resize(n, Stamp{0});
         }
     }
     auto mark(size_t i) -> void { epoch_[i] = cur_; }
     [[nodiscard]] auto is_marked(size_t i) const -> bool { return epoch_[i] == cur_; }
+    // One stamp per term, so this tracks the operator: it belongs in the memory breakdown even though it
+    // is propagator scratch rather than operator state. Reported as matched_scratch_bytes, which only
+    // MonomialPropagator can fill in because only it owns this object.
+    [[nodiscard]] auto memory_bytes() const -> size_t { return epoch_.capacity() * sizeof(Stamp); }
+    // Release geometric-growth overshoot at a quiescent point. Gated on a bounded margin (1/8 = 12.5%
+    // dead capacity) rather than unconditional: an exact fit taken while the operator is still growing
+    // reallocs the whole array on the very next gate that appends a term.
+    //
+    // CURRENTLY UNCALLED, on purpose. The only quiescent point that exists is
+    // initialize_operator_caches_(), which runs after every propagate() rather than once at the end of a
+    // run, and calling it there measured +0.2821 GiB of peak RSS. Kept because it is the correct shape
+    // for a caller that genuinely is at the end of a run; do not wire it back into the gate loop.
+    auto shrink_to_fit() -> void {
+        if (epoch_.capacity() > epoch_.size() + (epoch_.size() / 8)) {
+            epoch_.shrink_to_fit();
+        }
+    }
 };
 
 // A trivial aggregate on purpose — not std::pair — so DefaultInitVector can skip the zero-fill and lower

cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl

@@ -509,6 +509,25 @@ auto MonomialPropagator<NumModes>::initialize_operator_caches_() -> void {
     (void)mp_op_.inverted_index();
     mp_op_.op_coeffs.shrink_to_fit();
     mp_op_.shrink_state_to_fit();
+    // The follower-marking scratch is deliberately NOT shrunk here. An earlier revision of this change
+    // called matched_scratch_.shrink_to_fit() on this line, on the premise that the array "grows only
+    // when the term count does", which is false on Hubbard: this runs at the end of every gate loop --
+    // after EVERY build_graph() and propagate(), not once at the end of a run -- and Hubbard's 29 Trotter
+    // steps reach it 29 times with a term count that grows at each step, so the exact fit is voided and
+    // re-paid on a ~194 MB buffer. Its A/B read +0.2821 GiB of peak RSS (+3.12 B/term), 6/6 with
+    // non-overlapping distributions, at the 1x16 hubbard cell. The shrink also cannot free resident
+    // memory: resize(n, 0) never writes past n, so the released capacity was never faulted in.
+    //
+    // The operator's ROW STORE is left out for the same reason, and its slack is the larger of the two
+    // (measured 4.1-6.0 B/term, reported as d_terms_slack_bytes). Tried and measured on the branch that
+    // proposed it: `propagate` 1.030x slower on Hubbard c10, 6/6, p=0.031, against nothing resolvable on
+    // a single-propagate workload. Geometric growth plus shrink-to-fit churns whenever both run
+    // repeatedly, so collecting that slack needs a row store whose growth does not copy, not a
+    // better-timed shrink.
+    //
+    // Consequence for the memory breakdown: matched_scratch_bytes is capacity * sizeof(Stamp), so
+    // without the shrink it carries geometric-growth overshoot and no longer lands on exactly
+    // sizeof(Stamp) per term. The u16 width is untouched.
 }
 
 template <size_t NumModes>

cpp/monoprop/detail/operator/MPOperator.h

@@ -26,7 +26,7 @@
 #include <format>
 #include <print>
 
 #include "monoprop/TypeAliases.h"
 #include "monoprop/Utilities.h"
 #include "monoprop/detail/operator/InvertedIndex.h"
 #include "monoprop/detail/operator/OperatorIndex.h"
@@ -63,18 +63,18 @@
     // The store is non-copyable/non-movable, so it is heap-owned by unique_ptr (keeping MPOperator
     // itself cheaply movable). Always non-null.
     std::unique_ptr<OperatorIndex<NumModes>> store = std::make_unique<OperatorIndex<NumModes>>();
     VecD op_coeffs = {};
     // Only fully-paired terms score nonzero (see score_new_state_rows_), which on production models is
     // ~0.07% of the rows -- a dense vector here is 99.9% zeros. state_rows_ is strictly ascending: rows are
     // scored in ascending order and the set is only ever appended to.
     std::vector<TermIndex> state_rows_ = {};
     VecD state_vals_ = {};         // parallel to state_rows_; every entry is a unit phase (+-1), never 0
     size_t state_scored_rows_ = 0; // rows [0, state_scored_rows_) have been scored into state_rows_/state_vals_
     // The dense state: empty in Heisenberg unless a caller asks dense_state() to cache one; in Schrödinger
     // it is the live coefficient vector evolution mutates in place.
     VecD state_coeffs = {};
     MonomialMap<NumModes> init_op_map = {};
     VecZ initial_state = {};
     // Set once at propagator construction.
     Basis basis = Basis::Majorana;
     mutable std::optional<InvertedIndex<NumModes>> inverted_index_ = std::nullopt;
@@ -116,8 +116,18 @@
         return *inverted_index_;
     }
 
-    // Pending init_op_map terms are erased after the lookup loop: the flat_map is not iterable while
-    // mutating.
+    // Binds pending init_op_map terms to their rows and drops them from the map. Both loops below only
+    // READ the map, so the "not iterable while mutating" constraint that forced the old
+    // collect-then-erase still holds; what is gone is the vector of Monomial keys it collected into.
+    //
+    // The map is then RELEASED, not merely emptied: `erase` and `clear` both leave bucket_count(), and
+    // bucket_count() is what estimate_memory_usage() reports as init_operator_bytes, so a fully drained
+    // map used to hold its whole bucket array for the life of the run. Swap with an empty map.
+    //
+    // What that is worth tracks the OBSERVABLE, not the operator, and the two must be quoted separately:
+    // it is 1,189 bytes in total across the Hubbard and Pauli benchmark anchors, whose observables are
+    // single-site so the map never grows past its initial buckets, and 39.58 B/term — about 1.15 GB —
+    // on a workload with a 7M-term observable.
     auto get_operator() -> const VecD & {
         if (size() == op_coeffs.size()) {
             return op_coeffs;
@@ -129,18 +139,29 @@
             return op_coeffs;
         }
 
-        std::vector<Monomial<NumModes>> del;
+        size_t bound = 0;
         for (const auto &kv : init_op_map) {
-            const auto &mono = kv.first;
-            const auto coeff = kv.second;
-            if (const auto found = store->find(mono)) {
-                op_coeffs[*found] = coeff;
-                del.push_back(mono);
+            if (const auto found = store->find(kv.first)) {
+                op_coeffs[*found] = kv.second;
+                ++bound;
             }
         }
 
-        for (const auto &mono : del) {
-            init_op_map.erase(mono);
+        if (bound == init_op_map.size()) {
+            // Everything bound: drop the map outright, buckets included.
+            MonomialMap<NumModes>{}.swap(init_op_map);
+        }
+        else if (bound != 0) {
+            // Partial: rebuild from what is still pending, which sizes the buckets to it. The second
+            // find() pass is affordable because this path only runs while terms remain unbound.
+            MonomialMap<NumModes> pending;
+            pending.reserve(init_op_map.size() - bound);
+            for (const auto &kv : init_op_map) {
+                if (!store->find(kv.first)) {
+                    pending.emplace(kv.first, kv.second);
+                }
+            }
+            init_op_map.swap(pending);
         }
 
         return op_coeffs;
@@ -238,7 +259,7 @@
         }
 
         VecZ new_inds(size() - state_scored_rows_);
         std::iota(new_inds.begin(), new_inds.end(), state_scored_rows_);
 
         const auto paired_inds = is_fully_paired<NumModes>(new_inds, *store);
         state_rows_.reserve(state_rows_.size() + paired_inds.size());
@@ -279,7 +300,7 @@
 
 template <typename FlatMap>
 inline auto unordered_flat_map_storage_bytes(const FlatMap &map) -> size_t {
     return sizeof(FlatMap) + map.bucket_count() * (sizeof(typename FlatMap::value_type) + sizeof(unsigned char));
 }
 
 template <size_t NumModes>
@@ -291,6 +312,10 @@
     size_t init_operator_bytes = 0;
     size_t initial_state_bytes = 0;
     size_t inverted_index_bytes = 0;
+    // One epoch stamp per term (detail::MatchedEpochSet). It is propagator-owned rather than operator
+    // state, so only MonomialPropagator::operator_memory_usage() can fill it in and it stays 0 on a bare
+    // estimate_memory_usage() of an operator. It is real resident memory and IS part of total_bytes().
+    size_t matched_scratch_bytes = 0;
 
     // Diagnostics: breakdowns of the fields above, deliberately excluded from total_bytes() so they can
     // never double-count.
@@ -300,10 +325,19 @@
     size_t operator_terms_slack_bytes = 0; // of operator_terms_bytes: unused geometric-growth capacity
     // of state_coeffs_bytes: entries of the state that are not exactly 0.0
     size_t state_coeffs_nonzero = 0;
-
+    // Live entries behind init_operator_bytes. That field measures bucket_count(), not size(), and
+    // get_operator() drains the map — so a large byte count beside a zero entry count means retained
+    // buckets, not retained terms.
+    size_t init_operator_entries = 0;
+
+    // matched_scratch_bytes is part of this sum. A build without that field reports a total lower by
+    // roughly sizeof(MatchedEpochSet::Stamp) per term while holding the SAME or more resident memory —
+    // the array existed, it was simply counted nowhere. So an A/B or a benchmark series keyed on
+    // total_bytes across such a boundary reads a reporting improvement as a memory regression. Subtract
+    // matched_scratch_bytes before comparing, or re-measure the baseline.
     auto total_bytes() const -> size_t {
         return operator_terms_bytes + op_coeffs_bytes + state_coeffs_bytes + indexing_bytes + init_operator_bytes
-               + initial_state_bytes + inverted_index_bytes;
+               + initial_state_bytes + inverted_index_bytes + matched_scratch_bytes;
     }
 
     auto operator+=(const MPOperatorMemoryBreakdown &o) -> MPOperatorMemoryBreakdown & {
@@ -314,6 +348,8 @@
         init_operator_bytes += o.init_operator_bytes;
         initial_state_bytes += o.initial_state_bytes;
         inverted_index_bytes += o.inverted_index_bytes;
+        matched_scratch_bytes += o.matched_scratch_bytes;
+        init_operator_entries += o.init_operator_entries;
         inverted_index_dense_bytes += o.inverted_index_dense_bytes;
         inverted_index_sparse_bytes += o.inverted_index_sparse_bytes;
         inverted_index_dense_columns += o.inverted_index_dense_columns;
@@ -329,11 +365,12 @@
     breakdown.operator_terms_bytes = op.store->memory_bytes();
     breakdown.op_coeffs_bytes = op.op_coeffs.capacity() * sizeof(double);
     // Every representation of the state at once: the sparse scored set plus the dense vector.
     breakdown.state_coeffs_bytes = op.state_coeffs.capacity() * sizeof(double)
                                    + op.state_rows_.capacity() * sizeof(TermIndex)
                                    + op.state_vals_.capacity() * sizeof(double);
     breakdown.indexing_bytes = op.store->index_estimated_memory_bytes();
     breakdown.init_operator_bytes = unordered_flat_map_storage_bytes(op.init_op_map);
+    breakdown.init_operator_entries = op.init_op_map.size();
     breakdown.initial_state_bytes = op.initial_state.capacity() * sizeof(size_t);
     if (op.inverted_index_.has_value()) {
         breakdown.inverted_index_bytes = op.inverted_index_->memory_bytes();

cpp/tests/evolution_detail_tests.cpp

@@ -58,12 +58,19 @@ BOOST_AUTO_TEST_CASE(matched_epoch_tail_grow) {
     BOOST_TEST(!set.is_marked(3));
 }
 
-// When the epoch counter saturates uint32_t, begin_gate zero-fills and restarts so marks stay correct.
-BOOST_AUTO_TEST_CASE(matched_epoch_u32_wrap_resets) {
+// When the epoch counter saturates its stamp width, begin_gate zero-fills and restarts so marks stay
+// correct. LOAD-BEARING, together with the test below: the stamp is 16-bit (MatchedEpochSet::Stamp), so
+// this is a path real circuits take once per 65535 gates, not the effectively dead 2^32 path it was.
+// Neither test may be dropped in a refactor.
+//
+// This one reaches the boundary by ASSIGNING to cur_, which pins the branch condition and the counter
+// restart but cannot pin the std::fill: after the reset a slot holding the old maximum stamp differs
+// from the new cur_ whether or not the array was cleared. The test below is what pins the fill.
+BOOST_AUTO_TEST_CASE(matched_epoch_stamp_wrap_resets) {
     MatchedEpochSet set;
     set.begin_gate(4); // allocate the backing array
     // Force the counter to the wrap boundary; a stale slot still equals the pre-wrap counter.
-    set.cur_ = std::numeric_limits<uint32_t>::max();
+    set.cur_ = std::numeric_limits<MatchedEpochSet::Stamp>::max();
     set.mark(1);
     BOOST_TEST(set.is_marked(1));
 
@@ -74,6 +81,44 @@ BOOST_AUTO_TEST_CASE(matched_epoch_u32_wrap_resets) {
     BOOST_TEST(set.is_marked(2));
 }
 
+// Drive the wrap the way a circuit does — by counting gates — rather than by poking cur_. At a 16-bit
+// stamp the whole period is 65535 begin_gate() calls, so this is affordable; at 32 bits it was not,
+// which is why the branch never had such a test.
+//
+// LOAD-BEARING and not redundant with the test above. The mark is stamped at epoch 1 and the counter is
+// then walked all the way back round to 1, which is the ONLY arrangement in which a missing std::fill is
+// observable: the stale slot aliases exactly onto the new epoch and reads as marked. Remove the fill and
+// this test fails; remove this test and the fill becomes silently deletable.
+BOOST_AUTO_TEST_CASE(matched_epoch_stamp_wrap_reached_by_gate_count) {
+    constexpr auto kMaxStamp = std::numeric_limits<MatchedEpochSet::Stamp>::max();
+    constexpr size_t kPeriod = static_cast<size_t>(kMaxStamp); // gates between resets
+
+    MatchedEpochSet set;
+    set.begin_gate(4); // gate 1 -> cur_ == 1
+    BOOST_REQUIRE(set.cur_ == MatchedEpochSet::Stamp{1});
+    set.mark(1); // a mark stamped with epoch 1
+    BOOST_TEST(set.is_marked(1));
+
+    // Gates 2..kPeriod: one increment each, and no early wrap. Folded into a single assertion rather
+    // than 65534 of them so the loop stays cheap.
+    bool one_epoch_per_gate = true;
+    for (size_t k = 2; k <= kPeriod; ++k) {
+        set.begin_gate(4);
+        one_epoch_per_gate = one_epoch_per_gate && (static_cast<size_t>(set.cur_) == k);
+    }
+    BOOST_TEST(one_epoch_per_gate);
+    BOOST_TEST(set.cur_ == kMaxStamp); // boundary reached by counting, not by assignment
+
+    // Gate kPeriod+1 is the wrap. cur_ returns to 1 — exactly the stamp the surviving mark from gate 1
+    // carries — so is_marked(1) can only be false because the array was zero-filled.
+    set.begin_gate(4);
+    BOOST_TEST(set.cur_ == MatchedEpochSet::Stamp{1});
+    BOOST_TEST(!set.is_marked(1));
+    set.mark(2);
+    BOOST_TEST(set.is_marked(2));
+    BOOST_TEST(!set.is_marked(1));
+}
+
 BOOST_AUTO_TEST_CASE(cutoff_context_abs_coeff_for) {
     const VecD coeffs{-3.0, 2.0, 0.0};