From 213cc078b97ee96b581d165495289d143e15d580 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Mon, 17 Aug 2026 08:15:29 +0000 Subject: [PATCH 01/13] refactor(c++)!: :art: state each picture rule once in a stateless policy The Heisenberg/Schrodinger choice was one `bool schrodinger_`, derived from whether `schrodinger_cutoff` held a value, and read at 17 decision points inside `MonomialPropagator`. The gate-order rule alone was written four times, and `contract_partially` carried two near-identical 14-line arms. Introduce a picture policy pair, built exactly like the existing `MajoranaAlgebra`/`PauliAlgebra` pair: `HeisenbergPicture` and `SchrodingerPicture` in `cpp/monoprop/picture/Picture.h`, each stating one picture's whole rule set (gate traversal direction, applied-angle sign, `map_params` phase, the live coefficient vector, the contraction partner), and one `with_picture()` bridge from the runtime value to the policy type. No inheritance, no virtual member, and still one propagator class -- the bindings, `PartitionGroup` and the dispatch generators keep their shape. Every public entry point binds the policy once with `with_picture()`. The whole private layer is then templated on the policy and never re-tests which picture it is in; there is deliberately no runtime-dispatching helper layer. The constructor now takes `PictureSpec = std::variant` instead of `std::optional`, so a Heisenberg run cannot carry a state cutoff and a Schrodinger run cannot omit one. That removes a dead `value_or` fallback which could never fire. Python keeps the historical "a cutoff selects Schrodinger" spelling; the variant is built at the binding boundary. `ContractSink` and `apply_fused_contract` take the policy as a template parameter, so the two branches that sat in the hot loop become `if constexpr`. The policy is bound at the sink only: templating `build_layer` itself would multiply the `with_algebra` scan above it into four instantiations per mode width instead of two. `bool schrodinger` becomes `Picture` in `MPGraph`, `pare_graph`, `build_layer` and `MPOperator::update_initial_operator`, so call sites no longer read `false` with a comment to explain it. Also drops the two vestigial `virtual` keywords -- nothing overrode either. Results are bit-identical to the previous implementation, verified by sha256 over the exact IEEE-754 bytes of 95 expectation-value and gradient fingerprints across both pictures, both algebras and both partition counts. BREAKING CHANGE: the `MonomialPropagator` C++ constructor takes a `PictureSpec` where it took `std::optional schrodinger_cutoff`. The Python API is unchanged. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 8 + cpp/include/monoprop/MPFunctions.h | 2 +- cpp/include/monoprop/MPGraph.h | 18 +- cpp/include/monoprop/MonomialPropagator.h | 61 ++-- cpp/monoprop/CMakeLists.txt | 1 + cpp/monoprop/MPGraph.cpp | 6 +- cpp/monoprop/core/CMakeLists.txt | 1 + cpp/monoprop/core/Picture.h | 42 +++ .../detail/evolution/layer_build/Engine.h | 41 ++- .../detail/evolution/layer_build/FusedApply.h | 14 +- .../MonomialPropagator.inl | 302 +++++++++--------- cpp/monoprop/detail/operator/MPOperator.h | 5 +- cpp/monoprop/detail/pare/PareGraph.cpp | 6 +- cpp/monoprop/picture/CMakeLists.txt | 8 + cpp/monoprop/picture/Picture.h | 166 ++++++++++ cpp/tests/GraphBuildHarness.h | 4 +- cpp/tests/TestUtilities.h.in | 24 +- cpp/tests/build_graph_tests.cpp | 11 +- cpp/tests/ctor_validation_tests.cpp | 2 +- cpp/tests/exact_upper_atol_rescue.cpp | 2 +- cpp/tests/fused_cos_sweep_tests.cpp | 6 +- cpp/tests/gate_boundaries.cpp | 2 +- cpp/tests/mp_graph_tests.cpp | 18 +- cpp/tests/mp_operator_tests.cpp | 8 +- cpp/tests/mpfunctions.cpp | 8 +- .../mpi_distributed_layer_equivalence.cpp | 8 +- cpp/tests/mpi_fresh_insert_equivalence.cpp | 4 +- cpp/tests/pare_graph_tests.cpp | 2 +- cpp/tests/partition_equivalence_tests.cpp | 10 +- cpp/tests/pauli_build_layer_tests.cpp | 32 +- cpp/tests/update_initial_operator.cpp | 8 +- docs/content/docs/features/initialisation.mdx | 4 + docs/content/docs/testing.mdx | 5 +- src/monoprop/bindings/binder.h | 6 +- 34 files changed, 560 insertions(+), 285 deletions(-) create mode 100644 cpp/monoprop/core/Picture.h create mode 100644 cpp/monoprop/picture/CMakeLists.txt create mode 100644 cpp/monoprop/picture/Picture.h diff --git a/AGENTS.md b/AGENTS.md index a932c4ab..708274b5 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -67,6 +67,14 @@ Key files: (`MajoranaAlgebra`, `PauliAlgebra` in `algebra/Algebra.h`) over shared structural primitives (`algebra/AlgebraCommon.h`). The propagation backbone (the scan/fold in `detail/evolution/...`) is templated on the algebra policy and bound to a runtime `Basis` once, via `with_algebra`. +- **`Picture` / the picture policy** (`cpp/monoprop/core/Picture.h`, `cpp/monoprop/picture/Picture.h`): the two + simulation pictures are sibling models (`HeisenbergPicture`, `SchrodingerPicture`), built the same way as the + algebras. Each states one picture's whole rule set — gate traversal direction, applied-angle sign, + `map_params` phase, the live coefficient vector, the contraction partner — so no `if (schrodinger)` is + written twice. Cold sites reach a policy through the `picture_*()` helpers; the fused + `ContractSink`/`apply_fused_contract` pair takes the policy as a template parameter, bound once inside + `build_layer`. The picture is fixed at construction: the constructor takes a + `PictureSpec = std::variant`, so only a Schrodinger run carries a state cutoff. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index c7419ffe..d91a2097 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -110,7 +110,7 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, - bool schrodinger, + Picture picture, mpi::Comm comm, const std::function &full_cos_of_layer) -> MPGraph; } // namespace monoprop diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index e446c9d4..2f859447 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -21,6 +21,7 @@ #include #include +#include "monoprop/core/Picture.h" #include "monoprop/detail/graph/MPGraphLayers.h" #include "monoprop/detail/graph/MPGraphViews.h" #include "monoprop/monopropExport.h" @@ -33,7 +34,7 @@ class monoprop_EXPORT MPGraph { using LayerIterator = std::vector::iterator; using ConstLayerIterator = std::vector::const_iterator; - bool schrodinger_; + Picture picture_; std::vector layers_; size_t front_offset_ = 0; @@ -53,7 +54,10 @@ class monoprop_EXPORT MPGraph { auto active_end_iterator() const -> ConstLayerIterator { return layers_.end(); } - auto append_position() -> LayerIterator { return schrodinger_ ? active_begin_iterator() : active_end_iterator(); } + // Schrödinger prepends: it consumes the circuit front-to-back, so the newest gate is the earliest layer. + auto append_position() -> LayerIterator { + return is_schrodinger() ? active_begin_iterator() : active_end_iterator(); + } auto append_layer(Layer layer) -> void { layers_.emplace(append_position(), std::move(layer)); } @@ -65,11 +69,9 @@ class monoprop_EXPORT MPGraph { } public: - explicit MPGraph(bool schrodinger) : schrodinger_(schrodinger) {} + explicit MPGraph(Picture picture) : picture_(picture) {} - explicit MPGraph(bool schrodinger, std::vector layers) - : schrodinger_(schrodinger), - layers_(std::move(layers)) {} + explicit MPGraph(Picture picture, std::vector layers) : picture_(picture), layers_(std::move(layers)) {} /// Gate info (param_index, gen_coeff, gate_index) is written onto `storage` here while it is still /// mutable, before it is frozen into the Layer's shared const core. @@ -100,7 +102,9 @@ class monoprop_EXPORT MPGraph { /// Non-owning replay view over the active layers, in build order. auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } - auto is_schrodinger() const -> bool { return schrodinger_; } + auto picture() const -> Picture { return picture_; } + + auto is_schrodinger() const -> bool { return picture_ == Picture::Schrodinger; } /// A normally-built layer stores no cosine set, so the companion cosine-index count cannot come from /// the graph: only the operator's inverted index can supply it. diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index f80f3e87..a6151fbf 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -43,6 +43,7 @@ #include "monoprop/detail/evolution/CosineRecompute.h" #include "monoprop/detail/mpi/MPICompat.h" #include "monoprop/detail/mpi/MPIUtils.h" +#include "monoprop/picture/Picture.h" namespace monoprop { namespace detail { @@ -54,8 +55,7 @@ class PartitionGroup; } // namespace detail /// A propagator setting is out of range, or inconsistent with another setting. -// Covers a crossed atol pair and a logical width outside [1, NumModes]; also thrown from -// MonomialPropagatorImpl.h +// Covers a crossed atol pair and a logical width outside [1, NumModes]. class PropagatorConfigError : public std::runtime_error { public: using std::runtime_error::runtime_error; @@ -70,10 +70,11 @@ class MultiPartitionUnsupported : public std::runtime_error { template class MonomialPropagator { public: + /// `picture` selects Heisenberg or Schrodinger; only the Schrodinger arm carries a state cutoff. MonomialPropagator(const OperatorDict &initial_operator, unsigned int cutoff, const VecZ &initial_state, - std::optional schrodinger_cutoff, + const PictureSpec &picture, mpi::Comm comm, std::optional lower_atol = std::nullopt, std::optional upper_atol = std::nullopt, @@ -84,10 +85,10 @@ class MonomialPropagator { size_t partitions = 0); /// Out-of-line because partition_group_ is a unique_ptr to an incomplete type here. - virtual ~MonomialPropagator(); + ~MonomialPropagator(); /// Deep copy: clones the operator store, shares the immutable graph cores, and clones the whole - /// partition group on a facade. The virtual destructor suppresses implicit moves, so a "move" deep-copies. + /// partition group on a facade. Declaring it suppresses the implicit moves, so a "move" deep-copies. MonomialPropagator(const MonomialPropagator &other); auto operator=(const MonomialPropagator &) -> MonomialPropagator & = delete; @@ -210,7 +211,9 @@ class MonomialPropagator { }); } - auto schrodinger() const -> bool { return schrodinger_; } + auto picture() const -> Picture { return picture_; } + + auto schrodinger() const -> bool { return picture_ == Picture::Schrodinger; } auto basis() const -> Basis { return basis_; } @@ -274,7 +277,9 @@ class MonomialPropagator { auto evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>>; - virtual auto update_initial_operator(const OperatorDict &op_dict) -> void { apply_initial_operator_(op_dict); } + auto update_initial_operator(const OperatorDict &op_dict) -> void { + with_picture(picture_, [&]() { this->template apply_initial_operator_

(op_dict); }); + } protected: static inline const auto ev_fn = [](const EvalRequest &request, @@ -288,10 +293,11 @@ class MonomialPropagator { /// Distribute op_dict across ranks and apply this rank's share; returns its new (terms, coeffs) /// so caches can refresh. + template auto apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD>; - bool schrodinger_; - mpi::Comm comm_; // real MPI across nodes, or an in-process comm across partitions + Picture picture_; // immutable after construction: no path switches the picture mid-simulation + mpi::Comm comm_; // real MPI across nodes, or an in-process comm across partitions CutoffFn cutoff_fn_; detail::MPOperator mp_op_; MPGraph graph_; @@ -397,25 +403,33 @@ class MonomialPropagator { auto validate_cutoff_config_(CutoffType cutoff_type, const std::optional> &basis_change) const -> void; - auto initialize_operator_caches_() -> void; + // Everything below takes the picture as a policy type P, never as a runtime value: each public entry + // point binds it once with with_picture(), and the whole private layer is then written against one + // picture at a time. Nothing here re-tests which picture it is in. - auto current_picture_coeffs_() -> const VecD &; + template + auto initialize_operator_caches_() -> void; + // Grow `coeffs` to the operator's current term count, filling from the picture's live vector. + template auto extend_coeffs_from_current_picture_if_needed_(VecD &coeffs) -> void; + template auto evolve_mode_build_graph_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, const VecZ &gate_indices, std::optional only_rotate_len_k) -> void; - // Returns {build_angle, apply_angle}; apply is the build angle, negated in Schrödinger. - auto gate_angle_(const VecD &mapped_params, size_t i, size_t majoranas_size) const -> std::pair { - const size_t idx = schrodinger_ ? i : majoranas_size - 1 - i; - const double build_angle = mapped_params[idx]; - return {build_angle, schrodinger_ ? -build_angle : build_angle}; + // Returns {build_angle, apply_angle}; apply is the build angle, negated in Schrödinger. `slot` is the + // optimizer-order slot run_gate_loop_ resolved for this step -- the one place the order rule lives. + template + static auto gate_angle_(const VecD &mapped_params, size_t slot) -> std::pair { + const double build_angle = mapped_params[slot]; + return {build_angle, P::apply_sign * build_angle}; } + template auto evolve_mode_graph_with_coeffs_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, @@ -424,17 +438,23 @@ class MonomialPropagator { const VecD &operator_coeffs, std::optional only_rotate_len_k) -> void; + // build_layer resolves the same policy for its fused sink, so the cosine sweep and the apply agree. + template auto evolve_mode_contract_immediately_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, const VecD ¶meters, std::optional only_rotate_len_k) -> void; - template + // Walks the gates in simulation order, calling evolution_func(generator, only_rotate_len_k, slot). + // `slot` is the optimizer-order index the step consumes; resolving it here is what keeps the picture's + // traversal direction in a single place. + template auto run_gate_loop_(const std::vector &majoranas, std::optional only_rotate_len_k, EvolutionFunc evolution_func) -> void; + template auto propagate_one_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs = std::nullopt, @@ -445,6 +465,7 @@ class MonomialPropagator { // fused_scale_coeffs (ContractImmediately only): the picture's mutable coeff vector for the uncapped // fused cos sweep; the taken decision is reported via fused_scale so the apply matches. See build_layer. + template auto build_evolve_result_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs = std::nullopt, @@ -454,7 +475,11 @@ class MonomialPropagator { VecD *fused_scale_coeffs = nullptr, bool *fused_scale = nullptr) -> std::shared_ptr; - template + auto contract_partially_(const VecD ¶meters, bool inplace) -> VecD; + + template > auto make_functional_(Fn &&func, std::optional pare_threshold) -> std::function; diff --git a/cpp/monoprop/CMakeLists.txt b/cpp/monoprop/CMakeLists.txt index b18c8b33..46bbb0c2 100644 --- a/cpp/monoprop/CMakeLists.txt +++ b/cpp/monoprop/CMakeLists.txt @@ -113,6 +113,7 @@ target_link_libraries( add_subdirectory(algebra) add_subdirectory(core) add_subdirectory(detail) +add_subdirectory(picture) install( TARGETS diff --git a/cpp/monoprop/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index a425582d..2fd4237f 100644 --- a/cpp/monoprop/MPGraph.cpp +++ b/cpp/monoprop/MPGraph.cpp @@ -62,7 +62,7 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { const auto k = std::min(key, layers()); sliced_layers.reserve(k); - if (schrodinger_) { + if (is_schrodinger()) { const size_t active_end = active_end_index(); for (size_t i = 0; i < k; ++i) { sliced_layers.push_back(layers_[active_end - 1 - i]); @@ -84,12 +84,12 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { } } - return MPGraph(schrodinger_, std::move(sliced_layers)); + return MPGraph(picture_, std::move(sliced_layers)); } auto MPGraph::slice_view(size_t key) const -> MPGraphView { const auto k = std::min(key, layers()); - if (schrodinger_) { + if (is_schrodinger()) { return {layers_, active_end_index() - k, k, true}; } return {layers_, active_begin_index(), k, false}; diff --git a/cpp/monoprop/core/CMakeLists.txt b/cpp/monoprop/core/CMakeLists.txt index d9faa730..43c5a822 100644 --- a/cpp/monoprop/core/CMakeLists.txt +++ b/cpp/monoprop/core/CMakeLists.txt @@ -5,4 +5,5 @@ target_sources( TYPE HEADERS FILES "Monomial.h" + "Picture.h" ) diff --git a/cpp/monoprop/core/Picture.h b/cpp/monoprop/core/Picture.h new file mode 100644 index 00000000..654250a6 --- /dev/null +++ b/cpp/monoprop/core/Picture.h @@ -0,0 +1,42 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +// Deliberately dependency-free: exported headers (MPGraph.h, MPFunctions.h, MonomialPropagator.h) +// carry a Picture in their signatures. The behaviour behind each value lives in picture/Picture.h. + +#include +#include + +namespace monoprop { + +// Which object the circuit propagates. Fixed at propagator construction; no path switches it later. +enum class Picture : uint8_t { Heisenberg, Schrodinger }; + +// The picture selector as a constructor argument. Only the Schrödinger arm carries a cutoff, so a +// Heisenberg run cannot be given one, and a Schrödinger run cannot omit it. +struct Heisenberg {}; // propagate the observable backwards; the reference state is held fixed + +struct Schrodinger { + unsigned int state_cutoff; // bounds the monomial expansion of the state, read like `cutoff` +}; + +using PictureSpec = std::variant; + +inline auto kind_of(const PictureSpec &spec) -> Picture { + return std::holds_alternative(spec) ? Picture::Schrodinger : Picture::Heisenberg; +} + +} // namespace monoprop diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index e507e96d..d0c4a80e 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -34,6 +34,7 @@ #include "monoprop/detail/graph_encoding/MPGraphEncodingStorage.h" #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/picture/Picture.h" namespace monoprop::detail { @@ -182,7 +183,10 @@ struct GraphSink { // Fused ContractImmediately sink: applies each resolved rotation directly to op_coeffs via the // FusedContract record streams (no LayerCore — finalize returns nullptr). wants_values=true: the scan // captures the signed pre-cos v_src, and resolve reads v_tgt from op_coeffs (·inv_cos under the cos sweep). -template +// P is the picture policy (picture/Picture.h). It is a template parameter, not a field: the fresh +// cross-rank miss arm below sits in the per-query resolve loop, and the apply that drains these records +// (apply_fused_contract) must be specialized on the same policy. +template struct ContractSink { static constexpr bool wants_values = true; static constexpr size_t kStride = kQueryWordsFused; @@ -195,14 +199,13 @@ struct ContractSink { const VecD &op_coeffs; // the very array the scan read, not a copy bool fused_scale; // fused cos sweep active: hit v_tgt recovered as stored·inv_cos double inv_cos; - bool schrodinger; // fresh cross-rank miss coeff: 0 (Heisenberg) vs state-scored (Schrödinger) Basis basis; // Pauli vs Majorana state scoring of fresh cross-rank Schrödinger misses size_t def_base_ = 0; // deferred self-insert base into fc.inserts size_t cross_base_ = 0; // cross-rank resolver-half base into fc.cross_half - Monomial state_mask_{}; // Schrödinger fresh-insert scoring mask (empty in Heisenberg) + Monomial state_mask_{}; // Schrödinger fresh-insert scoring mask (unused in Heisenberg) - // No constructor on purpose: as an aggregate the call site names each field, so the two adjacent - // bools cannot be swapped silently. GraphSink keeps its ctor because it sizes `acc` from R. + // No constructor on purpose: as an aggregate the call site names each field. GraphSink keeps its ctor + // because it sizes `acc` from R. // Self-resolve hit (both endpoints local). always_inline: called once per surviving rotation in the // R=1 hot loop, where a real call is a measurable regression on the Pauli benches. @@ -233,7 +236,9 @@ struct ContractSink { size_t /*rank_count*/, MPOperator &op, const std::vector> & /*responses*/) -> void { - state_mask_ = schrodinger ? initial_state_mask(op.initial_state) : Monomial{}; + if constexpr (P::is_schrodinger) { + state_mask_ = initial_state_mask(op.initial_state); + } cross_base_ = fc.cross_half.size(); fc.cross_half.resize(cross_base_ + pr.nq_total); } @@ -247,7 +252,8 @@ struct ContractSink { if (ip < pr.base) { v_tgt = fused_scale ? op_coeffs[ip] * inv_cos : op_coeffs[ip]; } - else if (schrodinger) { + else if constexpr (P::is_schrodinger) { + // Fresh cross-rank insert: the state already carries a diagonal amplitude for a paired monomial. v_tgt = is_paired(pr.mono[g]) ? algebra_state_phase(basis, pr.mono[g], state_mask_) : 0.0; } @@ -532,7 +538,7 @@ auto build_layer(MPOperator &local_op, mpi::Comm comm, CosMask *out_cos = nullptr, FusedContract *fused_contract = nullptr, - bool schrodinger = false, + Picture picture = Picture::Heisenberg, VecD *fused_scale_coeffs = nullptr, bool *fused_scale_out = nullptr, Basis basis = Basis::Majorana) -> std::shared_ptr { @@ -612,14 +618,17 @@ auto build_layer(MPOperator &local_op, std::shared_ptr storage; if (use_fused) { const double inv_cos = fused_scale ? 1.0 / cos_build : 1.0; // pre-cos recovery factor for hit v_tgt - storage = run(ContractSink{.R = R, - .my_rank = my_rank, - .fc = *fused_contract, - .op_coeffs = coeffs, - .fused_scale = fused_scale, - .inv_cos = inv_cos, - .schrodinger = schrodinger, - .basis = basis}); + // The picture is bound here and nowhere higher: templating build_layer itself would multiply the + // with_algebra scan above into four instantiations per mode width instead of two. + storage = with_picture(picture, [&]() { + return run(ContractSink{.R = R, + .my_rank = my_rank, + .fc = *fused_contract, + .op_coeffs = coeffs, + .fused_scale = fused_scale, + .inv_cos = inv_cos, + .basis = basis}); + }); } else { storage = run(GraphSink{R, my_rank}); diff --git a/cpp/monoprop/detail/evolution/layer_build/FusedApply.h b/cpp/monoprop/detail/evolution/layer_build/FusedApply.h index fdd002de..e2be7bae 100644 --- a/cpp/monoprop/detail/evolution/layer_build/FusedApply.h +++ b/cpp/monoprop/detail/evolution/layer_build/FusedApply.h @@ -19,6 +19,7 @@ #include "monoprop/TypeAliases.h" #include "monoprop/detail/evolution/CosineRecompute.h" #include "monoprop/detail/evolution/layer_build/Common.h" +#include "monoprop/picture/Picture.h" namespace monoprop::detail { @@ -29,15 +30,14 @@ namespace monoprop::detail { // here; slots born after that sweep (fresh inserts) fold cos in via their apply arm below. // • two-pass (length cap / cos==0 fallback): scale_cos_mask runs here, then every arm is a plain add. // At R>1 each rank applies only the add to the slot it owns (half rotations in fc.cross_half). -inline auto apply_fused_contract(FusedContract &fc, - VecD &op_coeffs, - const CosMask &cos, - double param, - bool schrodinger, - bool fused_scale) -> void { +// P is the picture policy; it must be the one build_layer resolved for its ContractSink, or the records +// drained here were not built for this arithmetic. +template +auto apply_fused_contract(FusedContract &fc, VecD &op_coeffs, const CosMask &cos, double param, bool fused_scale) + -> void { // (1) insert records: v_tgt is the freshly-inserted term's pre-cos coeff, readable only now op_coeffs // is extended. Needed only in Schrödinger — a Heisenberg fresh insert has coeff 0, so skip the gather. - if (schrodinger) { + if constexpr (P::is_schrodinger) { for (size_t k = 0; k < fc.inserts.size(); ++k) { fc.inserts[k].v_tgt = op_coeffs[fc.inserts[k].tgt]; } diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 117590ec..00eba12b 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -74,7 +74,7 @@ template MonomialPropagator::MonomialPropagator(const OperatorDict &initial_operator, unsigned int cutoff, const VecZ &initial_state, - std::optional schrodinger_cutoff, + const PictureSpec &picture, mpi::Comm comm, std::optional lower_atol, std::optional upper_atol, @@ -83,10 +83,10 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope size_t logical_num_modes, Basis basis, size_t partitions) - : schrodinger_{schrodinger_cutoff.has_value()}, + : picture_{kind_of(picture)}, comm_{comm}, mp_op_{}, - graph_(schrodinger_cutoff.has_value()), + graph_(picture_), cutoff_{cutoff}, lower_atol_{lower_atol}, upper_atol_{upper_atol}, @@ -123,7 +123,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope return std::make_unique>(initial_operator, cutoff, initial_state, - schrodinger_cutoff, + picture, partition_comm, lower_atol, upper_atol, @@ -159,9 +159,17 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope } } - auto sc = schrodinger_cutoff.value_or(cutoff + 2); - sc = std::min(sc, static_cast(2 * logical_num_modes_)); - auto op = schrodinger_ ? generate_paired_op(sc / 2 + sc % 2, logical_num_modes_) : local_heisenberg_terms; + // The picture decides the initial monomial set: Heisenberg starts from the observable's own terms, + // Schrödinger from every paired monomial the state cutoff admits. The variant is what makes the state + // cutoff unconditionally present here -- the old optional needed a fallback that could never fire. + MonomialList op; + if (const auto *state = std::get_if(&picture)) { + const auto sc = std::min(state->state_cutoff, static_cast(2 * logical_num_modes_)); + op = generate_paired_op(sc / 2 + sc % 2, logical_num_modes_); + } + else { + op = std::move(local_heisenberg_terms); + } const size_t expected_local_terms = std::max(1, op.size() / std::max(1, num_ranks)); // Must run before the store: packed_inline_width_() derives the packed-row width from cutoff_fn_. @@ -184,7 +192,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope mp_op_.initial_state = initial_state; core_term_ = core_term; - initialize_operator_caches_(); + with_picture(picture_, [&]() { this->template initialize_operator_caches_

(); }); } template @@ -192,7 +200,7 @@ MonomialPropagator::~MonomialPropagator() = default; template MonomialPropagator::MonomialPropagator(const MonomialPropagator &other) - : schrodinger_(other.schrodinger_), + : picture_(other.picture_), comm_(other.comm_), cutoff_fn_(other.cutoff_fn_), mp_op_(other.mp_op_), @@ -335,20 +343,12 @@ auto MonomialPropagator::partitioned_graph_memory_usage_() const -> Gr template auto MonomialPropagator::packed_inline_width_() const -> size_t { - constexpr size_t kMax = detail::OperatorIndex::kMaxInlinePositions; - constexpr size_t kDefault = detail::OperatorIndex::kDefaultInlinePositions; - if (schrodinger_) { - return kDefault; - } - // The bound is already in physical slots (CutoffEvaluator::max_slot_bound), so nothing to scale. - const auto bound = detail::CutoffEvaluator(cutoff_fn_).max_slot_bound(); - if (!bound) { - return kDefault; - } - return std::min(*bound, kMax); + // A leaf: one policy question, one call site (the ctor, before any P is in scope), so it binds here. + return with_picture(picture_, [&]() { return P::template packed_inline_width(cutoff_fn_); }); } template +template auto MonomialPropagator::apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD> { ++initial_operator_epoch_; @@ -373,7 +373,7 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o } } - return mp_op_.update_initial_operator(new_op, schrodinger_); + return mp_op_.update_initial_operator(new_op, P::picture); } template @@ -496,33 +496,23 @@ auto MonomialPropagator::regenerate_cutoff_fn_() -> void { } template +template auto MonomialPropagator::initialize_operator_caches_() -> void { (void)mp_op_.get_operator(); - // Heisenberg warms the sparse state only; densifying here would defeat it. Schrödinger's dense vector - // IS the live evolved vector. - if (schrodinger_) { - (void)mp_op_.dense_state(); - } - else { - (void)mp_op_.sparse_state(); - } + P::warm_state(mp_op_); (void)mp_op_.inverted_index(); mp_op_.op_coeffs.shrink_to_fit(); mp_op_.shrink_state_to_fit(); } template -auto MonomialPropagator::current_picture_coeffs_() -> const VecD & { - return schrodinger_ ? mp_op_.dense_state() : mp_op_.get_operator(); -} - -template +template auto MonomialPropagator::extend_coeffs_from_current_picture_if_needed_(VecD &coeffs) -> void { if (coeffs.size() >= mp_op_.size()) { return; } - const auto ¤t = current_picture_coeffs_(); + const auto ¤t = P::live_coeffs(mp_op_); if (&coeffs == ¤t) { return; } @@ -534,29 +524,29 @@ auto MonomialPropagator::extend_coeffs_from_current_picture_if_needed_ } template +template auto MonomialPropagator::evolve_mode_build_graph_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, const VecZ &gate_indices, std::optional only_rotate_len_k) -> void { - const auto majoranas_size = majoranas.size(); - run_gate_loop_(majoranas, - only_rotate_len_k, - [this, ¶meter_mapping, &gen_coeffs, &gate_indices, majoranas_size](const VecZ &mono, - std::optional rot_len, - size_t i) { - const auto idx = !schrodinger_ ? majoranas_size - 1 - i : i; - propagate_one_(mono, - rot_len, - std::nullopt, - std::nullopt, - parameter_mapping[idx], - gen_coeffs[idx], - gate_indices[idx]); - }); + run_gate_loop_

(majoranas, + only_rotate_len_k, + [this, ¶meter_mapping, &gen_coeffs, &gate_indices](const VecZ &mono, + std::optional rot_len, + size_t slot) { + this->template propagate_one_

(mono, + rot_len, + std::nullopt, + std::nullopt, + parameter_mapping[slot], + gen_coeffs[slot], + gate_indices[slot]); + }); } template +template auto MonomialPropagator::evolve_mode_graph_with_coeffs_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, @@ -566,55 +556,62 @@ auto MonomialPropagator::evolve_mode_graph_with_coeffs_(const std::vec std::optional only_rotate_len_k) -> void { auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, 1.0); auto coeffs = operator_coeffs; - const auto majoranas_size = majoranas.size(); - - run_gate_loop_(majoranas, - only_rotate_len_k, - [this, ¶meter_mapping, &gen_coeffs, &gate_indices, &mapped_params, &coeffs, majoranas_size]( - const VecZ &mono, - std::optional rot_len, - size_t i) { - const auto idx = !schrodinger_ ? majoranas_size - 1 - i : i; - const auto [build_angle, apply_angle] = gate_angle_(mapped_params, i, majoranas_size); - // The cos word list is not persisted on the layer; the builder moves it out transiently. - auto cos = std::make_shared(); - auto storage = build_evolve_result_(mono, rot_len, std::cref(coeffs), build_angle, cos.get()); - graph_.append(storage, parameter_mapping[idx], gen_coeffs[idx], gate_indices[idx]); - - extend_coeffs_from_current_picture_if_needed_(coeffs); - - Layer layer(std::move(storage)); - detail::LayerCosScale cos_scale = [cos](size_t, double *c, double v) { - detail::scale_cos_mask(c, *cos, v); - }; - evolve_step(coeffs, layer, apply_angle, comm_, cos_scale); - }); + + run_gate_loop_

( + majoranas, + only_rotate_len_k, + [this, ¶meter_mapping, &gen_coeffs, &gate_indices, &mapped_params, &coeffs](const VecZ &mono, + std::optional rot_len, + size_t slot) { + const auto [build_angle, apply_angle] = gate_angle_

(mapped_params, slot); + // The cos word list is not persisted on the layer; the builder moves it out transiently. + auto cos = std::make_shared(); + auto storage = + this->template build_evolve_result_

(mono, rot_len, std::cref(coeffs), build_angle, cos.get()); + graph_.append(storage, parameter_mapping[slot], gen_coeffs[slot], gate_indices[slot]); + + this->template extend_coeffs_from_current_picture_if_needed_

(coeffs); + + Layer layer(std::move(storage)); + detail::LayerCosScale cos_scale = [cos](size_t, double *c, double v) { + detail::scale_cos_mask(c, *cos, v); + }; + evolve_step(coeffs, layer, apply_angle, comm_, cos_scale); + }); } template +template auto MonomialPropagator::evolve_mode_contract_immediately_(const std::vector &majoranas, const VecZ ¶meter_mapping, const VecD &gen_coeffs, const VecD ¶meters, std::optional only_rotate_len_k) -> void { auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, 1.0); - // Called for the side effect alone: it returns a reference to the very vector selected below. - (void)current_picture_coeffs_(); - VecD *op_coeffs = schrodinger_ ? &mp_op_.state_coeffs : &mp_op_.op_coeffs; - const auto majoranas_size = majoranas.size(); - run_gate_loop_( - majoranas, - only_rotate_len_k, - [this, &mapped_params, op_coeffs, majoranas_size](const VecZ &mono, std::optional rot_len, size_t i) { - const auto [build_angle, apply_angle] = gate_angle_(mapped_params, i, majoranas_size); - // extend_coeffs must run after build_evolve_result_'s self-rank grow and before the apply. - CosMask cos; - detail::FusedContract fc; - bool fused_scale = false; - build_evolve_result_(mono, rot_len, std::cref(*op_coeffs), build_angle, &cos, &fc, op_coeffs, &fused_scale); - extend_coeffs_from_current_picture_if_needed_(*op_coeffs); - detail::apply_fused_contract(fc, *op_coeffs, cos, apply_angle, schrodinger_, fused_scale); - }); + // Called for the side effect alone: it materializes the very vector the slot below points at. + (void)P::live_coeffs(mp_op_); + VecD *op_coeffs = &P::live_coeffs_slot(mp_op_); + run_gate_loop_

(majoranas, + only_rotate_len_k, + [this, &mapped_params, op_coeffs](const VecZ &mono, std::optional rot_len, size_t slot) { + const auto [build_angle, apply_angle] = gate_angle_

(mapped_params, slot); + // extend_coeffs must run after build_evolve_result_'s self-rank grow and before the apply. + CosMask cos; + detail::FusedContract fc; + bool fused_scale = false; + this->template build_evolve_result_

(mono, + rot_len, + std::cref(*op_coeffs), + build_angle, + &cos, + &fc, + op_coeffs, + &fused_scale); + this->template extend_coeffs_from_current_picture_if_needed_

(*op_coeffs); + // build_layer resolves P again from the same Picture value for its sink; the fused cosine + // sweep and this apply must agree. + detail::apply_fused_contract

(fc, *op_coeffs, cos, apply_angle, fused_scale); + }); } template @@ -650,10 +647,15 @@ auto MonomialPropagator::build_graph(const std::vector &majorana g += gate_offset; } - if (!parameters.has_value()) { - evolve_mode_build_graph_(majoranas, parameter_mapping, gen_coeffs, local_gates, only_rotate_len_k); - } - else { + with_picture(picture_, [&]() { + if (!parameters.has_value()) { + this->template evolve_mode_build_graph_

(majoranas, + parameter_mapping, + gen_coeffs, + local_gates, + only_rotate_len_k); + return; + } // map_params() indexes `parameters` by parameter_mapping, so a too-short vector reads out of bounds. validate_parameters_length(*parameters, parameter_mapping); // Coefficient-informed build: seed by contracting the existing graph so atol truncation sees @@ -674,19 +676,19 @@ auto MonomialPropagator::build_graph(const std::vector &majorana parameters->size())); } const VecD existing_params(parameters->begin(), parameters->begin() + static_cast(m)); - seed = contract_partially(existing_params, false); + seed = this->template contract_partially_

(existing_params, false); } else { - seed = current_picture_coeffs_(); + seed = P::live_coeffs(mp_op_); } - evolve_mode_graph_with_coeffs_(majoranas, - parameter_mapping, - gen_coeffs, - local_gates, - *parameters, - seed, - only_rotate_len_k); - } + this->template evolve_mode_graph_with_coeffs_

(majoranas, + parameter_mapping, + gen_coeffs, + local_gates, + *parameters, + seed, + only_rotate_len_k); + }); } template @@ -714,25 +716,31 @@ auto MonomialPropagator::propagate(const std::vector &majoranas, "build_graph() to extend it.", graph_layers())); } - evolve_mode_contract_immediately_(majoranas, parameter_mapping, gen_coeffs, parameters, only_rotate_len_k); + with_picture(picture_, [&]() { + this->template evolve_mode_contract_immediately_

(majoranas, + parameter_mapping, + gen_coeffs, + parameters, + only_rotate_len_k); + }); } template -template +template auto MonomialPropagator::run_gate_loop_(const std::vector &majoranas, std::optional only_rotate_len_k, EvolutionFunc evolution_func) -> void { // Serial per partition; parallelism comes from partitioning the operator across cores. for (size_t i = 0; i < majoranas.size(); ++i) { - const auto idx = !schrodinger_ ? majoranas.size() - 1 - i : i; - const auto &mono = majoranas[idx]; - evolution_func(mono, only_rotate_len_k, i); + const auto slot = P::gate_slot(i, majoranas.size()); + evolution_func(majoranas[slot], only_rotate_len_k, slot); } - initialize_operator_caches_(); + initialize_operator_caches_

(); } template +template auto MonomialPropagator::build_evolve_result_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs, @@ -758,13 +766,14 @@ auto MonomialPropagator::build_evolve_result_(const VecZ &gen_vec, comm_, out_cos, fused_contract, - schrodinger_, + P::picture, fused_scale_coeffs, fused_scale, basis_); } template +template auto MonomialPropagator::propagate_one_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs, @@ -772,7 +781,10 @@ auto MonomialPropagator::propagate_one_(const VecZ &gen_vec, size_t param_index, double gen_coeff, size_t gate_index) -> void { - graph_.append(build_evolve_result_(gen_vec, only_rotate_len_k, coeffs, param), param_index, gen_coeff, gate_index); + graph_.append(this->template build_evolve_result_

(gen_vec, only_rotate_len_k, coeffs, param), + param_index, + gen_coeff, + gate_index); } template @@ -916,7 +928,7 @@ auto build_cos_callbacks(const detail::InvertedIndex &inverted_index, } template -template +template auto MonomialPropagator::make_functional_(Fn &&func, std::optional pare_threshold) -> std::function { auto gate_arrays = graph_gate_arrays_(); @@ -930,13 +942,7 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optionalcore_term(); const auto comm = comm_; @@ -958,10 +964,8 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optional(combined); }; // Threshold the picture's driving vector: the Hamiltonian in Schrödinger, the state otherwise. - const auto keep = schrodinger_ ? indices_above(op, *pare_threshold) : state.indices_above(*pare_threshold); - const auto count = schrodinger_ ? op.size() : state.length(); - graph = - std::make_shared(pare_graph(graph_, keep, count, schrodinger_, comm_, full_cos_of_layer)); + const auto [keep, count] = P::pare_seed(state, op, *pare_threshold); + graph = std::make_shared(pare_graph(graph_, keep, count, P::picture, comm_, full_cos_of_layer)); } else { graph = std::shared_ptr(std::shared_ptr{}, &graph_); @@ -1013,7 +1017,8 @@ auto MonomialPropagator::expectation_value_functional(std::optional(r)](params); })[0]; }; } - return make_functional_(ev_fn, pare_threshold); + return with_picture(picture_, + [&]() { return this->template make_functional_

(ev_fn, pare_threshold); }); } template @@ -1028,7 +1033,9 @@ auto MonomialPropagator::expectation_value_and_gradient_functional(std [&](int r) { return (*fns)[static_cast(r)](params); })[0]; }; } - return make_functional_(ev_and_grad_fn, pare_threshold); + return with_picture(picture_, [&]() { + return this->template make_functional_

(ev_and_grad_fn, pare_threshold); + }); } template @@ -1054,46 +1061,41 @@ auto MonomialPropagator::contract_partially(const VecD ¶meters, bo if (partition_group_) { return concat_partitions_([&](MonomialPropagator &s) { return s.contract_partially(parameters, inplace); }); } + return with_picture(picture_, + [&]() { return this->template contract_partially_

(parameters, inplace); }); +} + +template +template +auto MonomialPropagator::contract_partially_(const VecD ¶meters, bool inplace) -> VecD { const auto gate_arrays = graph_gate_arrays_(); const auto ¶meter_mapping = gate_arrays.first; const auto &gen_coeffs = gate_arrays.second; validate_parameters_length(parameters, parameter_mapping); if (parameters.empty()) { - return current_picture_coeffs_(); + return P::live_coeffs(mp_op_); } const size_t num_majoranas = parameter_mapping.size(); + // The pictures differ in three values only: the source vector, the (phase, reverse) map_params pair, + // and the slot that receives an inplace result. Everything else -- and the order of every flop -- is shared. + const VecD &source = P::live_coeffs(mp_op_); + const auto mapped_params = + map_params(parameters, parameter_mapping, gen_coeffs, P::contract_phase, P::contract_reverse); + // Inplace slicing produces an owned MPGraph that must be bound to a named local before viewing // (never view a temporary); slice_view() views this graph's still-live layers directly. - if (schrodinger_) { - const auto &state = mp_op_.dense_state(); - const auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, -1.0); - VecD evolved_state; - if (inplace) { - const MPGraph sliced = graph_.slice_graph(num_majoranas, true); - evolved_state = evolve_operator_with_recompute_(VecD(state), sliced.replay_view(), mapped_params); - mp_op_.state_coeffs = evolved_state; - } - else { - evolved_state = - evolve_operator_with_recompute_(VecD(state), graph_.slice_view(num_majoranas), mapped_params); - } - return evolved_state; - } - - const auto &op = mp_op_.get_operator(); - const auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, 1.0, true); - VecD evolved_op; + VecD evolved; if (inplace) { const MPGraph sliced = graph_.slice_graph(num_majoranas, true); - evolved_op = evolve_operator_with_recompute_(VecD(op), sliced.replay_view(), mapped_params); - mp_op_.op_coeffs = evolved_op; + evolved = evolve_operator_with_recompute_(VecD(source), sliced.replay_view(), mapped_params); + P::live_coeffs_slot(mp_op_) = evolved; } else { - evolved_op = evolve_operator_with_recompute_(VecD(op), graph_.slice_view(num_majoranas), mapped_params); + evolved = evolve_operator_with_recompute_(VecD(source), graph_.slice_view(num_majoranas), mapped_params); } - return evolved_op; + return evolved; } template diff --git a/cpp/monoprop/detail/operator/MPOperator.h b/cpp/monoprop/detail/operator/MPOperator.h index 6d4e70da..06243885 100644 --- a/cpp/monoprop/detail/operator/MPOperator.h +++ b/cpp/monoprop/detail/operator/MPOperator.h @@ -28,6 +28,7 @@ #include "monoprop/TypeAliases.h" #include "monoprop/Utilities.h" +#include "monoprop/core/Picture.h" #include "monoprop/detail/operator/InvertedIndex.h" #include "monoprop/detail/operator/OperatorIndex.h" @@ -190,7 +191,7 @@ struct MPOperator { // Heisenberg rejects a term absent from both (new monomials may have no graph paths); Schrödinger // admits them freely (the state was already evolved). Returns the supplied terms with their encoded // coefficients, in order. - auto update_initial_operator(const OperatorDict &op_dict, bool schrodinger) + auto update_initial_operator(const OperatorDict &op_dict, Picture picture) -> std::pair, VecD> { MonomialMap new_op_map; std::pair, VecD> new_grad_op; @@ -203,7 +204,7 @@ struct MPOperator { const auto rank_init_op = init_op_map.find(mono); const auto coeff = algebra_encode_coeff(basis, v, mono); - if (!schrodinger) { + if (picture == Picture::Heisenberg) { if (rank_init_op != init_op_map.end()) { new_op_map[mono] = coeff; } diff --git a/cpp/monoprop/detail/pare/PareGraph.cpp b/cpp/monoprop/detail/pare/PareGraph.cpp index 4d6dccec..8ca382c5 100644 --- a/cpp/monoprop/detail/pare/PareGraph.cpp +++ b/cpp/monoprop/detail/pare/PareGraph.cpp @@ -103,7 +103,7 @@ auto mark_cross_rank_endpoints_kept(const LayerTraversal &layer, size_t my_rank, auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, - bool schrodinger, + Picture picture, mpi::Comm comm, const std::function &full_cos_of_layer) -> MPGraph { const size_t num_layers = graph.layers(); @@ -122,7 +122,7 @@ auto pare_graph(const MPGraph &graph, // mark_cross_rank_endpoints_kept), so nodes_to_keep stays consistent across ranks with no exchange. // Cross-rank lists are never pruned; the keep-set only has to be right so cos pruning stays exact. for (size_t iter = 0; iter < num_layers; ++iter) { - const size_t layer_idx = schrodinger ? iter : (num_layers - 1 - iter); + const size_t layer_idx = picture == Picture::Schrodinger ? iter : (num_layers - 1 - iter); const auto &layer = graph.get_layer(layer_idx); const auto lt = layer.traversal(); @@ -135,7 +135,7 @@ auto pare_graph(const MPGraph &graph, layers[layer_idx] = preserves ? Layer(layer.shared_core()) : Layer(layer.shared_core(), std::move(filtered)); } - return MPGraph(graph.is_schrodinger(), std::move(layers)); + return MPGraph(graph.picture(), std::move(layers)); } } // namespace monoprop diff --git a/cpp/monoprop/picture/CMakeLists.txt b/cpp/monoprop/picture/CMakeLists.txt new file mode 100644 index 00000000..3fac6f84 --- /dev/null +++ b/cpp/monoprop/picture/CMakeLists.txt @@ -0,0 +1,8 @@ +target_sources( + monoprop + PUBLIC + FILE_SET headers + TYPE HEADERS + FILES + "Picture.h" +) diff --git a/cpp/monoprop/picture/Picture.h b/cpp/monoprop/picture/Picture.h new file mode 100644 index 00000000..75758b1d --- /dev/null +++ b/cpp/monoprop/picture/Picture.h @@ -0,0 +1,166 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +// Each model answers the same fixed set of questions about a simulation picture: which way the gate +// loop walks the circuit, the sign an angle carries when it is applied, which coefficient vector the +// gates mutate, and which vector is the contraction partner. Sibling models, as MajoranaAlgebra and +// PauliAlgebra are -- see algebra/Algebra.h, whose with_algebra bridge this file mirrors. +// +// The models are not templates: NumModes is deduced per member, so with_picture needs no width and a +// call site reads P::live_coeffs(mp_op_). +// +// There is no runtime-dispatching helper layer on purpose. Each public entry point of MonomialPropagator +// binds the policy once with with_picture(); its whole private layer is then written against one picture +// and never re-tests which one it is. + +#include +#include +#include + +#include "monoprop/MPFunctions.h" +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/core/Picture.h" +#include "monoprop/detail/operator/MPOperator.h" + +namespace monoprop { + +struct HeisenbergPicture { + static constexpr Picture picture = Picture::Heisenberg; + static constexpr bool is_schrodinger = false; + + // Simulation step i consumes optimizer slot n-1-i: the observable walks the circuit backwards. + static auto gate_slot(size_t i, size_t n) -> size_t { return n - 1 - i; } + static constexpr double apply_sign = 1.0; // the applied angle is the build angle + + // map_params() arguments for contract_partially: forward phase, written in reverse. + static constexpr double contract_phase = 1.0; + static constexpr bool contract_reverse = true; + + // The live vector the gates mutate, and the slot it lives in. + template + static auto live_coeffs(detail::MPOperator &op) -> const VecD & { + return op.get_operator(); + } + template + static auto live_coeffs_slot(detail::MPOperator &op) -> VecD & { + return op.op_coeffs; + } + + // Warms the sparse state only; densifying here would defeat it. + template + static auto warm_state(detail::MPOperator &op) -> void { + (void)op.sparse_state(); + } + + // Energy only dots the state against the evolved operator, and the gradient scatters it into its + // own scratch, so the sparse scores are enough. + template + static auto eval_state(detail::MPOperator &op, size_t num_terms) -> EvalState { + const auto sparse = op.sparse_state(); + return EvalState::sparse(num_terms, sparse.rows, sparse.values); + } + + // The paring keep-set, thresholded on the picture's driving vector: (keep-set, local index count). + static auto pare_seed(const EvalState &state, const VecD & /*op*/, double threshold) -> std::pair { + return {state.indices_above(threshold), state.length()}; + } + + // A perf hint, never a correctness constraint: overflow spills losslessly. The bound is already in + // physical slots (CutoffEvaluator::max_slot_bound), so nothing to scale. NumModes is explicit + // everywhere it names a CutoffFn: 2*NumModes inside Monomial is a non-deduced context. + template + static auto packed_inline_width(const CutoffFn &cutoff_fn) -> size_t { + constexpr size_t kMax = detail::OperatorIndex::kMaxInlinePositions; + constexpr size_t kDefault = detail::OperatorIndex::kDefaultInlinePositions; + const auto bound = detail::CutoffEvaluator(cutoff_fn).max_slot_bound(); + if (!bound) { + return kDefault; + } + return std::min(*bound, kMax); + } +}; + +struct SchrodingerPicture { + static constexpr Picture picture = Picture::Schrodinger; + static constexpr bool is_schrodinger = true; + + // Simulation step i consumes optimizer slot i: the state walks the circuit front-to-back. + static auto gate_slot(size_t i, size_t /*n*/) -> size_t { return i; } + static constexpr double apply_sign = -1.0; // the applied angle is the negated build angle + + static constexpr double contract_phase = -1.0; + static constexpr bool contract_reverse = false; + + // The dense state IS the live evolved vector here, so it is both the source and the slot. + template + static auto live_coeffs(detail::MPOperator &op) -> const VecD & { + return op.dense_state(); + } + template + static auto live_coeffs_slot(detail::MPOperator &op) -> VecD & { + return op.state_coeffs; + } + + template + static auto warm_state(detail::MPOperator &op) -> void { + (void)op.dense_state(); + } + + // Snapshotted whole: dense_state() returns the live vector, which evolution then mutates. + template + static auto eval_state(detail::MPOperator &op, size_t /*num_terms*/) -> EvalState { + return EvalState::dense(op.dense_state()); + } + + static auto pare_seed(const EvalState & /*state*/, const VecD &op, double threshold) -> std::pair { + return {indices_above(op, threshold), op.size()}; + } + + // The state's monomials come from generate_paired_op(), not from cutoff_fn_, so the cutoff carries + // no structural bound on them. + template + static auto packed_inline_width(const CutoffFn & /*cutoff_fn*/) -> size_t { + return detail::OperatorIndex::kDefaultInlinePositions; + } +}; + +// Shape check only: the members the propagator actually calls are enforced by use, not by this concept. +template +concept PicturePolicy = requires { + { P::picture } -> std::convertible_to; + { P::is_schrodinger } -> std::convertible_to; + { P::apply_sign } -> std::convertible_to; + { P::contract_phase } -> std::convertible_to; + { P::contract_reverse } -> std::convertible_to; +}; + +static_assert(PicturePolicy); +static_assert(PicturePolicy); + +// The one runtime->policy branch, taken once per public call. decltype(auto), not auto, so a policy that +// hands back a reference into the operator does not decay to a copy; both arms must then deduce the same +// type. +template +auto with_picture(Picture picture, F &&f) -> decltype(auto) { + if (picture == Picture::Schrodinger) { + return std::forward(f).template operator()(); + } + return std::forward(f).template operator()(); +} + +} // namespace monoprop diff --git a/cpp/tests/GraphBuildHarness.h b/cpp/tests/GraphBuildHarness.h index 44b0f7de..5a862dee 100644 --- a/cpp/tests/GraphBuildHarness.h +++ b/cpp/tests/GraphBuildHarness.h @@ -35,8 +35,8 @@ inline auto layer_with_gate(std::size_t gate_index) -> monoprop::Layer { } // Built via append() so layer ordering (Heisenberg back-append, Schrödinger front-insert) is production's. -inline auto graph_with_gates(bool schrodinger, std::size_t n) -> monoprop::MPGraph { - monoprop::MPGraph graph(schrodinger); +inline auto graph_with_gates(monoprop::Picture picture, std::size_t n) -> monoprop::MPGraph { + monoprop::MPGraph graph(picture); for (std::size_t i = 0; i < n; ++i) { graph.append(std::make_shared(), /*param_index=*/0, /*gen_coeff=*/0.0, /*gate_index=*/i); } diff --git a/cpp/tests/TestUtilities.h.in b/cpp/tests/TestUtilities.h.in index 37f590ed..08d3277d 100644 --- a/cpp/tests/TestUtilities.h.in +++ b/cpp/tests/TestUtilities.h.in @@ -77,8 +77,13 @@ inline auto load_case_data(const std::string& filename) -> CaseData { return load_case(data_path); } +inline auto picture_label(const PictureSpec& picture) -> std::string { + const auto* state = std::get_if(&picture); + return state ? std::format("schrodinger({})", state->state_cutoff) : "heisenberg"; +} + struct SimulatorConfig { - std::optional schrodinger_cutoff = std::nullopt; + PictureSpec picture = Heisenberg{}; MPI_Comm comm = MPI_COMM_SELF; std::optional atol = std::nullopt; std::optional upper_atol = std::nullopt; @@ -92,7 +97,7 @@ inline auto build_simulator(const CaseData& data, const SimulatorConfig& cfg = { return MonomialPropagator(data.hamiltonian, cutoff, data.initial_state, - cfg.schrodinger_cutoff, + cfg.picture, cfg.comm, cfg.atol, cfg.upper_atol, @@ -130,8 +135,7 @@ inline auto test_evolve_build_graph(const CaseData& data, const SimulatorConfig& const std::optional pare_threshold = pare ? std::optional{1e-10} : std::nullopt; auto expval_fn = mp.expectation_value_functional(pare_threshold); double expval = expval_fn(data.parameters); - BOOST_TEST_CONTEXT("n_modes=" << n_modes << " pare=" << pare << " sch_cutoff=" - << (cfg.schrodinger_cutoff ? std::to_string(*cfg.schrodinger_cutoff) : "none")) { + BOOST_TEST_CONTEXT("n_modes=" << n_modes << " pare=" << pare << " picture=" << picture_label(cfg.picture)) { check_expval_close("Expectation Value Build Graph", expval, exact_expval); } } @@ -149,15 +153,14 @@ inline auto test_evolve_build_graph_with_coeffs(const CaseData& data, const std::optional pare_threshold = pare ? std::optional{1e-10} : std::nullopt; auto expval_fn = mp.expectation_value_functional(pare_threshold); double expval = expval_fn(data.parameters); - BOOST_TEST_CONTEXT("n_modes=" << n_modes << " pare=" << pare << " sch_cutoff=" - << (cfg.schrodinger_cutoff ? std::to_string(*cfg.schrodinger_cutoff) : "none")) { + BOOST_TEST_CONTEXT("n_modes=" << n_modes << " pare=" << pare << " picture=" << picture_label(cfg.picture)) { check_expval_close("Expectation Value Build Graph with coeffs", expval, exact_expval); } } // As above, but split across two build_graph calls so the second one lands on an already non-empty // graph -- the only path that reaches build_graph's contract_partially() seeding branch. Schrodinger -// only, so `cfg.schrodinger_cutoff` must be set: a Heisenberg build consumes each call back-to-front, +// only, so `cfg.picture` must hold a Schrodinger: a Heisenberg build consumes each call back-to-front, // so a forward split is not equivalent to one call. Each call's `parameters` covers the prefix its own // mapping reaches, which is what the seeding guard demands of the second call. template @@ -223,10 +226,11 @@ struct LihFixture { inline constexpr std::array ds_pare_values{false, true}; inline constexpr std::array ds_schrodinger_enabled{false, true}; -inline auto make_schrodinger_cutoff(bool enabled, int cutoff, int offset = 2) -> std::optional { +// The Schrodinger arm needs a looser state cutoff than the structural one to match Heisenberg. +inline auto make_picture(bool enabled, unsigned int cutoff, unsigned int offset = 2) -> PictureSpec { if (!enabled) { - return std::nullopt; + return Heisenberg{}; } - return cutoff + offset; + return Schrodinger{cutoff + offset}; } } // namespace test_utils diff --git a/cpp/tests/build_graph_tests.cpp b/cpp/tests/build_graph_tests.cpp index 83e36d2b..09f00b20 100644 --- a/cpp/tests/build_graph_tests.cpp +++ b/cpp/tests/build_graph_tests.cpp @@ -27,9 +27,8 @@ BOOST_DATA_TEST_CASE_F(ExampleDataFix, bdata::make(ds_pare_values) ^ bdata::make(ds_schrodinger_enabled), pare, sch_enabled) { - const auto schrodinger_cutoff = make_schrodinger_cutoff(sch_enabled, cutoff); SimulatorConfig cfg{ - .schrodinger_cutoff = schrodinger_cutoff ? std::optional(*schrodinger_cutoff) : std::nullopt, + .picture = make_picture(sch_enabled, cutoff), .cutoff_type = cutoff_type, .basis_change = basis_change, }; @@ -41,9 +40,8 @@ BOOST_DATA_TEST_CASE_F(ExampleDataFix, bdata::make(ds_pare_values) ^ bdata::make(ds_schrodinger_enabled), pare, sch_enabled) { - const auto schrodinger_cutoff = make_schrodinger_cutoff(sch_enabled, cutoff); SimulatorConfig cfg{ - .schrodinger_cutoff = schrodinger_cutoff ? std::optional(*schrodinger_cutoff) : std::nullopt, + .picture = make_picture(sch_enabled, cutoff), .cutoff_type = cutoff_type, .basis_change = basis_change, }; @@ -52,9 +50,8 @@ BOOST_DATA_TEST_CASE_F(ExampleDataFix, // Schrodinger-only by construction; the reason is on test_evolve_build_graph_with_coeffs_extend. BOOST_DATA_TEST_CASE_F(ExampleDataFix, build_graph_with_coeffs_extend_cases, bdata::make(ds_pare_values), pare) { - const auto schrodinger_cutoff = make_schrodinger_cutoff(/*enabled=*/true, cutoff); SimulatorConfig cfg{ - .schrodinger_cutoff = std::optional(*schrodinger_cutoff), + .picture = make_picture(/*enabled=*/true, cutoff), .cutoff_type = cutoff_type, .basis_change = basis_change, }; @@ -67,7 +64,7 @@ BOOST_AUTO_TEST_CASE(graph_size_reports_real_cosine_only_count) { const auto data = test_utils::load_case_data("random_exact.msgpack"); const auto sized = [&](unsigned int cutoff) { - auto sim = MonomialPropagator(data.hamiltonian, cutoff, data.initial_state, std::nullopt, MPI_COMM_SELF); + auto sim = MonomialPropagator(data.hamiltonian, cutoff, data.initial_state, Heisenberg{}, MPI_COMM_SELF); sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); return sim.graph_size(); }; diff --git a/cpp/tests/ctor_validation_tests.cpp b/cpp/tests/ctor_validation_tests.cpp index 2eddcd19..0557345d 100644 --- a/cpp/tests/ctor_validation_tests.cpp +++ b/cpp/tests/ctor_validation_tests.cpp @@ -46,7 +46,7 @@ auto make(const OperatorDict &op, return MP(op, cutoff, VecZ{}, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, lower_atol, upper_atol, diff --git a/cpp/tests/exact_upper_atol_rescue.cpp b/cpp/tests/exact_upper_atol_rescue.cpp index 3243d965..7462b372 100644 --- a/cpp/tests/exact_upper_atol_rescue.cpp +++ b/cpp/tests/exact_upper_atol_rescue.cpp @@ -40,7 +40,7 @@ auto build_zero_cutoff_full_rescue(const CaseData& data, MPI_Comm comm) -> Monom return MonomialPropagator(data.hamiltonian, /*cutoff=*/0U, data.initial_state, - /*schrodinger_cutoff=*/std::nullopt, + /*picture=*/Heisenberg{}, comm, /*atol=*/std::nullopt, /*upper_atol=*/std::optional{0.0}, diff --git a/cpp/tests/fused_cos_sweep_tests.cpp b/cpp/tests/fused_cos_sweep_tests.cpp index a8d36ea2..e00bd483 100644 --- a/cpp/tests/fused_cos_sweep_tests.cpp +++ b/cpp/tests/fused_cos_sweep_tests.cpp @@ -71,9 +71,11 @@ BOOST_FIXTURE_TEST_CASE(fused_sweep_matches_graph_replay_heisenberg_atol, Exampl // Schrödinger: fresh inserts are born after the sweep with a nonzero coeff, so the apply's insert arm // must fold the gate's cos into those slots itself. BOOST_FIXTURE_TEST_CASE(fused_sweep_matches_graph_replay_schrodinger, ExampleDataFix) { - check_agreement(data, SimulatorConfig{.schrodinger_cutoff = 2 * n_modes}, "schrodinger"); + check_agreement(data, SimulatorConfig{.picture = Schrodinger{2 * n_modes}}, "schrodinger"); } BOOST_FIXTURE_TEST_CASE(fused_sweep_matches_graph_replay_schrodinger_atol, ExampleDataFix) { - check_agreement(data, SimulatorConfig{.schrodinger_cutoff = 2 * n_modes, .atol = 1e-10}, "schrodinger atol=1e-10"); + check_agreement(data, + SimulatorConfig{.picture = Schrodinger{2 * n_modes}, .atol = 1e-10}, + "schrodinger atol=1e-10"); } diff --git a/cpp/tests/gate_boundaries.cpp b/cpp/tests/gate_boundaries.cpp index 0fa4b601..f2cc869a 100644 --- a/cpp/tests/gate_boundaries.cpp +++ b/cpp/tests/gate_boundaries.cpp @@ -35,7 +35,7 @@ auto make_sim() -> MonomialPropagator { return MonomialPropagator(ham, 2 * kModes, initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 2f7e6bb2..ee9cb422 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -29,7 +29,7 @@ using test_utils::graph_with_gates; using test_utils::layer_with_gate; BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_heisenberg_prefix_no_contract) { - auto graph = graph_with_gates(/*schrodinger=*/false, 5); // layers_ = [0,1,2,3,4] + auto graph = graph_with_gates(Picture::Heisenberg, 5); // layers_ = [0,1,2,3,4] auto sliced = graph.slice_graph(3, /*contract=*/false); BOOST_REQUIRE_EQUAL(sliced.layers(), 3U); @@ -43,7 +43,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_heisenberg_prefix_no_contract) { BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_schrodinger_contract_newest_first_copy_and_resize) { // Schrödinger stores newest-first: appending gates 0..4 gives layers_ = [4,3,2,1,0]. - auto graph = graph_with_gates(/*schrodinger=*/true, 5); + auto graph = graph_with_gates(Picture::Schrodinger, 5); auto sliced = graph.slice_graph(2, /*contract=*/true); // sliced = layers_[active_end-1-i] = layers_[4], layers_[3] = gates 0, 1 (oldest-first). @@ -59,7 +59,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_schrodinger_contract_newest_first_copy } BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_key_clamped_to_size) { - auto graph = graph_with_gates(/*schrodinger=*/false, 3); + auto graph = graph_with_gates(Picture::Heisenberg, 3); auto sliced = graph.slice_graph(100, /*contract=*/false); BOOST_CHECK_EQUAL(sliced.layers(), 3U); } @@ -67,7 +67,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_key_clamped_to_size) { // The maybe_compact_layers arms below are reached through Heisenberg slice_graph(contract=true). BOOST_AUTO_TEST_CASE(mp_graph_contract_clear_arm_when_prefix_covers_all) { - auto graph = graph_with_gates(/*schrodinger=*/false, 5); + auto graph = graph_with_gates(Picture::Heisenberg, 5); (void)graph.slice_graph(5, /*contract=*/true); // front_offset == size -> clear BOOST_CHECK_EQUAL(graph.layers(), 0U); // Graph is still usable after a full clear. @@ -77,7 +77,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_clear_arm_when_prefix_covers_all) { } BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { - auto graph = graph_with_gates(/*schrodinger=*/false, 100); + auto graph = graph_with_gates(Picture::Heisenberg, 100); (void)graph.slice_graph(3, /*contract=*/true); // front_offset 3 < 4096 -> no physical compaction BOOST_REQUIRE_EQUAL(graph.layers(), 97U); BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 3U); @@ -86,7 +86,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { BOOST_AUTO_TEST_CASE(mp_graph_contract_erase_arm_above_threshold) { // The erase arm fires only when front_offset >= 4096 AND 2*front_offset >= size. - auto graph = graph_with_gates(/*schrodinger=*/false, 8200); + auto graph = graph_with_gates(Picture::Heisenberg, 8200); auto sliced = graph.slice_graph(4100, /*contract=*/true); BOOST_CHECK_EQUAL(sliced.layers(), 4100U); BOOST_CHECK_EQUAL(sliced.get_layer_traversal(0).gate_index(), 0U); @@ -97,7 +97,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_erase_arm_above_threshold) { } BOOST_AUTO_TEST_CASE(mp_graph_slice_view_heisenberg_forward_window) { - auto graph = graph_with_gates(/*schrodinger=*/false, 5); + auto graph = graph_with_gates(Picture::Heisenberg, 5); auto view = graph.slice_view(3); BOOST_REQUIRE_EQUAL(view.layers(), 3U); BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); @@ -108,7 +108,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_view_heisenberg_forward_window) { BOOST_AUTO_TEST_CASE(mp_graph_slice_view_schrodinger_reversed_window) { // layers_ = [4,3,2,1,0]; slice_view(3) uses base=active_end-3=2, reverse=true. // get_layer_traversal(i) -> layers_[2 + (3-1-i)] -> gates 0,1,2 in replay order. - auto graph = graph_with_gates(/*schrodinger=*/true, 5); + auto graph = graph_with_gates(Picture::Schrodinger, 5); auto view = graph.slice_view(3); BOOST_REQUIRE_EQUAL(view.layers(), 3U); BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); @@ -133,7 +133,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { } BOOST_AUTO_TEST_CASE(mp_graph_get_layer_out_of_range_throws) { - auto graph = graph_with_gates(/*schrodinger=*/false, 3); + auto graph = graph_with_gates(Picture::Heisenberg, 3); BOOST_CHECK_NO_THROW((void)graph.get_layer(2)); BOOST_CHECK_THROW((void)graph.get_layer(3), std::out_of_range); // const overload takes the same guard. diff --git a/cpp/tests/mp_operator_tests.cpp b/cpp/tests/mp_operator_tests.cpp index d0a4698d..f6f20097 100644 --- a/cpp/tests/mp_operator_tests.cpp +++ b/cpp/tests/mp_operator_tests.cpp @@ -182,7 +182,7 @@ BOOST_AUTO_TEST_CASE(mp_operator_update_initial_operator_heisenberg_branches_pau dict[VecZ{0, 2}] = cd(1.5, 0.0); // present in store -> row coeff dict[VecZ{4, 6}] = cd(2.5, 0.0); // in init_op_map -> stays pending - const auto grad = op.update_initial_operator(dict, /*schrodinger=*/false); + const auto grad = op.update_initial_operator(dict, Picture::Heisenberg); BOOST_REQUIRE_EQUAL(op.op_coeffs.size(), 1U); BOOST_CHECK_EQUAL(op.op_coeffs[0], encode_pauli_coeff(cd(1.5, 0.0))); // Pauli encode path BOOST_CHECK(op.init_op_map.find(indices_to_bitset<8>({4, 6})) != op.init_op_map.end()); @@ -195,7 +195,7 @@ BOOST_AUTO_TEST_CASE(mp_operator_update_initial_operator_heisenberg_rejects_abse OperatorDict dict; dict[VecZ{1, 3, 5}] = cd(1.0, 0.0); // absent from both store and init_op_map - BOOST_CHECK_THROW(op.update_initial_operator(dict, /*schrodinger=*/false), std::runtime_error); + BOOST_CHECK_THROW(op.update_initial_operator(dict, Picture::Heisenberg), std::runtime_error); } BOOST_AUTO_TEST_CASE(mp_operator_update_initial_operator_schrodinger_admits_absent_term) { @@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(mp_operator_update_initial_operator_schrodinger_admits_abse OperatorDict dict; const auto fresh = indices_to_bitset<8>({1, 3, 5}); dict[VecZ{1, 3, 5}] = cd(4.0, 0.0); - op.update_initial_operator(dict, /*schrodinger=*/true); + op.update_initial_operator(dict, Picture::Schrodinger); BOOST_CHECK(op.init_op_map.find(fresh) != op.init_op_map.end()); } @@ -216,7 +216,7 @@ BOOST_AUTO_TEST_CASE(mp_operator_update_initial_operator_majorana_encode_identit OperatorDict dict; dict[VecZ{}] = cd(2.75, 0.0); - op.update_initial_operator(dict, /*schrodinger=*/false); + op.update_initial_operator(dict, Picture::Heisenberg); BOOST_REQUIRE_EQUAL(op.op_coeffs.size(), 1U); BOOST_CHECK_EQUAL(op.op_coeffs[0], algebra_encode_coeff<8>(Basis::Majorana, cd(2.75, 0.0), identity)); BOOST_CHECK_EQUAL(op.op_coeffs[0], 2.75); diff --git a/cpp/tests/mpfunctions.cpp b/cpp/tests/mpfunctions.cpp index b4ac2727..67249b43 100644 --- a/cpp/tests/mpfunctions.cpp +++ b/cpp/tests/mpfunctions.cpp @@ -264,9 +264,9 @@ BOOST_AUTO_TEST_CASE(sparse_energy_matches_the_dense_gradient_value_bit_exactly) constexpr size_t kNumModes = 8; const auto data = test_utils::load_case_data("random_exact.msgpack"); - for (const auto schrodinger_cutoff : {std::optional{}, std::optional{4}}) { - BOOST_TEST_CONTEXT("schrodinger_cutoff = " << (schrodinger_cutoff ? "4" : "none")) { - test_utils::SimulatorConfig cfg{.schrodinger_cutoff = schrodinger_cutoff, .comm = MPI_COMM_SELF}; + for (const PictureSpec picture : {PictureSpec{Heisenberg{}}, PictureSpec{Schrodinger{4U}}}) { + BOOST_TEST_CONTEXT("picture = " << test_utils::picture_label(picture)) { + test_utils::SimulatorConfig cfg{.picture = picture, .comm = MPI_COMM_SELF}; auto sim = test_utils::build_simulator(data, cfg); sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); @@ -291,7 +291,7 @@ BOOST_AUTO_TEST_CASE(interleaved_gradients_do_not_share_scratch_state) { auto build = [&data](unsigned int cutoff) { auto sim = - MonomialPropagator(data.hamiltonian, cutoff, data.initial_state, std::nullopt, MPI_COMM_SELF); + MonomialPropagator(data.hamiltonian, cutoff, data.initial_state, Heisenberg{}, MPI_COMM_SELF); sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); return sim; }; diff --git a/cpp/tests/mpi_distributed_layer_equivalence.cpp b/cpp/tests/mpi_distributed_layer_equivalence.cpp index 152eb2b7..6a71b49d 100644 --- a/cpp/tests/mpi_distributed_layer_equivalence.cpp +++ b/cpp/tests/mpi_distributed_layer_equivalence.cpp @@ -50,7 +50,7 @@ auto run_energy(const TestInputs& inputs, MPI_Comm comm) -> double { MonomialPropagator sim(inputs.data.hamiltonian, kCutoff, inputs.data.initial_state, - std::nullopt, + Heisenberg{}, comm, std::nullopt, std::nullopt, @@ -84,7 +84,7 @@ BOOST_AUTO_TEST_CASE(gradient_rank_count_within_fp_tolerance) { MonomialPropagator sim(inputs.data.hamiltonian, kCutoff, inputs.data.initial_state, - std::nullopt, + Heisenberg{}, comm, std::nullopt, std::nullopt, @@ -118,7 +118,7 @@ auto run_pauli_energy(MPI_Comm comm) -> double { MonomialPropagator sim(init, kPauliQ, VecZ{}, - std::nullopt, + Heisenberg{}, comm, 1e-12, std::nullopt, @@ -169,7 +169,7 @@ auto run_energy_partitioned(const TestInputs& inputs, MPI_Comm comm, size_t part MonomialPropagator sim(inputs.data.hamiltonian, kCutoff, inputs.data.initial_state, - std::nullopt, + Heisenberg{}, comm, std::nullopt, std::nullopt, diff --git a/cpp/tests/mpi_fresh_insert_equivalence.cpp b/cpp/tests/mpi_fresh_insert_equivalence.cpp index 2973131b..42543ca9 100644 --- a/cpp/tests/mpi_fresh_insert_equivalence.cpp +++ b/cpp/tests/mpi_fresh_insert_equivalence.cpp @@ -43,7 +43,7 @@ auto run_schrodinger_majorana(const CaseData& data, MPI_Comm comm) -> double { MonomialPropagator sim(data.hamiltonian, /*cutoff=*/2U, data.initial_state, - /*schrodinger_cutoff=*/std::optional{4U}, + /*picture=*/Schrodinger{4U}, comm, /*lower_atol=*/std::nullopt, /*upper_atol=*/std::optional{0.0}, @@ -76,7 +76,7 @@ auto run_schrodinger_pauli(MPI_Comm comm) -> double { MonomialPropagator sim(init, /*cutoff=*/2U, VecZ{}, - /*schrodinger_cutoff=*/std::optional{4U}, + /*picture=*/Schrodinger{4U}, comm, /*lower_atol=*/std::nullopt, /*upper_atol=*/std::optional{0.0}, diff --git a/cpp/tests/pare_graph_tests.cpp b/cpp/tests/pare_graph_tests.cpp index 9e3464d2..369f4b5c 100644 --- a/cpp/tests/pare_graph_tests.cpp +++ b/cpp/tests/pare_graph_tests.cpp @@ -105,7 +105,7 @@ BOOST_AUTO_TEST_CASE(pare_graph_emits_expected_layer_kinds) { seed.push_back(i); } - auto pared = pare_graph(graph, seed, local_index_count, /*schrodinger=*/false, MPI_COMM_SELF, provider); + auto pared = pare_graph(graph, seed, local_index_count, Picture::Heisenberg, MPI_COMM_SELF, provider); BOOST_REQUIRE_EQUAL(pared.layers(), graph.layers()); size_t pruned_count = 0; diff --git a/cpp/tests/partition_equivalence_tests.cpp b/cpp/tests/partition_equivalence_tests.cpp index 9c8a984e..72c0c937 100644 --- a/cpp/tests/partition_equivalence_tests.cpp +++ b/cpp/tests/partition_equivalence_tests.cpp @@ -45,7 +45,7 @@ auto majorana_sim(const CaseData &data, size_t partitions, std::optional return MonomialPropagator(data.hamiltonian, kCutoff, data.initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, lower_atol, std::nullopt, @@ -177,7 +177,7 @@ BOOST_AUTO_TEST_CASE(partition_setters_reach_every_partition) { MonomialPropagator sim(data.hamiltonian, cutoff, data.initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, @@ -223,7 +223,7 @@ auto pauli_sim(const std::map &obs, size_t partitions) -> M return MonomialPropagator(init, /*cutoff=*/kNq, /*initial_state=*/{}, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, /*lower_atol=*/1e-12, std::nullopt, @@ -286,7 +286,7 @@ BOOST_AUTO_TEST_CASE(partition_factory_exception_propagates_without_terminate) { BOOST_CHECK_THROW(MonomialPropagator(data.hamiltonian, kCutoff, data.initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, @@ -303,7 +303,7 @@ BOOST_AUTO_TEST_CASE(partition_factory_exception_propagates_without_terminate) { BOOST_CHECK_THROW(MonomialPropagator(bad_op, kCutoff, data.initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, diff --git a/cpp/tests/pauli_build_layer_tests.cpp b/cpp/tests/pauli_build_layer_tests.cpp index af069fcd..060efe1e 100644 --- a/cpp/tests/pauli_build_layer_tests.cpp +++ b/cpp/tests/pauli_build_layer_tests.cpp @@ -62,7 +62,7 @@ auto jw_basis_indices(size_t n) -> std::vector { template auto build_pauli_sim(const std::map &obs, unsigned int cutoff, - std::optional schrodinger_cutoff = std::nullopt, + PictureSpec picture = Heisenberg{}, const VecZ &initial_state = {}, std::optional lower_atol = std::nullopt) -> MonomialPropagator { OperatorDict init; @@ -72,7 +72,7 @@ auto build_pauli_sim(const std::map &obs, return MonomialPropagator(init, cutoff, initial_state, - schrodinger_cutoff, + picture, MPI_COMM_SELF, lower_atol, std::nullopt, @@ -232,7 +232,7 @@ auto jw_gate_arrays(const PauliCircuit &c) -> std::pair, VecD> template auto build_jw_sim(const std::map &obs, unsigned int cutoff, - std::optional schrodinger_cutoff = std::nullopt, + PictureSpec picture = Heisenberg{}, const VecZ &initial_state = {}, std::optional lower_atol = std::nullopt) -> MonomialPropagator { OperatorDict init; @@ -243,7 +243,7 @@ auto build_jw_sim(const std::map &obs, return MonomialPropagator(init, cutoff, initial_state, - schrodinger_cutoff, + picture, MPI_COMM_SELF, lower_atol, std::nullopt, @@ -342,21 +342,21 @@ BOOST_AUTO_TEST_CASE(pauli_build_layer_jw_isomorphism) { const auto [jw_majs, jw_gcs] = jw_gate_arrays(circ); struct Cfg { - std::optional sch; + PictureSpec picture; unsigned int cutoff; std::optional atol; const char *name; }; const std::vector cfgs{ - {std::nullopt, 3, std::nullopt, "heisenberg-full-cutoff"}, - {std::nullopt, 2, std::nullopt, "heisenberg-cutoff-2"}, - {std::nullopt, 3, std::optional(1e-6), "heisenberg-lower-atol"}, - {std::optional(5), 3, std::nullopt, "schrodinger-full-cutoff"}, - {std::optional(5), 3, std::optional(1e-6), "schrodinger-lower-atol"}, + {Heisenberg{}, 3, std::nullopt, "heisenberg-full-cutoff"}, + {Heisenberg{}, 2, std::nullopt, "heisenberg-cutoff-2"}, + {Heisenberg{}, 3, std::optional(1e-6), "heisenberg-lower-atol"}, + {Schrodinger{5}, 3, std::nullopt, "schrodinger-full-cutoff"}, + {Schrodinger{5}, 3, std::optional(1e-6), "schrodinger-lower-atol"}, }; for (const auto &cf : cfgs) { - auto nat = build_pauli_sim(obs, cf.cutoff, cf.sch, initial_state, cf.atol); - auto jw = build_jw_sim(obs, cf.cutoff, cf.sch, initial_state, cf.atol); + auto nat = build_pauli_sim(obs, cf.cutoff, cf.picture, initial_state, cf.atol); + auto jw = build_jw_sim(obs, cf.cutoff, cf.picture, initial_state, cf.atol); nat.build_graph(nat_monos, circ.param_map, nat_gcs); jw.build_graph(jw_majs, circ.param_map, jw_gcs); const double en = nat.expectation_value(circ.params); @@ -423,17 +423,17 @@ BOOST_AUTO_TEST_CASE(pauli_build_layer_replay_fold_consumers) { const auto [jw_majs, jw_gcs] = jw_gate_arrays(circ); // (a) fused contract-immediately propagate. - auto prop = build_pauli_sim(obs, 3, std::nullopt, initial_state); + auto prop = build_pauli_sim(obs, 3, Heisenberg{}, initial_state); prop.propagate(nat_monos, circ.param_map, nat_gcs, circ.params); const double e_prop = heisenberg_expval(prop); // (b) graph build + functional (expectation_value recomputes the cos from the fold). - auto grp = build_pauli_sim(obs, 3, std::nullopt, initial_state); + auto grp = build_pauli_sim(obs, 3, Heisenberg{}, initial_state); grp.build_graph(nat_monos, circ.param_map, nat_gcs); const double e_graph = grp.expectation_value(circ.params); // (c) contract_partially (evolve_operator_with_recompute — the same fold path, non-inplace). - auto ctr = build_pauli_sim(obs, 3, std::nullopt, initial_state); + auto ctr = build_pauli_sim(obs, 3, Heisenberg{}, initial_state); ctr.build_graph(nat_monos, circ.param_map, nat_gcs); const auto evolved = ctr.contract_partially(circ.params, /*inplace=*/false); const VecD st = ctr.mp_op().materialize_state(); @@ -444,7 +444,7 @@ BOOST_AUTO_TEST_CASE(pauli_build_layer_replay_fold_consumers) { const double e_contract = ctr.core_term() + s; // (d) JW-image Majorana reference. - auto jw = build_jw_sim(obs, 3, std::nullopt, initial_state); + auto jw = build_jw_sim(obs, 3, Heisenberg{}, initial_state); jw.build_graph(jw_majs, circ.param_map, jw_gcs); const double e_jw = jw.expectation_value(circ.params); diff --git a/cpp/tests/update_initial_operator.cpp b/cpp/tests/update_initial_operator.cpp index 76dc74ac..6c61b0fb 100644 --- a/cpp/tests/update_initial_operator.cpp +++ b/cpp/tests/update_initial_operator.cpp @@ -35,7 +35,7 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_updates_core_expval) { MonomialPropagator simulator(initial_ham, 2 * n_modes, initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, @@ -67,7 +67,7 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_invalidates_gradient_functional) { MonomialPropagator simulator(initial_ham, 2 * n_modes, initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, @@ -96,7 +96,7 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_throws_for_unknown_term_in_heisenbe MonomialPropagator simulator(initial_ham, 2 * n_modes, initial_state, - std::nullopt, + Heisenberg{}, MPI_COMM_SELF, std::nullopt, std::nullopt, @@ -121,7 +121,7 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_accepts_new_terms_in_schrodinger) { MonomialPropagator simulator(initial_ham, cutoff, initial_state, - cutoff, + Schrodinger{cutoff}, MPI_COMM_SELF, std::nullopt, std::nullopt, diff --git a/docs/content/docs/features/initialisation.mdx b/docs/content/docs/features/initialisation.mdx index 95dc79ab..c13a499f 100644 --- a/docs/content/docs/features/initialisation.mdx +++ b/docs/content/docs/features/initialisation.mdx @@ -16,6 +16,10 @@ Pass `schrodinger_cutoff` to switch from the default Heisenberg picture Schrödinger picture gates are applied to the state, whereas in the Heisenberg picture they are applied to the observable. +In C++ the choice is a `PictureSpec` — `Heisenberg{}` or `Schrodinger{state_cutoff}` — so only a +Schrödinger run carries a state cutoff. The Python constructors keep the `schrodinger_cutoff` +spelling and resolve it at the binding boundary. + `schrodinger_cutoff` is the truncation applied to the initial state. It must be set higher than the regular `cutoff` for the Schrödinger-picture result to match the Heisenberg-picture one. A good value of the parameter depends on the circuit diff --git a/docs/content/docs/testing.mdx b/docs/content/docs/testing.mdx index 2485508c..a97b76bf 100644 --- a/docs/content/docs/testing.mdx +++ b/docs/content/docs/testing.mdx @@ -196,11 +196,8 @@ BOOST_DATA_TEST_CASE_F(ExampleDataFix, bdata::make(ds_pare_values) ^ bdata::make(ds_schrodinger_enabled), pare, sch_enabled) { - const auto schrodinger_cutoff = make_schrodinger_cutoff(sch_enabled, cutoff); SimulatorConfig cfg{ - .schrodinger_cutoff = schrodinger_cutoff - ? std::optional(*schrodinger_cutoff) - : std::nullopt, + .picture = make_picture(sch_enabled, cutoff), .cutoff_type = cutoff_type, .basis_change = basis_change, }; diff --git a/src/monoprop/bindings/binder.h b/src/monoprop/bindings/binder.h index 5728aca5..3414a9bb 100644 --- a/src/monoprop/bindings/binder.h +++ b/src/monoprop/bindings/binder.h @@ -72,10 +72,14 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { size_t logical_num_modes, const std::string &basis, size_t partitions) { + // Python keeps the historical "a cutoff selects Schrodinger" spelling; the engine wants the + // choice as a variant, so it is resolved here, at the boundary. + const PictureSpec picture = schrodinger_cutoff.has_value() ? PictureSpec{Schrodinger{*schrodinger_cutoff}} + : PictureSpec{Heisenberg{}}; new (t) MonomialPropagator(initial_operator, cutoff, initial_state, - schrodinger_cutoff, + picture, get_mpi_comm(py_comm), lower_atol, upper_atol, From 1e0a160108392a2ac61a0229f255c507f5707ec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 08:32:42 +0000 Subject: [PATCH 02/13] docs: :memo: describe where the picture policy is actually bound The Picture bullet still said cold sites reach a policy through `picture_*()` helpers. Those helpers were deleted when every public entry point started binding the policy once with `with_picture`, so the sentence described a layer that no longer exists. Record the rule that replaced it, and why `build_layer` binds at the sink rather than templating itself. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 708274b5..6194541d 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -71,10 +71,13 @@ Key files: simulation pictures are sibling models (`HeisenbergPicture`, `SchrodingerPicture`), built the same way as the algebras. Each states one picture's whole rule set — gate traversal direction, applied-angle sign, `map_params` phase, the live coefficient vector, the contraction partner — so no `if (schrodinger)` is - written twice. Cold sites reach a policy through the `picture_*()` helpers; the fused - `ContractSink`/`apply_fused_contract` pair takes the policy as a template parameter, bound once inside - `build_layer`. The picture is fixed at construction: the constructor takes a - `PictureSpec = std::variant`, so only a Schrodinger run carries a state cutoff. + written twice. Every public entry point of `MonomialPropagator` binds the policy once with `with_picture`; + its whole private layer is templated on that policy and never re-tests which picture it is in, so there is + no runtime-dispatching helper layer. The fused `ContractSink`/`apply_fused_contract` pair likewise takes + the policy as a template parameter, bound once inside `build_layer` — at the sink only, or the + `with_algebra` scan above it would instantiate four times per mode width instead of two. The picture is + fixed at construction: the constructor takes a `PictureSpec = std::variant`, so + only a Schrodinger run carries a state cutoff. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` From 0bcadd9df79a814c3e625ee310fb3e00d87894c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 08:52:44 +0000 Subject: [PATCH 03/13] refactor(c++)!: :art: give `MPGraph` a layer order, not a picture `MPGraph` held a `Picture` and asked `is_schrodinger()` in three places, so an exported class that only ever needed to know which end of its layer vector a new gate attaches to carried a simulation concept instead. Replace it with a graph-local `LayerGrowth { Back, Front }`. `core/Picture.h` leaves the exported header, `picture()` and `is_schrodinger()` go (they had no callers outside `MPGraph.cpp` and one `pare_graph` line), and the three ordering-sensitive members ask the question in one spelling, `grows_at_front()`. `layer_growth_of` in `picture/Picture.h` is the single translation point, in the direction that keeps the dependency one-way. `pare_graph` loses its `Picture` parameter outright: it already had the graph, so it now reads the sweep direction off `graph.growth()` and cannot be handed a direction that disagrees with the graph it is paring. Write down the invariant that made this bounded, because it is easy to assume the opposite. Layers are stored in DESCENDING optimizer-slot order under both growth ends, and that mapping is depended on outside the class -- MPFunctions' `prepare_evolved_operator` hard-codes `fill_mapped_params(..., reverse=true)`, its gradient loop hard-codes `count-1-i`, and `graph_gate_arrays_` hard-codes the same. All three are picture-free only because `append_layer` normalizes two arrival orders into one storage order. `slice_graph` is the documented exception: its result is ordered for replay, which coincides with that mapping only for `LayerGrowth::Back`. That invariant is also why the direction bit itself cannot be removed. Storing layers by absolute optimizer slot is impossible -- `propagate()` appends one gate at a time, and in the back-growth order each new gate takes slot 0 and shifts every existing layer up -- and storing in arrival order instead would push the picture into `ev`/`ev_and_grad`, which take none today. Results are bit-identical over the same 95 expectation-value and gradient fingerprints as the picture refactor. BREAKING CHANGE: `MPGraph`'s constructors take `LayerGrowth` where they took `Picture`, `picture()` and `is_schrodinger()` are gone, and `pare_graph` no longer takes a `Picture`. The Python API is unchanged. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 5 ++- cpp/include/monoprop/MPFunctions.h | 2 +- cpp/include/monoprop/MPGraph.h | 41 +++++++++++++++---- cpp/monoprop/MPGraph.cpp | 8 ++-- .../MonomialPropagator.inl | 4 +- cpp/monoprop/detail/pare/PareGraph.cpp | 7 ++-- cpp/monoprop/picture/Picture.h | 7 ++++ cpp/tests/GraphBuildHarness.h | 6 +-- cpp/tests/mp_graph_tests.cpp | 32 ++++++++------- cpp/tests/pare_graph_tests.cpp | 2 +- 10 files changed, 76 insertions(+), 38 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 6194541d..5419f0e7 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -77,7 +77,10 @@ Key files: the policy as a template parameter, bound once inside `build_layer` — at the sink only, or the `with_algebra` scan above it would instantiate four times per mode width instead of two. The picture is fixed at construction: the constructor takes a `PictureSpec = std::variant`, so - only a Schrodinger run carries a state cutoff. + only a Schrodinger run carries a state cutoff. `MPGraph` deliberately knows no picture — it takes a + graph-local `LayerGrowth` saying which end a new layer attaches to, and `layer_growth_of` in + `picture/Picture.h` is the only translation. The graph's layer order itself is picture-independent and + load-bearing beyond the class; `MPGraph`'s own comment records which three call sites depend on it. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index d91a2097..0d4d9f7c 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -107,10 +107,10 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, const detail::CosCallbacks &cos = {}) -> std::pair; /// Prune `graph` to the subgraph reaching `nonzero_inds`; `full_cos_of_layer(i)` supplies layer i's full cosine set. +// The sweep direction comes from graph.growth(); there is no picture argument to disagree with it. monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, - Picture picture, mpi::Comm comm, const std::function &full_cos_of_layer) -> MPGraph; } // namespace monoprop diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 2f859447..b22d3691 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -15,26 +15,46 @@ #pragma once #include +#include #include #include #include #include #include -#include "monoprop/core/Picture.h" #include "monoprop/detail/graph/MPGraphLayers.h" #include "monoprop/detail/graph/MPGraphViews.h" #include "monoprop/monopropExport.h" namespace monoprop { +/// Which end of the layer store a newly appended gate attaches to. +// The graph needs this one bit and nothing more about the simulation: it does not name a picture, because +// the layer order it maintains is the same in both (see MPGraph). The two pictures assign optimizer slots +// to arriving gates in opposite directions, and that is all this distinguishes. +enum class LayerGrowth : uint8_t { + Back, ///< a new gate takes the lowest optimizer slot, so it attaches at the back + Front, ///< a new gate takes the highest optimizer slot, so it attaches at the front +}; + /// Ordered per-rank record of the evolution circuit, one Layer per generator. +// A graph built by append() stores its layers in DESCENDING optimizer-slot order under both growth ends: +// active layer i is optimizer slot layers()-1-i. That mapping is load-bearing outside this class, not just +// a local convention -- MPFunctions' prepare_evolved_operator hard-codes fill_mapped_params(..., reverse= +// true), its gradient loop hard-codes `count-1-i`, and MonomialPropagator::graph_gate_arrays_ hard-codes +// the same mapping. All three can stay free of any picture only because append_layer() normalizes the two +// arrival orders into this one storage order, so changing the storage order means changing all three with +// it. +// +// slice_graph() is the one exception: its result is ordered for replay, which coincides with the mapping +// above only for LayerGrowth::Back. Slices are for replay_view() and must not be fed to anything that +// reconstructs optimizer order. class monoprop_EXPORT MPGraph { private: using LayerIterator = std::vector::iterator; using ConstLayerIterator = std::vector::const_iterator; - Picture picture_; + LayerGrowth growth_; std::vector layers_; size_t front_offset_ = 0; @@ -54,9 +74,12 @@ class monoprop_EXPORT MPGraph { auto active_end_iterator() const -> ConstLayerIterator { return layers_.end(); } - // Schrödinger prepends: it consumes the circuit front-to-back, so the newest gate is the earliest layer. + // The growth end and the oldest end are opposite: the first gate appended sits at whichever end new + // layers do not take. Every ordering-sensitive member below asks the question in this one spelling. + auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } + auto append_position() -> LayerIterator { - return is_schrodinger() ? active_begin_iterator() : active_end_iterator(); + return grows_at_front() ? active_begin_iterator() : active_end_iterator(); } auto append_layer(Layer layer) -> void { layers_.emplace(append_position(), std::move(layer)); } @@ -69,9 +92,9 @@ class monoprop_EXPORT MPGraph { } public: - explicit MPGraph(Picture picture) : picture_(picture) {} + explicit MPGraph(LayerGrowth growth) : growth_(growth) {} - explicit MPGraph(Picture picture, std::vector layers) : picture_(picture), layers_(std::move(layers)) {} + explicit MPGraph(LayerGrowth growth, std::vector layers) : growth_(growth), layers_(std::move(layers)) {} /// Gate info (param_index, gen_coeff, gate_index) is written onto `storage` here while it is still /// mutable, before it is frozen into the Layer's shared const core. @@ -87,6 +110,7 @@ class monoprop_EXPORT MPGraph { /// Slice the graph at `key` (the number of earliest operations to include); `contract` also removes /// the sliced part from this graph. + // Layers come back in application order, earliest first -- see the note on this class. auto slice_graph(size_t key, bool contract = false) -> MPGraph; auto slice_view(size_t key) const -> MPGraphView; @@ -102,9 +126,8 @@ class monoprop_EXPORT MPGraph { /// Non-owning replay view over the active layers, in build order. auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } - auto picture() const -> Picture { return picture_; } - - auto is_schrodinger() const -> bool { return picture_ == Picture::Schrodinger; } + /// Carried so a derived graph (a slice, a pared copy) keeps its source's layer order. + auto growth() const -> LayerGrowth { return growth_; } /// A normally-built layer stores no cosine set, so the companion cosine-index count cannot come from /// the graph: only the operator's inverted index can supply it. diff --git a/cpp/monoprop/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index 2fd4237f..79c9f116 100644 --- a/cpp/monoprop/MPGraph.cpp +++ b/cpp/monoprop/MPGraph.cpp @@ -62,7 +62,7 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { const auto k = std::min(key, layers()); sliced_layers.reserve(k); - if (is_schrodinger()) { + if (grows_at_front()) { const size_t active_end = active_end_index(); for (size_t i = 0; i < k; ++i) { sliced_layers.push_back(layers_[active_end - 1 - i]); @@ -84,12 +84,14 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { } } - return MPGraph(picture_, std::move(sliced_layers)); + return MPGraph(growth_, std::move(sliced_layers)); } auto MPGraph::slice_view(size_t key) const -> MPGraphView { const auto k = std::min(key, layers()); - if (is_schrodinger()) { + // Oldest end at the back: the window is the tail, walked back-to-front so the view still yields the + // earliest operation first. + if (grows_at_front()) { return {layers_, active_end_index() - k, k, true}; } return {layers_, active_begin_index(), k, false}; diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 00eba12b..38c85a76 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -86,7 +86,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope : picture_{kind_of(picture)}, comm_{comm}, mp_op_{}, - graph_(picture_), + graph_(layer_growth_of(picture_)), cutoff_{cutoff}, lower_atol_{lower_atol}, upper_atol_{upper_atol}, @@ -965,7 +965,7 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optional(pare_graph(graph_, keep, count, P::picture, comm_, full_cos_of_layer)); + graph = std::make_shared(pare_graph(graph_, keep, count, comm_, full_cos_of_layer)); } else { graph = std::shared_ptr(std::shared_ptr{}, &graph_); diff --git a/cpp/monoprop/detail/pare/PareGraph.cpp b/cpp/monoprop/detail/pare/PareGraph.cpp index 8ca382c5..458f49e2 100644 --- a/cpp/monoprop/detail/pare/PareGraph.cpp +++ b/cpp/monoprop/detail/pare/PareGraph.cpp @@ -103,7 +103,6 @@ auto mark_cross_rank_endpoints_kept(const LayerTraversal &layer, size_t my_rank, auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, - Picture picture, mpi::Comm comm, const std::function &full_cos_of_layer) -> MPGraph { const size_t num_layers = graph.layers(); @@ -122,7 +121,9 @@ auto pare_graph(const MPGraph &graph, // mark_cross_rank_endpoints_kept), so nodes_to_keep stays consistent across ranks with no exchange. // Cross-rank lists are never pruned; the keep-set only has to be right so cos pruning stays exact. for (size_t iter = 0; iter < num_layers; ++iter) { - const size_t layer_idx = picture == Picture::Schrodinger ? iter : (num_layers - 1 - iter); + // Backward means last-applied first, which is the growth end: new layers attach where the + // circuit's latest operations are. + const size_t layer_idx = graph.growth() == LayerGrowth::Front ? iter : (num_layers - 1 - iter); const auto &layer = graph.get_layer(layer_idx); const auto lt = layer.traversal(); @@ -135,7 +136,7 @@ auto pare_graph(const MPGraph &graph, layers[layer_idx] = preserves ? Layer(layer.shared_core()) : Layer(layer.shared_core(), std::move(filtered)); } - return MPGraph(graph.picture(), std::move(layers)); + return MPGraph(graph.growth(), std::move(layers)); } } // namespace monoprop diff --git a/cpp/monoprop/picture/Picture.h b/cpp/monoprop/picture/Picture.h index 75758b1d..501139d8 100644 --- a/cpp/monoprop/picture/Picture.h +++ b/cpp/monoprop/picture/Picture.h @@ -152,6 +152,13 @@ concept PicturePolicy = requires { static_assert(PicturePolicy); static_assert(PicturePolicy); +// The only place a picture becomes a graph layer order. MPGraph deliberately knows nothing about pictures, +// and the translation belongs on this side: Heisenberg gives each arriving gate a lower optimizer slot than +// the last, Schrödinger a higher one. +inline auto layer_growth_of(Picture picture) -> LayerGrowth { + return picture == Picture::Schrodinger ? LayerGrowth::Front : LayerGrowth::Back; +} + // The one runtime->policy branch, taken once per public call. decltype(auto), not auto, so a policy that // hands back a reference into the operator does not decay to a copy; both arms must then deduce the same // type. diff --git a/cpp/tests/GraphBuildHarness.h b/cpp/tests/GraphBuildHarness.h index 5a862dee..a45a7998 100644 --- a/cpp/tests/GraphBuildHarness.h +++ b/cpp/tests/GraphBuildHarness.h @@ -34,9 +34,9 @@ inline auto layer_with_gate(std::size_t gate_index) -> monoprop::Layer { return monoprop::Layer(core_with_gate(gate_index)); } -// Built via append() so layer ordering (Heisenberg back-append, Schrödinger front-insert) is production's. -inline auto graph_with_gates(monoprop::Picture picture, std::size_t n) -> monoprop::MPGraph { - monoprop::MPGraph graph(picture); +// Built via append() so the layer ordering is production's for the growth end under test. +inline auto graph_with_gates(monoprop::LayerGrowth growth, std::size_t n) -> monoprop::MPGraph { + monoprop::MPGraph graph(growth); for (std::size_t i = 0; i < n; ++i) { graph.append(std::make_shared(), /*param_index=*/0, /*gen_coeff=*/0.0, /*gate_index=*/i); } diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index ee9cb422..669ee4cd 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -28,8 +28,9 @@ using test_utils::core_with_gate; using test_utils::graph_with_gates; using test_utils::layer_with_gate; -BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_heisenberg_prefix_no_contract) { - auto graph = graph_with_gates(Picture::Heisenberg, 5); // layers_ = [0,1,2,3,4] +// LayerGrowth::Back is the Heisenberg order: gate 0 arrives first and stays at the front. +BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_growth_back_prefix_no_contract) { + auto graph = graph_with_gates(LayerGrowth::Back, 5); // layers_ = [0,1,2,3,4] auto sliced = graph.slice_graph(3, /*contract=*/false); BOOST_REQUIRE_EQUAL(sliced.layers(), 3U); @@ -41,9 +42,10 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_heisenberg_prefix_no_contract) { BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 0U); } -BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_schrodinger_contract_newest_first_copy_and_resize) { - // Schrödinger stores newest-first: appending gates 0..4 gives layers_ = [4,3,2,1,0]. - auto graph = graph_with_gates(Picture::Schrodinger, 5); +// LayerGrowth::Front is the Schrödinger order. +BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_growth_front_contract_newest_first_copy_and_resize) { + // Front growth stores newest-first: appending gates 0..4 gives layers_ = [4,3,2,1,0]. + auto graph = graph_with_gates(LayerGrowth::Front, 5); auto sliced = graph.slice_graph(2, /*contract=*/true); // sliced = layers_[active_end-1-i] = layers_[4], layers_[3] = gates 0, 1 (oldest-first). @@ -59,15 +61,15 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_schrodinger_contract_newest_first_copy } BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_key_clamped_to_size) { - auto graph = graph_with_gates(Picture::Heisenberg, 3); + auto graph = graph_with_gates(LayerGrowth::Back, 3); auto sliced = graph.slice_graph(100, /*contract=*/false); BOOST_CHECK_EQUAL(sliced.layers(), 3U); } -// The maybe_compact_layers arms below are reached through Heisenberg slice_graph(contract=true). +// The maybe_compact_layers arms below are reached through back-growth slice_graph(contract=true). BOOST_AUTO_TEST_CASE(mp_graph_contract_clear_arm_when_prefix_covers_all) { - auto graph = graph_with_gates(Picture::Heisenberg, 5); + auto graph = graph_with_gates(LayerGrowth::Back, 5); (void)graph.slice_graph(5, /*contract=*/true); // front_offset == size -> clear BOOST_CHECK_EQUAL(graph.layers(), 0U); // Graph is still usable after a full clear. @@ -77,7 +79,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_clear_arm_when_prefix_covers_all) { } BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { - auto graph = graph_with_gates(Picture::Heisenberg, 100); + auto graph = graph_with_gates(LayerGrowth::Back, 100); (void)graph.slice_graph(3, /*contract=*/true); // front_offset 3 < 4096 -> no physical compaction BOOST_REQUIRE_EQUAL(graph.layers(), 97U); BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 3U); @@ -86,7 +88,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { BOOST_AUTO_TEST_CASE(mp_graph_contract_erase_arm_above_threshold) { // The erase arm fires only when front_offset >= 4096 AND 2*front_offset >= size. - auto graph = graph_with_gates(Picture::Heisenberg, 8200); + auto graph = graph_with_gates(LayerGrowth::Back, 8200); auto sliced = graph.slice_graph(4100, /*contract=*/true); BOOST_CHECK_EQUAL(sliced.layers(), 4100U); BOOST_CHECK_EQUAL(sliced.get_layer_traversal(0).gate_index(), 0U); @@ -96,8 +98,8 @@ BOOST_AUTO_TEST_CASE(mp_graph_contract_erase_arm_above_threshold) { BOOST_CHECK_EQUAL(graph.get_layer_traversal(4099).gate_index(), 8199U); } -BOOST_AUTO_TEST_CASE(mp_graph_slice_view_heisenberg_forward_window) { - auto graph = graph_with_gates(Picture::Heisenberg, 5); +BOOST_AUTO_TEST_CASE(mp_graph_slice_view_growth_back_forward_window) { + auto graph = graph_with_gates(LayerGrowth::Back, 5); auto view = graph.slice_view(3); BOOST_REQUIRE_EQUAL(view.layers(), 3U); BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); @@ -105,10 +107,10 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_view_heisenberg_forward_window) { BOOST_CHECK_EQUAL(view.get_layer_traversal(2).gate_index(), 2U); } -BOOST_AUTO_TEST_CASE(mp_graph_slice_view_schrodinger_reversed_window) { +BOOST_AUTO_TEST_CASE(mp_graph_slice_view_growth_front_reversed_window) { // layers_ = [4,3,2,1,0]; slice_view(3) uses base=active_end-3=2, reverse=true. // get_layer_traversal(i) -> layers_[2 + (3-1-i)] -> gates 0,1,2 in replay order. - auto graph = graph_with_gates(Picture::Schrodinger, 5); + auto graph = graph_with_gates(LayerGrowth::Front, 5); auto view = graph.slice_view(3); BOOST_REQUIRE_EQUAL(view.layers(), 3U); BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); @@ -133,7 +135,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { } BOOST_AUTO_TEST_CASE(mp_graph_get_layer_out_of_range_throws) { - auto graph = graph_with_gates(Picture::Heisenberg, 3); + auto graph = graph_with_gates(LayerGrowth::Back, 3); BOOST_CHECK_NO_THROW((void)graph.get_layer(2)); BOOST_CHECK_THROW((void)graph.get_layer(3), std::out_of_range); // const overload takes the same guard. diff --git a/cpp/tests/pare_graph_tests.cpp b/cpp/tests/pare_graph_tests.cpp index 369f4b5c..2464b2d0 100644 --- a/cpp/tests/pare_graph_tests.cpp +++ b/cpp/tests/pare_graph_tests.cpp @@ -105,7 +105,7 @@ BOOST_AUTO_TEST_CASE(pare_graph_emits_expected_layer_kinds) { seed.push_back(i); } - auto pared = pare_graph(graph, seed, local_index_count, Picture::Heisenberg, MPI_COMM_SELF, provider); + auto pared = pare_graph(graph, seed, local_index_count, MPI_COMM_SELF, provider); BOOST_REQUIRE_EQUAL(pared.layers(), graph.layers()); size_t pruned_count = 0; From 82af677f65a6b75297245a3f92af94961fa8176b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 09:13:25 +0000 Subject: [PATCH 04/13] refactor(c++): :broom: keep the layer-growth bit in one place Cleanup pass over the LayerGrowth change. The growth end is `gate_slot`'s slope, so it belongs on the picture policies beside it rather than in a second free function that re-branches on `Picture`. `HeisenbergPicture`/`SchrodingerPicture` now carry `layer_growth`, the `PicturePolicy` concept requires it, and `layer_growth_of` is gone -- which makes `picture/Picture.h`'s own "no runtime-dispatching helper layer on purpose" true again. The one caller binds it through `with_picture` like every other policy read. `MPGraph::grows_at_front()` becomes public so `pare_graph` asks the growth question in the same spelling as the class's own members instead of comparing enumerators in another translation unit. `slice_view`'s two returns collapse into one now that the reverse flag is visibly the same predicate, and `slice_graph`'s hand-rolled reverse-index loop becomes the reverse-iterator range it was spelling out. Trim the documentation to what the code cannot say: drop the paraphrase of three remote call sites in favour of naming them, drop the enum note that restated its own enumerator docs, and drop a comment that described `pare_graph`'s removed parameter. Correct two comments this refactor falsified -- `core/Picture.h` named two exported headers that no longer carry a `Picture`, and `replay_view()` plus `graph_gate_arrays_` still framed storage as "build"/"simulation" order, which holds only for `LayerGrowth::Back`. Also add the `MPGraph.h` include that `picture/Picture.h` was getting transitively through `MPFunctions.h`. Still bit-identical over the same 95 fingerprints. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 6 ++-- cpp/include/monoprop/MPFunctions.h | 1 - cpp/include/monoprop/MPGraph.h | 35 ++++++++----------- cpp/monoprop/MPGraph.cpp | 19 +++++----- cpp/monoprop/core/Picture.h | 4 +-- .../MonomialPropagator.inl | 4 +-- cpp/monoprop/detail/pare/PareGraph.cpp | 2 +- cpp/monoprop/picture/Picture.h | 13 ++++--- cpp/tests/mp_graph_tests.cpp | 2 +- 9 files changed, 38 insertions(+), 48 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 5419f0e7..8f9caec2 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -78,9 +78,9 @@ Key files: `with_algebra` scan above it would instantiate four times per mode width instead of two. The picture is fixed at construction: the constructor takes a `PictureSpec = std::variant`, so only a Schrodinger run carries a state cutoff. `MPGraph` deliberately knows no picture — it takes a - graph-local `LayerGrowth` saying which end a new layer attaches to, and `layer_growth_of` in - `picture/Picture.h` is the only translation. The graph's layer order itself is picture-independent and - load-bearing beyond the class; `MPGraph`'s own comment records which three call sites depend on it. + graph-local `LayerGrowth` saying which end a new layer attaches to, and each policy carries its own + `layer_growth` beside `gate_slot`. The graph's layer order itself is picture-independent and load-bearing + beyond the class; `MPGraph`'s own comment records what depends on it. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index 0d4d9f7c..087fd0a8 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -107,7 +107,6 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, const detail::CosCallbacks &cos = {}) -> std::pair; /// Prune `graph` to the subgraph reaching `nonzero_inds`; `full_cos_of_layer(i)` supplies layer i's full cosine set. -// The sweep direction comes from graph.growth(); there is no picture argument to disagree with it. monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index b22d3691..552b09af 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -29,26 +29,18 @@ namespace monoprop { /// Which end of the layer store a newly appended gate attaches to. -// The graph needs this one bit and nothing more about the simulation: it does not name a picture, because -// the layer order it maintains is the same in both (see MPGraph). The two pictures assign optimizer slots -// to arriving gates in opposite directions, and that is all this distinguishes. +// This one bit is all the graph needs to know about the simulation that drives it. enum class LayerGrowth : uint8_t { Back, ///< a new gate takes the lowest optimizer slot, so it attaches at the back Front, ///< a new gate takes the highest optimizer slot, so it attaches at the front }; /// Ordered per-rank record of the evolution circuit, one Layer per generator. -// A graph built by append() stores its layers in DESCENDING optimizer-slot order under both growth ends: -// active layer i is optimizer slot layers()-1-i. That mapping is load-bearing outside this class, not just -// a local convention -- MPFunctions' prepare_evolved_operator hard-codes fill_mapped_params(..., reverse= -// true), its gradient loop hard-codes `count-1-i`, and MonomialPropagator::graph_gate_arrays_ hard-codes -// the same mapping. All three can stay free of any picture only because append_layer() normalizes the two -// arrival orders into this one storage order, so changing the storage order means changing all three with -// it. -// -// slice_graph() is the one exception: its result is ordered for replay, which coincides with the mapping -// above only for LayerGrowth::Back. Slices are for replay_view() and must not be fed to anything that -// reconstructs optimizer order. +// A graph built by append() stores its layers in DESCENDING optimizer-slot order under either growth end: +// active layer i is optimizer slot layers()-1-i. append_layer() normalizes the two arrival orders into that +// one storage order, which is what lets everything reconstructing optimizer order from a graph -- the +// evolved-operator setup and the gradient loop in MPFunctions, MonomialPropagator::graph_gate_arrays_ -- do +// it without knowing which growth end this graph has. Changing the storage order means changing those too. class monoprop_EXPORT MPGraph { private: using LayerIterator = std::vector::iterator; @@ -74,10 +66,6 @@ class monoprop_EXPORT MPGraph { auto active_end_iterator() const -> ConstLayerIterator { return layers_.end(); } - // The growth end and the oldest end are opposite: the first gate appended sits at whichever end new - // layers do not take. Every ordering-sensitive member below asks the question in this one spelling. - auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } - auto append_position() -> LayerIterator { return grows_at_front() ? active_begin_iterator() : active_end_iterator(); } @@ -110,7 +98,9 @@ class monoprop_EXPORT MPGraph { /// Slice the graph at `key` (the number of earliest operations to include); `contract` also removes /// the sliced part from this graph. - // Layers come back in application order, earliest first -- see the note on this class. + // Layers come back in application order, earliest first, which matches this class's storage order only + // for LayerGrowth::Back. A slice is for replay_view() and must not be fed to anything that + // reconstructs optimizer order from it. auto slice_graph(size_t key, bool contract = false) -> MPGraph; auto slice_view(size_t key) const -> MPGraphView; @@ -123,12 +113,17 @@ class monoprop_EXPORT MPGraph { auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } - /// Non-owning replay view over the active layers, in build order. + /// Non-owning replay view over the active layers, in stored (descending optimizer-slot) order. auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } /// Carried so a derived graph (a slice, a pared copy) keeps its source's layer order. auto growth() const -> LayerGrowth { return growth_; } + /// Whether new layers attach at the front, so the oldest operation is at the back. + // The one spelling of the growth question: every ordering-sensitive site, here and in pare_graph, asks + // it this way rather than comparing enumerators. + auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } + /// A normally-built layer stores no cosine set, so the companion cosine-index count cannot come from /// the graph: only the operator's inverted index can supply it. auto total_cycles() const -> size_t; diff --git a/cpp/monoprop/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index 79c9f116..dbee5bd4 100644 --- a/cpp/monoprop/MPGraph.cpp +++ b/cpp/monoprop/MPGraph.cpp @@ -63,13 +63,12 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { sliced_layers.reserve(k); if (grows_at_front()) { - const size_t active_end = active_end_index(); - for (size_t i = 0; i < k; ++i) { - sliced_layers.push_back(layers_[active_end - 1 - i]); - } + // Oldest end at the back, so the slice is the tail read backwards; active_end_index() is + // layers_.size(), which is what makes rbegin() the right starting point. + sliced_layers.insert(sliced_layers.end(), layers_.rbegin(), layers_.rbegin() + static_cast(k)); if (contract && k != 0) { - layers_.resize(active_end - k); + layers_.resize(active_end_index() - k); } } else { @@ -89,12 +88,10 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { auto MPGraph::slice_view(size_t key) const -> MPGraphView { const auto k = std::min(key, layers()); - // Oldest end at the back: the window is the tail, walked back-to-front so the view still yields the - // earliest operation first. - if (grows_at_front()) { - return {layers_, active_end_index() - k, k, true}; - } - return {layers_, active_begin_index(), k, false}; + // The window sits at the oldest end and is walked away from it, so the view yields the earliest + // operation first either way -- which makes the reverse flag exactly the growth predicate. + const bool front = grows_at_front(); + return {layers_, front ? active_end_index() - k : active_begin_index(), k, front}; } auto MPGraph::total_cycles() const -> size_t { diff --git a/cpp/monoprop/core/Picture.h b/cpp/monoprop/core/Picture.h index 654250a6..221bda74 100644 --- a/cpp/monoprop/core/Picture.h +++ b/cpp/monoprop/core/Picture.h @@ -14,8 +14,8 @@ #pragma once -// Deliberately dependency-free: exported headers (MPGraph.h, MPFunctions.h, MonomialPropagator.h) -// carry a Picture in their signatures. The behaviour behind each value lives in picture/Picture.h. +// Deliberately dependency-free: MonomialPropagator.h carries a Picture in its signatures, and MPOperator.h +// in a parameter. The behaviour behind each value lives in picture/Picture.h. #include #include diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 38c85a76..37e69eb4 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -86,7 +86,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope : picture_{kind_of(picture)}, comm_{comm}, mp_op_{}, - graph_(layer_growth_of(picture_)), + graph_(with_picture(picture_, []() { return P::layer_growth; })), cutoff_{cutoff}, lower_atol_{lower_atol}, upper_atol_{upper_atol}, @@ -871,7 +871,7 @@ auto MonomialPropagator::graph_gate_arrays_() const -> std::pair #include "monoprop/MPFunctions.h" +#include "monoprop/MPGraph.h" #include "monoprop/TypeAliases.h" #include "monoprop/algebra/AlgebraCommon.h" #include "monoprop/core/Monomial.h" @@ -45,6 +46,8 @@ struct HeisenbergPicture { // Simulation step i consumes optimizer slot n-1-i: the observable walks the circuit backwards. static auto gate_slot(size_t i, size_t n) -> size_t { return n - 1 - i; } + // gate_slot's slope, as MPGraph needs it: descending, so each arriving gate takes a lower slot. + static constexpr LayerGrowth layer_growth = LayerGrowth::Back; static constexpr double apply_sign = 1.0; // the applied angle is the build angle // map_params() arguments for contract_partially: forward phase, written in reverse. @@ -101,6 +104,8 @@ struct SchrodingerPicture { // Simulation step i consumes optimizer slot i: the state walks the circuit front-to-back. static auto gate_slot(size_t i, size_t /*n*/) -> size_t { return i; } + // gate_slot's slope, as MPGraph needs it: ascending, so each arriving gate takes a higher slot. + static constexpr LayerGrowth layer_growth = LayerGrowth::Front; static constexpr double apply_sign = -1.0; // the applied angle is the negated build angle static constexpr double contract_phase = -1.0; @@ -147,18 +152,12 @@ concept PicturePolicy = requires { { P::apply_sign } -> std::convertible_to; { P::contract_phase } -> std::convertible_to; { P::contract_reverse } -> std::convertible_to; + { P::layer_growth } -> std::convertible_to; }; static_assert(PicturePolicy); static_assert(PicturePolicy); -// The only place a picture becomes a graph layer order. MPGraph deliberately knows nothing about pictures, -// and the translation belongs on this side: Heisenberg gives each arriving gate a lower optimizer slot than -// the last, Schrödinger a higher one. -inline auto layer_growth_of(Picture picture) -> LayerGrowth { - return picture == Picture::Schrodinger ? LayerGrowth::Front : LayerGrowth::Back; -} - // The one runtime->policy branch, taken once per public call. decltype(auto), not auto, so a policy that // hands back a reference into the operator does not decay to a copy; both arms must then deduce the same // type. diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 669ee4cd..1dd99139 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -28,7 +28,7 @@ using test_utils::core_with_gate; using test_utils::graph_with_gates; using test_utils::layer_with_gate; -// LayerGrowth::Back is the Heisenberg order: gate 0 arrives first and stays at the front. +// LayerGrowth::Back is the Heisenberg order. BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_growth_back_prefix_no_contract) { auto graph = graph_with_gates(LayerGrowth::Back, 5); // layers_ = [0,1,2,3,4] auto sliced = graph.slice_graph(3, /*contract=*/false); From 8cce64920fd70855b341eada71ddd1ace7f43ff7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 10:11:18 +0000 Subject: [PATCH 05/13] refactor(c++): :broom: give MPGraph a contraction view and a drain contract_partially() is the only production caller of slice_graph()/slice_view(), and it always passes key == layers(): graph_gate_arrays_() sizes its arrays from graph_.layers(). The two functions therefore mean one view and one drain, so name them that way. contraction_view() needs no arms. active_end_index() - layers() equals active_begin_index(), so the window is the same under either growth end and only the reverse flag differs -- the key parameter was all that ever separated them. The inplace path also stops copying every Layer. It copied them out into an owned graph so the source could be drained before the evolution; with the order settled it can evolve over the graph's own layers and drain afterwards. slice_graph()/slice_view() stay for now, with their tests, so nothing loses coverage in this commit. Assisted-by: ClaudeCode:claude-opus-5 --- cpp/include/monoprop/MPGraph.h | 12 +++++ cpp/monoprop/MPGraph.cpp | 6 +++ .../MonomialPropagator.inl | 13 ++--- cpp/tests/mp_graph_tests.cpp | 47 +++++++++++++++++++ 4 files changed, 69 insertions(+), 9 deletions(-) diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 552b09af..bc65c9dd 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -96,6 +96,12 @@ class monoprop_EXPORT MPGraph { append_layer(Layer(std::move(storage))); } + /// Drop every layer. The graph stays usable, and a later append() starts from an empty store. + auto clear() -> void { + layers_.clear(); + front_offset_ = 0; + } + /// Slice the graph at `key` (the number of earliest operations to include); `contract` also removes /// the sliced part from this graph. // Layers come back in application order, earliest first, which matches this class's storage order only @@ -116,6 +122,12 @@ class monoprop_EXPORT MPGraph { /// Non-owning replay view over the active layers, in stored (descending optimizer-slot) order. auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } + /// The active layers in the order this graph's own build walked them, which is the order a + /// contraction must replay them in. + // Not replay_view(): a contraction drives the live coefficient vector, so it follows the simulation + // direction, where the picture-free evaluation order is always the stored one. + auto contraction_view() const -> MPGraphView; + /// Carried so a derived graph (a slice, a pared copy) keeps its source's layer order. auto growth() const -> LayerGrowth { return growth_; } diff --git a/cpp/monoprop/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index dbee5bd4..18140304 100644 --- a/cpp/monoprop/MPGraph.cpp +++ b/cpp/monoprop/MPGraph.cpp @@ -86,6 +86,12 @@ auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { return MPGraph(growth_, std::move(sliced_layers)); } +auto MPGraph::contraction_view() const -> MPGraphView { + // One window, two directions: active_end_index() - layers() is active_begin_index(), so only the + // reverse flag distinguishes the growth ends. + return {layers_, active_begin_index(), layers(), grows_at_front()}; +} + auto MPGraph::slice_view(size_t key) const -> MPGraphView { const auto k = std::min(key, layers()); // The window sits at the oldest end and is walked away from it, so the view yields the earliest diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 37e69eb4..58cc1b72 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -1077,23 +1077,18 @@ auto MonomialPropagator::contract_partially_(const VecD ¶meters, b return P::live_coeffs(mp_op_); } - const size_t num_majoranas = parameter_mapping.size(); // The pictures differ in three values only: the source vector, the (phase, reverse) map_params pair, // and the slot that receives an inplace result. Everything else -- and the order of every flop -- is shared. const VecD &source = P::live_coeffs(mp_op_); const auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, P::contract_phase, P::contract_reverse); - // Inplace slicing produces an owned MPGraph that must be bound to a named local before viewing - // (never view a temporary); slice_view() views this graph's still-live layers directly. - VecD evolved; + VecD evolved = evolve_operator_with_recompute_(VecD(source), graph_.contraction_view(), mapped_params); if (inplace) { - const MPGraph sliced = graph_.slice_graph(num_majoranas, true); - evolved = evolve_operator_with_recompute_(VecD(source), sliced.replay_view(), mapped_params); + // Drained only after the evolution: the cosine callbacks hold pointers into the layers' stored cos + // sets, and they die inside evolve_operator_with_recompute_. P::live_coeffs_slot(mp_op_) = evolved; - } - else { - evolved = evolve_operator_with_recompute_(VecD(source), graph_.slice_view(num_majoranas), mapped_params); + graph_.clear(); } return evolved; } diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 1dd99139..dfcedd51 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -118,6 +118,53 @@ BOOST_AUTO_TEST_CASE(mp_graph_slice_view_growth_front_reversed_window) { BOOST_CHECK_EQUAL(view.get_layer_traversal(2).gate_index(), 2U); } +// replay_view() is the evaluation order: stored order, so the growth end shows through it. +BOOST_AUTO_TEST_CASE(mp_graph_replay_view_is_stored_order) { + auto back = graph_with_gates(LayerGrowth::Back, 4); // layers_ = [0,1,2,3] + auto front = graph_with_gates(LayerGrowth::Front, 4); // layers_ = [3,2,1,0] + const auto back_view = back.replay_view(); + const auto front_view = front.replay_view(); + for (std::size_t i = 0; i < 4; ++i) { + BOOST_CHECK_EQUAL(back_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(front_view.get_layer_traversal(i).gate_index(), 3U - i); + } +} + +// contraction_view() is the build order, so it yields the same sequence under either growth end. That +// equality is the point: a contraction replays the circuit the way the build walked it. +BOOST_AUTO_TEST_CASE(mp_graph_contraction_view_is_build_order_under_either_growth) { + auto back = graph_with_gates(LayerGrowth::Back, 4); + auto front = graph_with_gates(LayerGrowth::Front, 4); + const auto back_view = back.contraction_view(); + const auto front_view = front.contraction_view(); + BOOST_REQUIRE_EQUAL(back_view.layers(), 4U); + BOOST_REQUIRE_EQUAL(front_view.layers(), 4U); + for (std::size_t i = 0; i < 4; ++i) { + BOOST_CHECK_EQUAL(back_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(front_view.get_layer_traversal(i).gate_index(), i); + } +} + +BOOST_AUTO_TEST_CASE(mp_graph_clear_empties_and_leaves_the_graph_usable) { + auto graph = graph_with_gates(LayerGrowth::Front, 5); + graph.clear(); + BOOST_REQUIRE_EQUAL(graph.layers(), 0U); + graph.append(std::make_shared(), 0, 0.0, /*gate_index=*/7); + BOOST_REQUIRE_EQUAL(graph.layers(), 1U); + BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 7U); +} + +// clear() must also drop the dead front prefix a contracting slice leaves behind. +BOOST_AUTO_TEST_CASE(mp_graph_clear_resets_the_dead_front_prefix) { + auto graph = graph_with_gates(LayerGrowth::Back, 100); + (void)graph.slice_graph(3, /*contract=*/true); // front_offset 3, below the compaction threshold + graph.clear(); + BOOST_REQUIRE_EQUAL(graph.layers(), 0U); + graph.append(std::make_shared(), 0, 0.0, /*gate_index=*/11); + BOOST_REQUIRE_EQUAL(graph.layers(), 1U); + BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 11U); +} + BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { std::vector layers; for (std::size_t g = 10; g < 14; ++g) { From f919da261f006aa8a683922d412c561a442717d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 10:17:02 +0000 Subject: [PATCH 06/13] refactor(c++)!: :fire: drop MPGraph's unused partial-slice machinery front_offset_ is always 0 in production. slice_graph()'s only caller passed key == layers(), so both arms emptied the store and maybe_compact_layers() only ever reached its clear() branch. The 4096 threshold, the half-size test and the erase arm were reachable from mp_graph_tests.cpp alone, which is why their tests go with them. With the dead prefix gone, layers() is layers_.size(), the layer offset is the layer index, and both views base at 0. contraction_view() becomes a one-liner in the header. This also discharges the slice-metadata problem: slice_graph() stamped its source's growth onto a slice whose layer order contradicted it under LayerGrowth::Front, and a prose warning was the only guard. pare_graph is now the sole caller of MPGraph(growth, layers), and it writes each layer back at the index it read it from. BREAKING CHANGE: MPGraph::slice_graph() and MPGraph::slice_view() are removed. Use contraction_view() for the replay and clear() for the drain. Assisted-by: ClaudeCode:claude-opus-5 --- cpp/include/monoprop/MPGraph.h | 69 +++++++--------------- cpp/monoprop/MPGraph.cpp | 79 ++----------------------- cpp/tests/mp_graph_tests.cpp | 105 +-------------------------------- 3 files changed, 29 insertions(+), 224 deletions(-) diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index bc65c9dd..21b991e2 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -37,46 +37,25 @@ enum class LayerGrowth : uint8_t { /// Ordered per-rank record of the evolution circuit, one Layer per generator. // A graph built by append() stores its layers in DESCENDING optimizer-slot order under either growth end: -// active layer i is optimizer slot layers()-1-i. append_layer() normalizes the two arrival orders into that -// one storage order, which is what lets everything reconstructing optimizer order from a graph -- the +// layer i is optimizer slot layers()-1-i. append_layer() normalizes the two arrival orders into that one +// storage order, which is what lets everything reconstructing optimizer order from a graph -- the // evolved-operator setup and the gradient loop in MPFunctions, MonomialPropagator::graph_gate_arrays_ -- do // it without knowing which growth end this graph has. Changing the storage order means changing those too. class monoprop_EXPORT MPGraph { private: - using LayerIterator = std::vector::iterator; - using ConstLayerIterator = std::vector::const_iterator; - LayerGrowth growth_; std::vector layers_; - size_t front_offset_ = 0; - - auto active_begin_index() const -> size_t { return front_offset_; } - - auto active_end_index() const -> size_t { return layers_.size(); } - - auto active_begin_iterator() -> LayerIterator { - return layers_.begin() + static_cast(active_begin_index()); - } - - auto active_end_iterator() -> LayerIterator { return layers_.end(); } - - auto active_begin_iterator() const -> ConstLayerIterator { - return layers_.begin() + static_cast(active_begin_index()); - } - - auto active_end_iterator() const -> ConstLayerIterator { return layers_.end(); } - auto append_position() -> LayerIterator { - return grows_at_front() ? active_begin_iterator() : active_end_iterator(); + auto append_position() -> std::vector::iterator { + return grows_at_front() ? layers_.begin() : layers_.end(); } auto append_layer(Layer layer) -> void { layers_.emplace(append_position(), std::move(layer)); } - auto checked_layer_offset(size_t layer_idx) const -> size_t { + auto check_layer_index_(size_t layer_idx) const -> void { if (layer_idx >= layers()) { throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers())); } - return active_begin_index() + layer_idx; } public: @@ -97,38 +76,32 @@ class monoprop_EXPORT MPGraph { } /// Drop every layer. The graph stays usable, and a later append() starts from an empty store. - auto clear() -> void { - layers_.clear(); - front_offset_ = 0; - } - - /// Slice the graph at `key` (the number of earliest operations to include); `contract` also removes - /// the sliced part from this graph. - // Layers come back in application order, earliest first, which matches this class's storage order only - // for LayerGrowth::Back. A slice is for replay_view() and must not be fed to anything that - // reconstructs optimizer order from it. - auto slice_graph(size_t key, bool contract = false) -> MPGraph; + auto clear() -> void { layers_.clear(); } - auto slice_view(size_t key) const -> MPGraphView; + auto layers() const -> size_t { return layers_.size(); } - auto layers() const -> size_t { return active_end_index() - active_begin_index(); } - - auto get_layer(size_t layer_idx) -> Layer& { return layers_[checked_layer_offset(layer_idx)]; } + auto get_layer(size_t layer_idx) -> Layer& { + check_layer_index_(layer_idx); + return layers_[layer_idx]; + } - auto get_layer(size_t layer_idx) const -> const Layer& { return layers_[checked_layer_offset(layer_idx)]; } + auto get_layer(size_t layer_idx) const -> const Layer& { + check_layer_index_(layer_idx); + return layers_[layer_idx]; + } auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } - /// Non-owning replay view over the active layers, in stored (descending optimizer-slot) order. - auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } + /// Non-owning replay view over the layers, in stored (descending optimizer-slot) order. + auto replay_view() const -> MPGraphView { return {layers_, 0, layers(), false}; } - /// The active layers in the order this graph's own build walked them, which is the order a - /// contraction must replay them in. + /// The layers in the order this graph's own build walked them, which is the order a contraction must + /// replay them in. // Not replay_view(): a contraction drives the live coefficient vector, so it follows the simulation // direction, where the picture-free evaluation order is always the stored one. - auto contraction_view() const -> MPGraphView; + auto contraction_view() const -> MPGraphView { return {layers_, 0, layers(), grows_at_front()}; } - /// Carried so a derived graph (a slice, a pared copy) keeps its source's layer order. + /// Carried so a pared copy keeps its source's layer order. auto growth() const -> LayerGrowth { return growth_; } /// Whether new layers attach at the front, so the oldest operation is at the back. diff --git a/cpp/monoprop/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index 18140304..95be3767 100644 --- a/cpp/monoprop/MPGraph.cpp +++ b/cpp/monoprop/MPGraph.cpp @@ -14,14 +14,8 @@ #include "monoprop/MPGraph.h" -#include -#include -#include #include -#include - -#include -#include +#include #include "monoprop/TypeAliases.h" @@ -29,24 +23,6 @@ namespace monoprop { namespace { -// Erase the dead front prefix only once it is both large and >= half the vector, to bound amortized cost. -auto maybe_compact_layers(std::vector &layers, size_t &front_offset) -> void { - if (front_offset == 0) { - return; - } - - if (front_offset >= layers.size()) { - layers.clear(); - front_offset = 0; - return; - } - - if (front_offset >= 4096 && front_offset * 2 >= layers.size()) { - layers.erase(layers.begin(), layers.begin() + static_cast(front_offset)); - front_offset = 0; - } -} - auto layer_storage_memory_usage(const LayerCore &storage) -> GraphMemoryBreakdown { GraphMemoryBreakdown breakdown; breakdown.layer_storage_object_bytes = sizeof(LayerCore); @@ -57,53 +33,10 @@ auto layer_storage_memory_usage(const LayerCore &storage) -> GraphMemoryBreakdow } // namespace -auto MPGraph::slice_graph(size_t key, bool contract) -> MPGraph { - std::vector sliced_layers; - const auto k = std::min(key, layers()); - sliced_layers.reserve(k); - - if (grows_at_front()) { - // Oldest end at the back, so the slice is the tail read backwards; active_end_index() is - // layers_.size(), which is what makes rbegin() the right starting point. - sliced_layers.insert(sliced_layers.end(), layers_.rbegin(), layers_.rbegin() + static_cast(k)); - - if (contract && k != 0) { - layers_.resize(active_end_index() - k); - } - } - else { - const size_t active_begin = active_begin_index(); - const size_t slice_end = active_begin + k; - const auto begin = active_begin_iterator(); - sliced_layers.insert(sliced_layers.end(), begin, begin + static_cast(k)); - - if (contract && k != 0) { - front_offset_ = slice_end; - maybe_compact_layers(layers_, front_offset_); - } - } - - return MPGraph(growth_, std::move(sliced_layers)); -} - -auto MPGraph::contraction_view() const -> MPGraphView { - // One window, two directions: active_end_index() - layers() is active_begin_index(), so only the - // reverse flag distinguishes the growth ends. - return {layers_, active_begin_index(), layers(), grows_at_front()}; -} - -auto MPGraph::slice_view(size_t key) const -> MPGraphView { - const auto k = std::min(key, layers()); - // The window sits at the oldest end and is walked away from it, so the view yields the earliest - // operation first either way -- which makes the reverse flag exactly the growth predicate. - const bool front = grows_at_front(); - return {layers_, front ? active_end_index() - k : active_begin_index(), k, front}; -} - auto MPGraph::total_cycles() const -> size_t { size_t total = 0; - for (auto it = active_begin_iterator(); it != active_end_iterator(); ++it) { - total += it->traversal().total_cycles(); + for (const auto &layer : layers_) { + total += layer.traversal().total_cycles(); } return total; } @@ -113,12 +46,12 @@ auto MPGraph::storage_memory_usage() const -> GraphMemoryBreakdown { breakdown.layer_descriptor_bytes = layers_.capacity() * sizeof(Layer); std::unordered_set seen_storage; - for (auto it = active_begin_iterator(); it != active_end_iterator(); ++it) { - if (const auto storage = it->shared_core(); storage != nullptr && seen_storage.insert(storage.get()).second) { + for (const auto &layer : layers_) { + if (const auto storage = layer.shared_core(); storage != nullptr && seen_storage.insert(storage.get()).second) { breakdown += layer_storage_memory_usage(*storage); } // Pruned cos is owned per-layer, not by the shared core, so it accumulates without the dedup. - if (const CosMask *cos = it->pruned_cos(); cos != nullptr) { + if (const CosMask *cos = layer.pruned_cos(); cos != nullptr) { breakdown.cos_data_bytes += cos->blocks.capacity() * sizeof(std::pair); } } diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index dfcedd51..64d1c727 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -12,8 +12,8 @@ // See the License for the specific language governing permissions and // limitations under the License. -// White-box tests for MPGraph transforms and MPGraphView, built by direct Layer construction -// (GraphBuildHarness). Each layer's distinct gate_index is the oracle for slice / view ordering. +// White-box tests for MPGraph's views and MPGraphView, built by direct Layer construction +// (GraphBuildHarness). Each layer's distinct gate_index is the oracle for view ordering. #include @@ -28,96 +28,6 @@ using test_utils::core_with_gate; using test_utils::graph_with_gates; using test_utils::layer_with_gate; -// LayerGrowth::Back is the Heisenberg order. -BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_growth_back_prefix_no_contract) { - auto graph = graph_with_gates(LayerGrowth::Back, 5); // layers_ = [0,1,2,3,4] - auto sliced = graph.slice_graph(3, /*contract=*/false); - - BOOST_REQUIRE_EQUAL(sliced.layers(), 3U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(0).gate_index(), 0U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(1).gate_index(), 1U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(2).gate_index(), 2U); - // Non-contracting slice leaves the source untouched. - BOOST_CHECK_EQUAL(graph.layers(), 5U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 0U); -} - -// LayerGrowth::Front is the Schrödinger order. -BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_growth_front_contract_newest_first_copy_and_resize) { - // Front growth stores newest-first: appending gates 0..4 gives layers_ = [4,3,2,1,0]. - auto graph = graph_with_gates(LayerGrowth::Front, 5); - auto sliced = graph.slice_graph(2, /*contract=*/true); - - // sliced = layers_[active_end-1-i] = layers_[4], layers_[3] = gates 0, 1 (oldest-first). - BOOST_REQUIRE_EQUAL(sliced.layers(), 2U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(0).gate_index(), 0U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(1).gate_index(), 1U); - - // Contract resized layers_ to the newest 3 (gates 4,3,2, still newest-first). - BOOST_REQUIRE_EQUAL(graph.layers(), 3U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 4U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(1).gate_index(), 3U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(2).gate_index(), 2U); -} - -BOOST_AUTO_TEST_CASE(mp_graph_slice_graph_key_clamped_to_size) { - auto graph = graph_with_gates(LayerGrowth::Back, 3); - auto sliced = graph.slice_graph(100, /*contract=*/false); - BOOST_CHECK_EQUAL(sliced.layers(), 3U); -} - -// The maybe_compact_layers arms below are reached through back-growth slice_graph(contract=true). - -BOOST_AUTO_TEST_CASE(mp_graph_contract_clear_arm_when_prefix_covers_all) { - auto graph = graph_with_gates(LayerGrowth::Back, 5); - (void)graph.slice_graph(5, /*contract=*/true); // front_offset == size -> clear - BOOST_CHECK_EQUAL(graph.layers(), 0U); - // Graph is still usable after a full clear. - graph.append(std::make_shared(), 0, 0.0, /*gate_index=*/42); - BOOST_REQUIRE_EQUAL(graph.layers(), 1U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 42U); -} - -BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { - auto graph = graph_with_gates(LayerGrowth::Back, 100); - (void)graph.slice_graph(3, /*contract=*/true); // front_offset 3 < 4096 -> no physical compaction - BOOST_REQUIRE_EQUAL(graph.layers(), 97U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 3U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(96).gate_index(), 99U); -} - -BOOST_AUTO_TEST_CASE(mp_graph_contract_erase_arm_above_threshold) { - // The erase arm fires only when front_offset >= 4096 AND 2*front_offset >= size. - auto graph = graph_with_gates(LayerGrowth::Back, 8200); - auto sliced = graph.slice_graph(4100, /*contract=*/true); - BOOST_CHECK_EQUAL(sliced.layers(), 4100U); - BOOST_CHECK_EQUAL(sliced.get_layer_traversal(0).gate_index(), 0U); - - BOOST_REQUIRE_EQUAL(graph.layers(), 4100U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 4100U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(4099).gate_index(), 8199U); -} - -BOOST_AUTO_TEST_CASE(mp_graph_slice_view_growth_back_forward_window) { - auto graph = graph_with_gates(LayerGrowth::Back, 5); - auto view = graph.slice_view(3); - BOOST_REQUIRE_EQUAL(view.layers(), 3U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(1).gate_index(), 1U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(2).gate_index(), 2U); -} - -BOOST_AUTO_TEST_CASE(mp_graph_slice_view_growth_front_reversed_window) { - // layers_ = [4,3,2,1,0]; slice_view(3) uses base=active_end-3=2, reverse=true. - // get_layer_traversal(i) -> layers_[2 + (3-1-i)] -> gates 0,1,2 in replay order. - auto graph = graph_with_gates(LayerGrowth::Front, 5); - auto view = graph.slice_view(3); - BOOST_REQUIRE_EQUAL(view.layers(), 3U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(0).gate_index(), 0U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(1).gate_index(), 1U); - BOOST_CHECK_EQUAL(view.get_layer_traversal(2).gate_index(), 2U); -} - // replay_view() is the evaluation order: stored order, so the growth end shows through it. BOOST_AUTO_TEST_CASE(mp_graph_replay_view_is_stored_order) { auto back = graph_with_gates(LayerGrowth::Back, 4); // layers_ = [0,1,2,3] @@ -154,17 +64,6 @@ BOOST_AUTO_TEST_CASE(mp_graph_clear_empties_and_leaves_the_graph_usable) { BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 7U); } -// clear() must also drop the dead front prefix a contracting slice leaves behind. -BOOST_AUTO_TEST_CASE(mp_graph_clear_resets_the_dead_front_prefix) { - auto graph = graph_with_gates(LayerGrowth::Back, 100); - (void)graph.slice_graph(3, /*contract=*/true); // front_offset 3, below the compaction threshold - graph.clear(); - BOOST_REQUIRE_EQUAL(graph.layers(), 0U); - graph.append(std::make_shared(), 0, 0.0, /*gate_index=*/11); - BOOST_REQUIRE_EQUAL(graph.layers(), 1U); - BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 11U); -} - BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { std::vector layers; for (std::size_t g = 10; g < 14; ++g) { From 3be42bbd0a6088e74306f3b3905d543926c98b76 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 10:21:43 +0000 Subject: [PATCH 07/13] refactor(c++)!: :broom: give MPGraphView a span instead of a container pointer The (container pointer, base, count) trio was reimplementing std::span. With the dead front prefix gone the window is always the whole store, so base_ has nothing left to carry and the caller folds the window into the span. sizeof(MPGraphView) drops from 40 to 24 bytes and the type becomes trivially copyable; it is copied by value into an EvalRequest on every evaluation. More importantly the class stops naming std::vector, so a later change to how MPGraph addresses its store does not have to reach into the view. Also records why this is not a std::ranges adaptor: the consumers index it rather than iterate it (it has no begin()/end()), one type has to serve both directions and std::ranges::any_view is C++26, and the type sits in exported signatures where a ranges adaptor would mangle library-version-dependently. BREAKING CHANGE: MPGraphView's constructor takes (std::span, bool) instead of (const std::vector &, size_t, size_t, bool). Assisted-by: ClaudeCode:claude-opus-5 --- cpp/include/monoprop/MPGraph.h | 4 +-- cpp/monoprop/detail/graph/MPGraphViews.h | 35 ++++++++++++++---------- cpp/tests/mp_graph_tests.cpp | 4 +-- 3 files changed, 24 insertions(+), 19 deletions(-) diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 21b991e2..b3ac71cf 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -93,13 +93,13 @@ class monoprop_EXPORT MPGraph { auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } /// Non-owning replay view over the layers, in stored (descending optimizer-slot) order. - auto replay_view() const -> MPGraphView { return {layers_, 0, layers(), false}; } + auto replay_view() const -> MPGraphView { return {layers_, false}; } /// The layers in the order this graph's own build walked them, which is the order a contraction must /// replay them in. // Not replay_view(): a contraction drives the live coefficient vector, so it follows the simulation // direction, where the picture-free evaluation order is always the stored one. - auto contraction_view() const -> MPGraphView { return {layers_, 0, layers(), grows_at_front()}; } + auto contraction_view() const -> MPGraphView { return {layers_, grows_at_front()}; } /// Carried so a pared copy keeps its source's layer order. auto growth() const -> LayerGrowth { return growth_; } diff --git a/cpp/monoprop/detail/graph/MPGraphViews.h b/cpp/monoprop/detail/graph/MPGraphViews.h index 0af7d7a2..3128a082 100644 --- a/cpp/monoprop/detail/graph/MPGraphViews.h +++ b/cpp/monoprop/detail/graph/MPGraphViews.h @@ -15,6 +15,7 @@ #pragma once #include +#include #include #include #include @@ -56,34 +57,38 @@ struct GraphMemoryBreakdown final { } }; -// `reverse` traverses the window newest-first (Schrödinger replay order). Non-owning — the layer vector -// must outlive the view. +// A bounds-checked, optionally reversed index space over a window of layers. `reverse` traverses it +// newest-first. Non-owning: the layer storage must outlive the view. +// +// Not a std::ranges adaptor, and not iterable, for three reasons: +// - Index i addresses three parallel things at once -- the layer, params[i] in evolve_operator(), and +// the recipe cache in build_cos_callbacks(). ev_and_grad() passes a computed index, not an iteration +// step. The consumers need an index space, not a sequence. +// - One type must serve both directions. views::counted(...) and views::reverse(views::counted(...)) +// are different types, and std::ranges::any_view is C++26. +// - The type appears in the exported signatures of ev(), ev_and_grad(), evolve_operator() and +// state_operator_derivative_local(). A ranges adaptor would put a library-version-dependent +// template soup into their mangled names. class MPGraphView { public: - MPGraphView(const std::vector &layers, size_t base, size_t count, bool reverse) - : layers_(&layers), - base_(base), - count_(count), - reverse_(reverse) {} + MPGraphView(std::span layers, bool reverse) : layers_(layers), reverse_(reverse) {} - auto layers() const -> size_t { return count_; } + auto layers() const -> size_t { return layers_.size(); } - auto get_layer(size_t layer_idx) const -> const Layer & { return (*layers_)[checked_layer_offset(layer_idx)]; } + auto get_layer(size_t layer_idx) const -> const Layer & { return layers_[checked_layer_offset(layer_idx)]; } auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } private: auto checked_layer_offset(size_t layer_idx) const -> size_t { - if (layer_idx >= count_) { - throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, count_)); + if (layer_idx >= layers_.size()) { + throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers_.size())); } - return base_ + (reverse_ ? count_ - 1 - layer_idx : layer_idx); + return reverse_ ? layers_.size() - 1 - layer_idx : layer_idx; } - const std::vector *layers_ = nullptr; - size_t base_ = 0; - size_t count_ = 0; + std::span layers_; bool reverse_ = false; }; diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 64d1c727..20991590 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -70,8 +70,8 @@ BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { layers.push_back(layer_with_gate(g)); // [10,11,12,13] } - const MPGraphView fwd(layers, /*base=*/0, /*count=*/4, /*reverse=*/false); - const MPGraphView rev(layers, /*base=*/0, /*count=*/4, /*reverse=*/true); + const MPGraphView fwd(layers, /*reverse=*/false); + const MPGraphView rev(layers, /*reverse=*/true); for (std::size_t i = 0; i < 4; ++i) { BOOST_CHECK_EQUAL(fwd.get_layer_traversal(i).gate_index(), 10U + i); BOOST_CHECK_EQUAL(rev.get_layer_traversal(i).gate_index(), 13U - i); From 34a5f39577af7d21cf785db933c6e46e934b9260 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 10:25:54 +0000 Subject: [PATCH 08/13] refactor(c++): :broom: name MPGraph's store-to-slot map Four sites open-coded `n - 1 - x` to convert between store order and optimizer order, each deriving MPGraph's storage invariant from a prose comment. They now call slot_of_layer(), which is its own inverse, so one name serves both directions and two static_asserts pin it. Two other sites share the arithmetic but not the meaning, and keep their own names: Picture::gate_slot maps a simulation step to a slot (it coincides only under Heisenberg, where the sibling policy returns `i`), and MPGraphView::checked_layer_offset reverses a window. Also corrects fill_mapped_params' comment, which said the graph is traversed in simulation order. That holds for contract_partially's view, not for the evaluation path's, which is always the stored order. Assisted-by: ClaudeCode:claude-opus-5 --- cpp/include/monoprop/MPGraph.h | 10 ++++++++++ cpp/monoprop/MPFunctions.cpp | 8 ++++---- .../detail/monomial_propagator/MonomialPropagator.inl | 5 ++--- 3 files changed, 16 insertions(+), 7 deletions(-) diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index b3ac71cf..37f4ec6f 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -35,6 +35,16 @@ enum class LayerGrowth : uint8_t { Front, ///< a new gate takes the highest optimizer slot, so it attaches at the front }; +/// The optimizer slot of stored layer `layer_idx` in an `n`-layer graph. +// Its own inverse, so the slot-to-layer direction calls it too. The one spelling of MPGraph's storage +// invariant: every conversion between store order and optimizer order goes through it. +constexpr auto slot_of_layer(size_t layer_idx, size_t n) -> size_t { + return n - 1 - layer_idx; +} + +static_assert(slot_of_layer(0, 4) == 3); +static_assert(slot_of_layer(slot_of_layer(1, 4), 4) == 1); + /// Ordered per-rank record of the evolution circuit, one Layer per generator. // A graph built by append() stores its layers in DESCENDING optimizer-slot order under either growth end: // layer i is optimizer slot layers()-1-i. append_layer() normalizes the two arrival orders into that one diff --git a/cpp/monoprop/MPFunctions.cpp b/cpp/monoprop/MPFunctions.cpp index b716aad2..0456fc3d 100644 --- a/cpp/monoprop/MPFunctions.cpp +++ b/cpp/monoprop/MPFunctions.cpp @@ -52,8 +52,8 @@ auto eval_scratch() -> EvalScratch & { return scratch; } -// Graph is traversed in simulation order but parameter_mapping is stored in optimizer order; write the -// mapped coefficients forward or reversed accordingly. +// parameter_mapping is indexed by optimizer slot; `result` is indexed by the position the caller's graph +// view traverses. `reverse` says the two disagree, and disagreement is exactly slot_of_layer. auto fill_mapped_params(VecD &result, const VecD ¶meters, const VecZ ¶meter_mapping, @@ -63,7 +63,7 @@ auto fill_mapped_params(VecD &result, const size_t count = parameter_mapping.size(); result.resize(count); for (size_t i = 0; i < count; ++i) { - const size_t dst = reverse ? (count - 1 - i) : i; + const size_t dst = reverse ? slot_of_layer(i, count) : i; result[dst] = phase * parameters[parameter_mapping[i]] * gen_coeffs[i]; } } @@ -229,7 +229,7 @@ auto ev_and_grad(const EvalRequest &request, mpi::Comm comm, const detail::CosCa const auto ¶meter_mapping = request.parameter_mapping; scratch.gradient.assign(request.params.size(), 0.0); for (size_t i = 0; i < parameter_mapping.size(); ++i) { - const auto idx = parameter_mapping.size() - 1 - i; + const auto idx = slot_of_layer(i, parameter_mapping.size()); const auto param_ind = parameter_mapping[i]; scratch.gradient[param_ind] += state_operator_derivative_local(state_, diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 58cc1b72..8e04832b 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -845,7 +845,7 @@ auto MonomialPropagator::set_parameter_mapping(const VecZ ¶meter_m if (parameter_mapping.size() == count) { // Per-layer mapping in optimizer order. for (size_t layer = 0; layer < count; ++layer) { - relabel(layer, parameter_mapping[count - 1 - layer]); + relabel(layer, parameter_mapping[slot_of_layer(layer, count)]); } } else if (parameter_mapping.size() == gates) { @@ -871,10 +871,9 @@ auto MonomialPropagator::graph_gate_arrays_() const -> std::pair Date: Tue, 18 Aug 2026 10:30:43 +0000 Subject: [PATCH 09/13] refactor(c++)!: :broom: let pare_graph name its own sweep direction pare_graph read graph.grows_at_front() to choose where its reachability sweep starts. That is not a layer-order question. The sweep must start at the end of the replay where the paring seed sits: Heisenberg seeds from the state, which contracts against the replay's last layer, and Schrodinger seeds from the operator, which enters at the first. The growth bit answered it only because each picture's build direction and its seed vector happen to agree. It now takes a PareSweep, supplied by the policy as `pare_sweep` beside `pare_seed`, which is what it describes. The old comment there was also wrong for LayerGrowth::Back: it claimed new layers attach where the circuit's latest operations are, but under back growth the last-appended layer is the circuit's earliest gate. grows_at_front() has no caller outside MPGraph again, so it goes back to private. BREAKING CHANGE: pare_graph() takes a trailing PareSweep argument. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 4 ++-- cpp/include/monoprop/MPFunctions.h | 12 +++++++++++- cpp/include/monoprop/MPGraph.h | 9 ++++----- .../monomial_propagator/MonomialPropagator.inl | 3 ++- cpp/monoprop/detail/pare/PareGraph.cpp | 9 +++++---- cpp/monoprop/picture/Picture.h | 5 +++++ cpp/tests/pare_graph_tests.cpp | 3 ++- 7 files changed, 31 insertions(+), 14 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 8f9caec2..9daf96c3 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -70,8 +70,8 @@ Key files: - **`Picture` / the picture policy** (`cpp/monoprop/core/Picture.h`, `cpp/monoprop/picture/Picture.h`): the two simulation pictures are sibling models (`HeisenbergPicture`, `SchrodingerPicture`), built the same way as the algebras. Each states one picture's whole rule set — gate traversal direction, applied-angle sign, - `map_params` phase, the live coefficient vector, the contraction partner — so no `if (schrodinger)` is - written twice. Every public entry point of `MonomialPropagator` binds the policy once with `with_picture`; + `map_params` phase, the live coefficient vector, the contraction partner, the end its paring sweep starts + from — so no `if (schrodinger)` is written twice. Every public entry point of `MonomialPropagator` binds the policy once with `with_picture`; its whole private layer is templated on that policy and never re-tests which picture it is in, so there is no runtime-dispatching helper layer. The fused `ContractSink`/`apply_fused_contract` pair likewise takes the policy as a template parameter, bound once inside `build_layer` — at the sink only, or the diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index 087fd0a8..521af005 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include #include @@ -106,10 +107,19 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, mpi::Comm comm = MPI_COMM_WORLD, const detail::CosCallbacks &cos = {}) -> std::pair; +/// Which end of the replay holds the paring seed, so which end the reachability sweep starts from. +// Not the graph's layer order: it is a property of the seed vector, and the two only happen to agree +// because each picture's build direction and its seed vector do. +enum class PareSweep : uint8_t { + FromOutput, ///< the seed contracts against the replay's last layer (a Heisenberg state) + FromInput, ///< the seed enters at the replay's first layer (a Schrödinger operator) +}; + /// Prune `graph` to the subgraph reaching `nonzero_inds`; `full_cos_of_layer(i)` supplies layer i's full cosine set. monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, mpi::Comm comm, - const std::function &full_cos_of_layer) -> MPGraph; + const std::function &full_cos_of_layer, + PareSweep sweep) -> MPGraph; } // namespace monoprop diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 37f4ec6f..946e36e8 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -56,6 +56,10 @@ class monoprop_EXPORT MPGraph { LayerGrowth growth_; std::vector layers_; + // The one spelling of the growth question: the ordering-sensitive members ask it this way rather + // than comparing enumerators. + auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } + auto append_position() -> std::vector::iterator { return grows_at_front() ? layers_.begin() : layers_.end(); } @@ -114,11 +118,6 @@ class monoprop_EXPORT MPGraph { /// Carried so a pared copy keeps its source's layer order. auto growth() const -> LayerGrowth { return growth_; } - /// Whether new layers attach at the front, so the oldest operation is at the back. - // The one spelling of the growth question: every ordering-sensitive site, here and in pare_graph, asks - // it this way rather than comparing enumerators. - auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } - /// A normally-built layer stores no cosine set, so the companion cosine-index count cannot come from /// the graph: only the operator's inverted index can supply it. auto total_cycles() const -> size_t; diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 8e04832b..e64ac4e8 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -964,7 +964,8 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optional(pare_graph(graph_, keep, count, comm_, full_cos_of_layer)); + graph = + std::make_shared(pare_graph(graph_, keep, count, comm_, full_cos_of_layer, P::pare_sweep)); } else { graph = std::shared_ptr(std::shared_ptr{}, &graph_); diff --git a/cpp/monoprop/detail/pare/PareGraph.cpp b/cpp/monoprop/detail/pare/PareGraph.cpp index 1403d241..743e4c17 100644 --- a/cpp/monoprop/detail/pare/PareGraph.cpp +++ b/cpp/monoprop/detail/pare/PareGraph.cpp @@ -104,7 +104,8 @@ auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, mpi::Comm comm, - const std::function &full_cos_of_layer) -> MPGraph { + const std::function &full_cos_of_layer, + PareSweep sweep) -> MPGraph { const size_t num_layers = graph.layers(); const auto my_rank = static_cast(mpi::rank(comm)); @@ -121,9 +122,9 @@ auto pare_graph(const MPGraph &graph, // mark_cross_rank_endpoints_kept), so nodes_to_keep stays consistent across ranks with no exchange. // Cross-rank lists are never pruned; the keep-set only has to be right so cos pruning stays exact. for (size_t iter = 0; iter < num_layers; ++iter) { - // Backward means last-applied first, which is the growth end: new layers attach where the - // circuit's latest operations are. - const size_t layer_idx = graph.grows_at_front() ? iter : (num_layers - 1 - iter); + // The sweep walks away from the seed end, so reachability reaches every kept node before the + // cosine filter runs on the layer that produced it. + const size_t layer_idx = (sweep == PareSweep::FromInput) ? iter : slot_of_layer(iter, num_layers); const auto &layer = graph.get_layer(layer_idx); const auto lt = layer.traversal(); diff --git a/cpp/monoprop/picture/Picture.h b/cpp/monoprop/picture/Picture.h index 208fda70..aa841951 100644 --- a/cpp/monoprop/picture/Picture.h +++ b/cpp/monoprop/picture/Picture.h @@ -82,6 +82,8 @@ struct HeisenbergPicture { static auto pare_seed(const EvalState &state, const VecD & /*op*/, double threshold) -> std::pair { return {state.indices_above(threshold), state.length()}; } + // The state contracts against the evolved operator, so the seed sits at the replay's output end. + static constexpr PareSweep pare_sweep = PareSweep::FromOutput; // A perf hint, never a correctness constraint: overflow spills losslessly. The bound is already in // physical slots (CutoffEvaluator::max_slot_bound), so nothing to scale. NumModes is explicit @@ -135,6 +137,8 @@ struct SchrodingerPicture { static auto pare_seed(const EvalState & /*state*/, const VecD &op, double threshold) -> std::pair { return {indices_above(op, threshold), op.size()}; } + // `op` is what the replay evolves, so the seed sits at the replay's input end. + static constexpr PareSweep pare_sweep = PareSweep::FromInput; // The state's monomials come from generate_paired_op(), not from cutoff_fn_, so the cutoff carries // no structural bound on them. @@ -153,6 +157,7 @@ concept PicturePolicy = requires { { P::contract_phase } -> std::convertible_to; { P::contract_reverse } -> std::convertible_to; { P::layer_growth } -> std::convertible_to; + { P::pare_sweep } -> std::convertible_to; }; static_assert(PicturePolicy); diff --git a/cpp/tests/pare_graph_tests.cpp b/cpp/tests/pare_graph_tests.cpp index 2464b2d0..b826cf07 100644 --- a/cpp/tests/pare_graph_tests.cpp +++ b/cpp/tests/pare_graph_tests.cpp @@ -105,7 +105,8 @@ BOOST_AUTO_TEST_CASE(pare_graph_emits_expected_layer_kinds) { seed.push_back(i); } - auto pared = pare_graph(graph, seed, local_index_count, MPI_COMM_SELF, provider); + // Heisenberg's sweep: the simulator built above is a Heisenberg one. + auto pared = pare_graph(graph, seed, local_index_count, MPI_COMM_SELF, provider, PareSweep::FromOutput); BOOST_REQUIRE_EQUAL(pared.layers(), graph.layers()); size_t pruned_count = 0; From 4759fa8dd6181cc0d57e32157ad878b52e296f6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 10:40:19 +0000 Subject: [PATCH 10/13] refactor(c++)!: :zap: store MPGraph layers in arrival order, reverse on read layers_ grew at either end so that layer indices could stay in descending optimizer-slot order under both build directions. An ascending-slot build therefore inserted at layers_.begin() once per gate, which is O(n) per append. The store now only ever push_backs, and get_layer() maps a layer index onto arrival order instead. The layer order is unchanged, so nothing downstream moves; measured on this machine, appending 80k layers goes from 2863 ms to 15 ms (20k/40k/80k: 177/716/2863 ms before, 4.6/7.5/14.7 ms after -- quadratic to linear). Two consequences make the class simpler rather than cleverer: - contraction_view() is never reversed now, under either build direction: arrival order IS build order. replay_view() carries the flag instead, since the evaluation order is the layer order. - MPGraph(order, std::vector) is gone. pare_graph copies the graph and calls the new replace_layer() instead of default-building a vector of n Layers, each of which allocated a LayerCore the sweep immediately overwrote. One constructor also means no object can carry metadata that contradicts its own layout. LayerGrowth becomes ArrivalOrder{DescendingSlot, AscendingSlot} and the policy member becomes arrival_order, which say what gate_slot's slope means rather than describing a store layout that no longer exists. Nothing reads the bit back, so the accessor is gone too: it enters at construction and stays inside the class. BREAKING CHANGE: LayerGrowth is renamed to ArrivalOrder, with its enumerators renamed to DescendingSlot/AscendingSlot; MPGraph::growth() and the MPGraph(ArrivalOrder, std::vector) constructor are removed. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 7 +- cpp/include/monoprop/MPGraph.h | 69 +++++++++---------- .../MonomialPropagator.inl | 2 +- cpp/monoprop/detail/pare/PareGraph.cpp | 9 ++- cpp/monoprop/picture/Picture.h | 10 +-- cpp/tests/GraphBuildHarness.h | 6 +- cpp/tests/mp_graph_tests.cpp | 69 +++++++++++++------ 7 files changed, 100 insertions(+), 72 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 9daf96c3..4794e4cf 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -78,9 +78,10 @@ Key files: `with_algebra` scan above it would instantiate four times per mode width instead of two. The picture is fixed at construction: the constructor takes a `PictureSpec = std::variant`, so only a Schrodinger run carries a state cutoff. `MPGraph` deliberately knows no picture — it takes a - graph-local `LayerGrowth` saying which end a new layer attaches to, and each policy carries its own - `layer_growth` beside `gate_slot`. The graph's layer order itself is picture-independent and load-bearing - beyond the class; `MPGraph`'s own comment records what depends on it. + graph-local `ArrivalOrder` naming the slope of the optimizer slots a build hands to `append()`, and each + policy carries its own `arrival_order` beside `gate_slot`. Nothing reads it back: the bit enters at + construction and stays inside the class. Layer indices are picture-independent and load-bearing beyond + `MPGraph`; its own comment records what depends on them, and `slot_of_layer` is their one spelling. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 946e36e8..8ea224e7 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -28,11 +28,11 @@ namespace monoprop { -/// Which end of the layer store a newly appended gate attaches to. +/// The slope of the optimizer slots a build hands to append(), which is Picture::gate_slot's slope. // This one bit is all the graph needs to know about the simulation that drives it. -enum class LayerGrowth : uint8_t { - Back, ///< a new gate takes the lowest optimizer slot, so it attaches at the back - Front, ///< a new gate takes the highest optimizer slot, so it attaches at the front +enum class ArrivalOrder : uint8_t { + DescendingSlot, ///< each arriving gate takes a lower slot than the last + AscendingSlot, ///< each arriving gate takes a higher slot than the last }; /// The optimizer slot of stored layer `layer_idx` in an `n`-layer graph. @@ -46,26 +46,28 @@ static_assert(slot_of_layer(0, 4) == 3); static_assert(slot_of_layer(slot_of_layer(1, 4), 4) == 1); /// Ordered per-rank record of the evolution circuit, one Layer per generator. -// A graph built by append() stores its layers in DESCENDING optimizer-slot order under either growth end: -// layer i is optimizer slot layers()-1-i. append_layer() normalizes the two arrival orders into that one -// storage order, which is what lets everything reconstructing optimizer order from a graph -- the -// evolved-operator setup and the gradient loop in MPFunctions, MonomialPropagator::graph_gate_arrays_ -- do -// it without knowing which growth end this graph has. Changing the storage order means changing those too. +// LAYER INDICES ARE IN DESCENDING OPTIMIZER-SLOT ORDER under either arrival order: layer i is optimizer +// slot slot_of_layer(i, layers()). That one order is what lets everything reconstructing optimizer order +// from a graph -- the evolved-operator setup and the gradient loop in MPFunctions, +// MonomialPropagator::graph_gate_arrays_ -- do it without knowing how this graph was built. Changing it +// means changing those too. +// +// layers_ holds arrival order, so a push_back is all an append costs; reverse_indexing_() maps a layer +// index onto it. An AscendingSlot build therefore stores its layers backwards, and get_layer() is the only +// place that knows. class monoprop_EXPORT MPGraph { private: - LayerGrowth growth_; + ArrivalOrder arrival_; std::vector layers_; - // The one spelling of the growth question: the ordering-sensitive members ask it this way rather - // than comparing enumerators. - auto grows_at_front() const -> bool { return growth_ == LayerGrowth::Front; } + // True when arrival order runs against layer order, so layer 0 is the last element. + auto reverse_indexing_() const -> bool { return arrival_ == ArrivalOrder::AscendingSlot; } - auto append_position() -> std::vector::iterator { - return grows_at_front() ? layers_.begin() : layers_.end(); + auto stored_offset_(size_t layer_idx) const -> size_t { + check_layer_index_(layer_idx); + return reverse_indexing_() ? slot_of_layer(layer_idx, layers_.size()) : layer_idx; } - auto append_layer(Layer layer) -> void { layers_.emplace(append_position(), std::move(layer)); } - auto check_layer_index_(size_t layer_idx) const -> void { if (layer_idx >= layers()) { throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers())); @@ -73,9 +75,7 @@ class monoprop_EXPORT MPGraph { } public: - explicit MPGraph(LayerGrowth growth) : growth_(growth) {} - - explicit MPGraph(LayerGrowth growth, std::vector layers) : growth_(growth), layers_(std::move(layers)) {} + explicit MPGraph(ArrivalOrder arrival) : arrival_(arrival) {} /// Gate info (param_index, gen_coeff, gate_index) is written onto `storage` here while it is still /// mutable, before it is frozen into the Layer's shared const core. @@ -86,37 +86,32 @@ class monoprop_EXPORT MPGraph { storage->param_index = param_index; storage->gen_coeff = gen_coeff; storage->gate_index = gate_index; - append_layer(Layer(std::move(storage))); + layers_.push_back(Layer(std::move(storage))); } + /// Swap in a rebuilt layer, addressed the same way get_layer() addresses it. + auto replace_layer(size_t layer_idx, Layer layer) -> void { layers_[stored_offset_(layer_idx)] = std::move(layer); } + /// Drop every layer. The graph stays usable, and a later append() starts from an empty store. auto clear() -> void { layers_.clear(); } auto layers() const -> size_t { return layers_.size(); } - auto get_layer(size_t layer_idx) -> Layer& { - check_layer_index_(layer_idx); - return layers_[layer_idx]; - } + auto get_layer(size_t layer_idx) -> Layer& { return layers_[stored_offset_(layer_idx)]; } - auto get_layer(size_t layer_idx) const -> const Layer& { - check_layer_index_(layer_idx); - return layers_[layer_idx]; - } + auto get_layer(size_t layer_idx) const -> const Layer& { return layers_[stored_offset_(layer_idx)]; } auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } - /// Non-owning replay view over the layers, in stored (descending optimizer-slot) order. - auto replay_view() const -> MPGraphView { return {layers_, false}; } + /// Non-owning replay view over the layers, in layer (descending optimizer-slot) order. + auto replay_view() const -> MPGraphView { return {layers_, reverse_indexing_()}; } /// The layers in the order this graph's own build walked them, which is the order a contraction must /// replay them in. - // Not replay_view(): a contraction drives the live coefficient vector, so it follows the simulation - // direction, where the picture-free evaluation order is always the stored one. - auto contraction_view() const -> MPGraphView { return {layers_, grows_at_front()}; } - - /// Carried so a pared copy keeps its source's layer order. - auto growth() const -> LayerGrowth { return growth_; } + // Never reversed, whichever way the build ran: arrival order IS build order, and a contraction drives + // the live coefficient vector, so it follows the simulation. The picture-free evaluation order is the + // layer order instead, which is why replay_view() is the one that carries the flag. + auto contraction_view() const -> MPGraphView { return {layers_, false}; } /// A normally-built layer stores no cosine set, so the companion cosine-index count cannot come from /// the graph: only the operator's inverted index can supply it. diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index e64ac4e8..3be3a0f3 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -86,7 +86,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope : picture_{kind_of(picture)}, comm_{comm}, mp_op_{}, - graph_(with_picture(picture_, []() { return P::layer_growth; })), + graph_(with_picture(picture_, []() { return P::arrival_order; })), cutoff_{cutoff}, lower_atol_{lower_atol}, upper_atol_{upper_atol}, diff --git a/cpp/monoprop/detail/pare/PareGraph.cpp b/cpp/monoprop/detail/pare/PareGraph.cpp index 743e4c17..114aca01 100644 --- a/cpp/monoprop/detail/pare/PareGraph.cpp +++ b/cpp/monoprop/detail/pare/PareGraph.cpp @@ -116,7 +116,9 @@ auto pare_graph(const MPGraph &graph, } } - std::vector layers(num_layers); + // Copied, not default-built: a fresh std::vector(num_layers) allocates a LayerCore per layer + // that the sweep then overwrites, and the copy keeps the source's arrival order for free. + MPGraph pared = graph; // Single backward sweep, entirely rank-local: every cross-rank endpoint is force-kept (see // mark_cross_rank_endpoints_kept), so nodes_to_keep stays consistent across ranks with no exchange. @@ -134,10 +136,11 @@ auto pare_graph(const MPGraph &graph, const CosMask full = full_cos_of_layer(layer_idx); auto [filtered, preserves] = filter_layer_cosine_data(full, nodes_to_keep); - layers[layer_idx] = preserves ? Layer(layer.shared_core()) : Layer(layer.shared_core(), std::move(filtered)); + pared.replace_layer(layer_idx, + preserves ? Layer(layer.shared_core()) : Layer(layer.shared_core(), std::move(filtered))); } - return MPGraph(graph.growth(), std::move(layers)); + return pared; } } // namespace monoprop diff --git a/cpp/monoprop/picture/Picture.h b/cpp/monoprop/picture/Picture.h index aa841951..030f7a39 100644 --- a/cpp/monoprop/picture/Picture.h +++ b/cpp/monoprop/picture/Picture.h @@ -46,8 +46,8 @@ struct HeisenbergPicture { // Simulation step i consumes optimizer slot n-1-i: the observable walks the circuit backwards. static auto gate_slot(size_t i, size_t n) -> size_t { return n - 1 - i; } - // gate_slot's slope, as MPGraph needs it: descending, so each arriving gate takes a lower slot. - static constexpr LayerGrowth layer_growth = LayerGrowth::Back; + // gate_slot's slope, as MPGraph needs it. + static constexpr ArrivalOrder arrival_order = ArrivalOrder::DescendingSlot; static constexpr double apply_sign = 1.0; // the applied angle is the build angle // map_params() arguments for contract_partially: forward phase, written in reverse. @@ -106,8 +106,8 @@ struct SchrodingerPicture { // Simulation step i consumes optimizer slot i: the state walks the circuit front-to-back. static auto gate_slot(size_t i, size_t /*n*/) -> size_t { return i; } - // gate_slot's slope, as MPGraph needs it: ascending, so each arriving gate takes a higher slot. - static constexpr LayerGrowth layer_growth = LayerGrowth::Front; + // gate_slot's slope, as MPGraph needs it. + static constexpr ArrivalOrder arrival_order = ArrivalOrder::AscendingSlot; static constexpr double apply_sign = -1.0; // the applied angle is the negated build angle static constexpr double contract_phase = -1.0; @@ -156,7 +156,7 @@ concept PicturePolicy = requires { { P::apply_sign } -> std::convertible_to; { P::contract_phase } -> std::convertible_to; { P::contract_reverse } -> std::convertible_to; - { P::layer_growth } -> std::convertible_to; + { P::arrival_order } -> std::convertible_to; { P::pare_sweep } -> std::convertible_to; }; diff --git a/cpp/tests/GraphBuildHarness.h b/cpp/tests/GraphBuildHarness.h index a45a7998..5602d5b0 100644 --- a/cpp/tests/GraphBuildHarness.h +++ b/cpp/tests/GraphBuildHarness.h @@ -34,9 +34,9 @@ inline auto layer_with_gate(std::size_t gate_index) -> monoprop::Layer { return monoprop::Layer(core_with_gate(gate_index)); } -// Built via append() so the layer ordering is production's for the growth end under test. -inline auto graph_with_gates(monoprop::LayerGrowth growth, std::size_t n) -> monoprop::MPGraph { - monoprop::MPGraph graph(growth); +// Built via append() so the layer ordering is production's for the arrival order under test. +inline auto graph_with_gates(monoprop::ArrivalOrder arrival, std::size_t n) -> monoprop::MPGraph { + monoprop::MPGraph graph(arrival); for (std::size_t i = 0; i < n; ++i) { graph.append(std::make_shared(), /*param_index=*/0, /*gen_coeff=*/0.0, /*gate_index=*/i); } diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 20991590..4c5c7093 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -28,35 +28,64 @@ using test_utils::core_with_gate; using test_utils::graph_with_gates; using test_utils::layer_with_gate; -// replay_view() is the evaluation order: stored order, so the growth end shows through it. -BOOST_AUTO_TEST_CASE(mp_graph_replay_view_is_stored_order) { - auto back = graph_with_gates(LayerGrowth::Back, 4); // layers_ = [0,1,2,3] - auto front = graph_with_gates(LayerGrowth::Front, 4); // layers_ = [3,2,1,0] - const auto back_view = back.replay_view(); - const auto front_view = front.replay_view(); +// graph_with_gates tags gate_index by arrival, so an AscendingSlot graph reads back in reverse. Every +// case below reads that way, which is the whole storage invariant in one line. +BOOST_AUTO_TEST_CASE(mp_graph_layer_index_reverses_an_ascending_slot_arrival) { + auto descending = graph_with_gates(ArrivalOrder::DescendingSlot, 5); + auto ascending = graph_with_gates(ArrivalOrder::AscendingSlot, 5); + for (std::size_t i = 0; i < 5; ++i) { + BOOST_CHECK_EQUAL(descending.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(ascending.get_layer_traversal(i).gate_index(), 4U - i); + } + // Appending again keeps the mapping: the new layer is the last arrival either way. + descending.append(std::make_shared(), 0, 0.0, /*gate_index=*/5); + ascending.append(std::make_shared(), 0, 0.0, /*gate_index=*/5); + BOOST_CHECK_EQUAL(descending.get_layer_traversal(5).gate_index(), 5U); + BOOST_CHECK_EQUAL(ascending.get_layer_traversal(0).gate_index(), 5U); +} + +// replay_view() is the evaluation order: layer order, so the arrival order shows through it. +BOOST_AUTO_TEST_CASE(mp_graph_replay_view_is_layer_order) { + auto descending = graph_with_gates(ArrivalOrder::DescendingSlot, 4); + auto ascending = graph_with_gates(ArrivalOrder::AscendingSlot, 4); + const auto descending_view = descending.replay_view(); + const auto ascending_view = ascending.replay_view(); for (std::size_t i = 0; i < 4; ++i) { - BOOST_CHECK_EQUAL(back_view.get_layer_traversal(i).gate_index(), i); - BOOST_CHECK_EQUAL(front_view.get_layer_traversal(i).gate_index(), 3U - i); + BOOST_CHECK_EQUAL(descending_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(ascending_view.get_layer_traversal(i).gate_index(), 3U - i); } } -// contraction_view() is the build order, so it yields the same sequence under either growth end. That +// contraction_view() is the build order, so it yields the same sequence under either arrival order. That // equality is the point: a contraction replays the circuit the way the build walked it. -BOOST_AUTO_TEST_CASE(mp_graph_contraction_view_is_build_order_under_either_growth) { - auto back = graph_with_gates(LayerGrowth::Back, 4); - auto front = graph_with_gates(LayerGrowth::Front, 4); - const auto back_view = back.contraction_view(); - const auto front_view = front.contraction_view(); - BOOST_REQUIRE_EQUAL(back_view.layers(), 4U); - BOOST_REQUIRE_EQUAL(front_view.layers(), 4U); +BOOST_AUTO_TEST_CASE(mp_graph_contraction_view_is_build_order_under_either_arrival) { + auto descending = graph_with_gates(ArrivalOrder::DescendingSlot, 4); + auto ascending = graph_with_gates(ArrivalOrder::AscendingSlot, 4); + const auto descending_view = descending.contraction_view(); + const auto ascending_view = ascending.contraction_view(); + BOOST_REQUIRE_EQUAL(descending_view.layers(), 4U); + BOOST_REQUIRE_EQUAL(ascending_view.layers(), 4U); for (std::size_t i = 0; i < 4; ++i) { - BOOST_CHECK_EQUAL(back_view.get_layer_traversal(i).gate_index(), i); - BOOST_CHECK_EQUAL(front_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(descending_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(ascending_view.get_layer_traversal(i).gate_index(), i); + } +} + +// replace_layer() addresses layers the way get_layer() does, which is what keeps pare_graph's sweep from +// writing its filtered layer onto the mirror-image one. +BOOST_AUTO_TEST_CASE(mp_graph_replace_layer_addresses_the_same_layer_as_get_layer) { + for (const auto arrival : {ArrivalOrder::DescendingSlot, ArrivalOrder::AscendingSlot}) { + auto graph = graph_with_gates(arrival, 4); + graph.replace_layer(1, layer_with_gate(99)); + BOOST_CHECK_EQUAL(graph.get_layer_traversal(1).gate_index(), 99U); + // The others are untouched, so nothing was written through the mirror index. + BOOST_CHECK_NE(graph.get_layer_traversal(2).gate_index(), 99U); + BOOST_CHECK_EQUAL(graph.layers(), 4U); } } BOOST_AUTO_TEST_CASE(mp_graph_clear_empties_and_leaves_the_graph_usable) { - auto graph = graph_with_gates(LayerGrowth::Front, 5); + auto graph = graph_with_gates(ArrivalOrder::AscendingSlot, 5); graph.clear(); BOOST_REQUIRE_EQUAL(graph.layers(), 0U); graph.append(std::make_shared(), 0, 0.0, /*gate_index=*/7); @@ -81,7 +110,7 @@ BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { } BOOST_AUTO_TEST_CASE(mp_graph_get_layer_out_of_range_throws) { - auto graph = graph_with_gates(LayerGrowth::Back, 3); + auto graph = graph_with_gates(ArrivalOrder::DescendingSlot, 3); BOOST_CHECK_NO_THROW((void)graph.get_layer(2)); BOOST_CHECK_THROW((void)graph.get_layer(3), std::out_of_range); // const overload takes the same guard. From d5b5bff25dad2fd579ae5dce3ebcad391eeec502 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 12:10:01 +0000 Subject: [PATCH 11/13] refactor(c++)!: :broom: give each order bit and index map one home Cleanup over the picture-policy and MPGraph work, from a review for reuse, simplification, efficiency and altitude. No functional change: all 95 bit-identity fingerprints are unchanged. Index mapping - slot_of_layer, and a new bounds-checked checked_layer_offset, move down into MPGraphViews.h so a graph and the views over it share one spelling of the store<->layer map and one out-of-range diagnostic. Both had a copy, with the same format string, and only MPGraph's went through slot_of_layer -- which its own comment claimed was the single spelling. - HeisenbergPicture::gate_slot calls slot_of_layer instead of open-coding n - 1 - i. Derivable policy bits - is_schrodinger and contract_reverse are derived, not hand-set. The PicturePolicy concept shape-checks exactly these constants and none of the behavioural members, so a picture whose bits disagreed still compiled. - contract_phase is gone: it was apply_sign under a second name. Both are the sign an applied angle carries, and contract_partially replays the very angles the build applied, so they cannot differ. - pare_sweep and PareSweep are gone; MPGraph::layer_of_unbuild_step answers it instead. The sweep starts at the last gate the build applied under either picture -- Heisenberg reaches it via layer n-1 (slot 0, its last step), Schrodinger via layer 0 (slot n-1, also its last step) -- because the seed is the result of the whole evolution, and reachability runs backwards from a result through the circuit that produced it. The removed comment called that agreement a coincidence; it is a derivation, and it was the parameter's only reason to exist. Wasted work - apply_initial_operator_, build_evolve_result_ and propagate_one_ used their policy parameter for one thing: reading P::picture back out as the runtime value their callee already takes. Each duplicated a whole body per mode width for no specialization. - build_graph's seed vector is moved into evolve_mode_graph_with_coeffs_ rather than copied; it is dead at the call site. - FusedApply.h no longer includes picture/Picture.h. apply_fused_contract only tests a dependent P::is_schrodinger, which needs no definition of P, so the include put the evaluation API and the graph behind a hot-path header. Also drops the dead picture() accessor, includes orphaned by the above, and the stale slice_graph/front_offset references the earlier removal left behind in the test README and the graph harness. Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 14 +++-- cpp/include/monoprop/MPFunctions.h | 14 +---- cpp/include/monoprop/MPGraph.h | 36 +++++------ cpp/include/monoprop/MonomialPropagator.h | 17 ++---- .../detail/evolution/layer_build/FusedApply.h | 1 - cpp/monoprop/detail/graph/MPGraphViews.h | 35 +++++++---- .../MonomialPropagator.inl | 59 ++++++++----------- cpp/monoprop/detail/pare/PareGraph.cpp | 9 ++- cpp/monoprop/picture/Picture.h | 22 +++---- cpp/tests/GraphBuildHarness.h | 4 +- cpp/tests/README.md | 5 +- cpp/tests/mp_graph_tests.cpp | 16 ++++- cpp/tests/pare_graph_tests.cpp | 2 +- 13 files changed, 114 insertions(+), 120 deletions(-) diff --git a/AGENTS.md b/AGENTS.md index 4794e4cf..114efd52 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -69,9 +69,9 @@ Key files: templated on the algebra policy and bound to a runtime `Basis` once, via `with_algebra`. - **`Picture` / the picture policy** (`cpp/monoprop/core/Picture.h`, `cpp/monoprop/picture/Picture.h`): the two simulation pictures are sibling models (`HeisenbergPicture`, `SchrodingerPicture`), built the same way as the - algebras. Each states one picture's whole rule set — gate traversal direction, applied-angle sign, - `map_params` phase, the live coefficient vector, the contraction partner, the end its paring sweep starts - from — so no `if (schrodinger)` is written twice. Every public entry point of `MonomialPropagator` binds the policy once with `with_picture`; + algebras. Each states one picture's whole rule set — gate traversal direction, the sign an applied angle + carries (`apply_sign`, which is also the phase a contraction's `map_params` replays it with), the live + coefficient vector, the contraction partner — so no `if (schrodinger)` is written twice. Every public entry point of `MonomialPropagator` binds the policy once with `with_picture`; its whole private layer is templated on that policy and never re-tests which picture it is in, so there is no runtime-dispatching helper layer. The fused `ContractSink`/`apply_fused_contract` pair likewise takes the policy as a template parameter, bound once inside `build_layer` — at the sink only, or the @@ -79,9 +79,11 @@ Key files: fixed at construction: the constructor takes a `PictureSpec = std::variant`, so only a Schrodinger run carries a state cutoff. `MPGraph` deliberately knows no picture — it takes a graph-local `ArrivalOrder` naming the slope of the optimizer slots a build hands to `append()`, and each - policy carries its own `arrival_order` beside `gate_slot`. Nothing reads it back: the bit enters at - construction and stays inside the class. Layer indices are picture-independent and load-bearing beyond - `MPGraph`; its own comment records what depends on them, and `slot_of_layer` is their one spelling. + policy carries its own `arrival_order` beside `gate_slot`. The bit enters at construction and never leaves: + a caller that needs a traversal order asks the graph for it (`replay_view`, `contraction_view`, + `layer_of_unbuild_step`) rather than re-deriving one from the picture — which is why `pare_graph` takes no + sweep direction. Layer indices are picture-independent and load-bearing beyond `MPGraph`; its own comment + records what depends on them, and `slot_of_layer` is their one spelling. - **The partition facade**: `partitions > 1` makes a `MonomialPropagator` a facade over S single-partition propagators, one hash partition each. Every method that fans out must use the private partition vocabulary declared in `MonomialPropagator.h` (`for_each_partition_`, `map_partitions_`, `concat_partitions_` diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index 521af005..8940683a 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -14,7 +14,6 @@ #pragma once -#include #include #include #include @@ -107,19 +106,12 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, mpi::Comm comm = MPI_COMM_WORLD, const detail::CosCallbacks &cos = {}) -> std::pair; -/// Which end of the replay holds the paring seed, so which end the reachability sweep starts from. -// Not the graph's layer order: it is a property of the seed vector, and the two only happen to agree -// because each picture's build direction and its seed vector do. -enum class PareSweep : uint8_t { - FromOutput, ///< the seed contracts against the replay's last layer (a Heisenberg state) - FromInput, ///< the seed enters at the replay's first layer (a Schrödinger operator) -}; - /// Prune `graph` to the subgraph reaching `nonzero_inds`; `full_cos_of_layer(i)` supplies layer i's full cosine set. +// Takes no sweep direction: the seed is the result of the whole evolution under either picture, so the +// sweep is always the graph's own build reversed. MPGraph::layer_of_unbuild_step is that order. monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, mpi::Comm comm, - const std::function &full_cos_of_layer, - PareSweep sweep) -> MPGraph; + const std::function &full_cos_of_layer) -> MPGraph; } // namespace monoprop diff --git a/cpp/include/monoprop/MPGraph.h b/cpp/include/monoprop/MPGraph.h index 8ea224e7..4d74e6af 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -16,8 +16,6 @@ #include #include -#include -#include #include #include #include @@ -35,16 +33,6 @@ enum class ArrivalOrder : uint8_t { AscendingSlot, ///< each arriving gate takes a higher slot than the last }; -/// The optimizer slot of stored layer `layer_idx` in an `n`-layer graph. -// Its own inverse, so the slot-to-layer direction calls it too. The one spelling of MPGraph's storage -// invariant: every conversion between store order and optimizer order goes through it. -constexpr auto slot_of_layer(size_t layer_idx, size_t n) -> size_t { - return n - 1 - layer_idx; -} - -static_assert(slot_of_layer(0, 4) == 3); -static_assert(slot_of_layer(slot_of_layer(1, 4), 4) == 1); - /// Ordered per-rank record of the evolution circuit, one Layer per generator. // LAYER INDICES ARE IN DESCENDING OPTIMIZER-SLOT ORDER under either arrival order: layer i is optimizer // slot slot_of_layer(i, layers()). That one order is what lets everything reconstructing optimizer order @@ -64,14 +52,7 @@ class monoprop_EXPORT MPGraph { auto reverse_indexing_() const -> bool { return arrival_ == ArrivalOrder::AscendingSlot; } auto stored_offset_(size_t layer_idx) const -> size_t { - check_layer_index_(layer_idx); - return reverse_indexing_() ? slot_of_layer(layer_idx, layers_.size()) : layer_idx; - } - - auto check_layer_index_(size_t layer_idx) const -> void { - if (layer_idx >= layers()) { - throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers())); - } + return checked_layer_offset(layer_idx, layers_.size(), reverse_indexing_()); } public: @@ -86,7 +67,7 @@ class monoprop_EXPORT MPGraph { storage->param_index = param_index; storage->gen_coeff = gen_coeff; storage->gate_index = gate_index; - layers_.push_back(Layer(std::move(storage))); + layers_.emplace_back(std::move(storage)); } /// Swap in a rebuilt layer, addressed the same way get_layer() addresses it. @@ -103,6 +84,19 @@ class monoprop_EXPORT MPGraph { auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } + /// The layer that step `step` of an unbuild traversal addresses -- the reverse of the order this graph's + /// own build walked, so step 0 is the last gate the build applied. + // Picture-free, and that is a derivation rather than a coincidence: whichever way the build ran, the + // last arrival is the last gate applied. A reachability sweep seeded on the result of the whole + // evolution (pare_graph) therefore always walks this way, because reachability propagates backwards + // from a result through the circuit that produced it. + // Arrival step -> layer is the inverse of layer -> store offset, and slot_of_layer is its own inverse, + // so this is stored_offset_'s map with the flag negated. Walking it from 0 walks the store backwards + // under either arrival order. + auto layer_of_unbuild_step(size_t step) const -> size_t { + return checked_layer_offset(step, layers_.size(), !reverse_indexing_()); + } + /// Non-owning replay view over the layers, in layer (descending optimizer-slot) order. auto replay_view() const -> MPGraphView { return {layers_, reverse_indexing_()}; } diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index a6151fbf..a32e0c3b 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -211,8 +211,6 @@ class MonomialPropagator { }); } - auto picture() const -> Picture { return picture_; } - auto schrodinger() const -> bool { return picture_ == Picture::Schrodinger; } auto basis() const -> Basis { return basis_; } @@ -277,9 +275,7 @@ class MonomialPropagator { auto evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>>; - auto update_initial_operator(const OperatorDict &op_dict) -> void { - with_picture(picture_, [&]() { this->template apply_initial_operator_

(op_dict); }); - } + auto update_initial_operator(const OperatorDict &op_dict) -> void { apply_initial_operator_(op_dict); } protected: static inline const auto ev_fn = [](const EvalRequest &request, @@ -293,7 +289,6 @@ class MonomialPropagator { /// Distribute op_dict across ranks and apply this rank's share; returns its new (terms, coeffs) /// so caches can refresh. - template auto apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD>; Picture picture_; // immutable after construction: no path switches the picture mid-simulation @@ -403,9 +398,9 @@ class MonomialPropagator { auto validate_cutoff_config_(CutoffType cutoff_type, const std::optional> &basis_change) const -> void; - // Everything below takes the picture as a policy type P, never as a runtime value: each public entry - // point binds it once with with_picture(), and the whole private layer is then written against one - // picture at a time. Nothing here re-tests which picture it is in. + // Each public entry point binds the picture once with with_picture(), and everything below that has a + // picture rule of its own takes it as the policy type P. Nothing here re-tests which picture it is in; + // the few members that carry no rule read picture_ only to forward it to a runtime-valued callee. template auto initialize_operator_caches_() -> void; @@ -435,7 +430,7 @@ class MonomialPropagator { const VecD &gen_coeffs, const VecZ &gate_indices, const VecD ¶meters, - const VecD &operator_coeffs, + VecD operator_coeffs, // by value: the caller's seed is dead after the call std::optional only_rotate_len_k) -> void; // build_layer resolves the same policy for its fused sink, so the cosine sweep and the apply agree. @@ -454,7 +449,6 @@ class MonomialPropagator { std::optional only_rotate_len_k, EvolutionFunc evolution_func) -> void; - template auto propagate_one_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs = std::nullopt, @@ -465,7 +459,6 @@ class MonomialPropagator { // fused_scale_coeffs (ContractImmediately only): the picture's mutable coeff vector for the uncapped // fused cos sweep; the taken decision is reported via fused_scale so the apply matches. See build_layer. - template auto build_evolve_result_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs = std::nullopt, diff --git a/cpp/monoprop/detail/evolution/layer_build/FusedApply.h b/cpp/monoprop/detail/evolution/layer_build/FusedApply.h index e2be7bae..4eb78ddf 100644 --- a/cpp/monoprop/detail/evolution/layer_build/FusedApply.h +++ b/cpp/monoprop/detail/evolution/layer_build/FusedApply.h @@ -19,7 +19,6 @@ #include "monoprop/TypeAliases.h" #include "monoprop/detail/evolution/CosineRecompute.h" #include "monoprop/detail/evolution/layer_build/Common.h" -#include "monoprop/picture/Picture.h" namespace monoprop::detail { diff --git a/cpp/monoprop/detail/graph/MPGraphViews.h b/cpp/monoprop/detail/graph/MPGraphViews.h index 3128a082..4fcbd16e 100644 --- a/cpp/monoprop/detail/graph/MPGraphViews.h +++ b/cpp/monoprop/detail/graph/MPGraphViews.h @@ -18,7 +18,6 @@ #include #include #include -#include #include #include @@ -33,6 +32,28 @@ class LayerIndexOutOfRange : public std::out_of_range { using std::out_of_range::out_of_range; }; +/// The optimizer slot of stored layer `layer_idx` in an `n`-layer graph. +// Its own inverse, so the slot-to-layer direction calls it too. The one spelling of MPGraph's storage +// invariant: every conversion between store order and optimizer order goes through it. It lives here, +// below MPGraph, so the views over a graph reach the same spelling the graph itself uses. +constexpr auto slot_of_layer(size_t layer_idx, size_t n) -> size_t { + return n - 1 - layer_idx; +} + +static_assert(slot_of_layer(0, 4) == 3); +static_assert(slot_of_layer(slot_of_layer(1, 4), 4) == 1); + +// The store offset of layer `layer_idx` in a `count`-layer store. `reverse` means the store runs against +// layer order, so layer 0 is its last element. A graph and the views over it both index through here, which +// is what keeps them agreeing on the mapping and on the diagnostic. +inline auto checked_layer_offset(size_t layer_idx, size_t count, bool reverse) -> size_t { + if (layer_idx >= count) { + throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, count)); + } + + return reverse ? slot_of_layer(layer_idx, count) : layer_idx; +} + // One rank's own graph memory only. struct GraphMemoryBreakdown final { size_t layer_descriptor_bytes = 0; @@ -75,19 +96,13 @@ class MPGraphView { auto layers() const -> size_t { return layers_.size(); } - auto get_layer(size_t layer_idx) const -> const Layer & { return layers_[checked_layer_offset(layer_idx)]; } + auto get_layer(size_t layer_idx) const -> const Layer & { + return layers_[checked_layer_offset(layer_idx, layers_.size(), reverse_)]; + } auto get_layer_traversal(size_t layer_idx) const -> LayerTraversal { return get_layer(layer_idx).traversal(); } private: - auto checked_layer_offset(size_t layer_idx) const -> size_t { - if (layer_idx >= layers_.size()) { - throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers_.size())); - } - - return reverse_ ? layers_.size() - 1 - layer_idx : layer_idx; - } - std::span layers_; bool reverse_ = false; }; diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 3be3a0f3..da647dce 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -348,7 +348,6 @@ auto MonomialPropagator::packed_inline_width_() const -> size_t { } template -template auto MonomialPropagator::apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD> { ++initial_operator_epoch_; @@ -373,7 +372,7 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o } } - return mp_op_.update_initial_operator(new_op, P::picture); + return mp_op_.update_initial_operator(new_op, picture_); } template @@ -535,13 +534,13 @@ auto MonomialPropagator::evolve_mode_build_graph_(const std::vector rot_len, size_t slot) { - this->template propagate_one_

(mono, - rot_len, - std::nullopt, - std::nullopt, - parameter_mapping[slot], - gen_coeffs[slot], - gate_indices[slot]); + this->propagate_one_(mono, + rot_len, + std::nullopt, + std::nullopt, + parameter_mapping[slot], + gen_coeffs[slot], + gate_indices[slot]); }); } @@ -552,10 +551,10 @@ auto MonomialPropagator::evolve_mode_graph_with_coeffs_(const std::vec const VecD &gen_coeffs, const VecZ &gate_indices, const VecD ¶meters, - const VecD &operator_coeffs, + VecD operator_coeffs, std::optional only_rotate_len_k) -> void { auto mapped_params = map_params(parameters, parameter_mapping, gen_coeffs, 1.0); - auto coeffs = operator_coeffs; + auto coeffs = std::move(operator_coeffs); run_gate_loop_

( majoranas, @@ -566,8 +565,7 @@ auto MonomialPropagator::evolve_mode_graph_with_coeffs_(const std::vec const auto [build_angle, apply_angle] = gate_angle_

(mapped_params, slot); // The cos word list is not persisted on the layer; the builder moves it out transiently. auto cos = std::make_shared(); - auto storage = - this->template build_evolve_result_

(mono, rot_len, std::cref(coeffs), build_angle, cos.get()); + auto storage = this->build_evolve_result_(mono, rot_len, std::cref(coeffs), build_angle, cos.get()); graph_.append(storage, parameter_mapping[slot], gen_coeffs[slot], gate_indices[slot]); this->template extend_coeffs_from_current_picture_if_needed_

(coeffs); @@ -599,17 +597,15 @@ auto MonomialPropagator::evolve_mode_contract_immediately_(const std:: CosMask cos; detail::FusedContract fc; bool fused_scale = false; - this->template build_evolve_result_

(mono, - rot_len, - std::cref(*op_coeffs), - build_angle, - &cos, - &fc, - op_coeffs, - &fused_scale); + this->build_evolve_result_(mono, + rot_len, + std::cref(*op_coeffs), + build_angle, + &cos, + &fc, + op_coeffs, + &fused_scale); this->template extend_coeffs_from_current_picture_if_needed_

(*op_coeffs); - // build_layer resolves P again from the same Picture value for its sink; the fused cosine - // sweep and this apply must agree. detail::apply_fused_contract

(fc, *op_coeffs, cos, apply_angle, fused_scale); }); } @@ -686,7 +682,7 @@ auto MonomialPropagator::build_graph(const std::vector &majorana gen_coeffs, local_gates, *parameters, - seed, + std::move(seed), only_rotate_len_k); }); } @@ -740,7 +736,6 @@ auto MonomialPropagator::run_gate_loop_(const std::vector &major } template -template auto MonomialPropagator::build_evolve_result_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs, @@ -766,14 +761,13 @@ auto MonomialPropagator::build_evolve_result_(const VecZ &gen_vec, comm_, out_cos, fused_contract, - P::picture, + picture_, fused_scale_coeffs, fused_scale, basis_); } template -template auto MonomialPropagator::propagate_one_(const VecZ &gen_vec, std::optional only_rotate_len_k, std::optional> coeffs, @@ -781,10 +775,7 @@ auto MonomialPropagator::propagate_one_(const VecZ &gen_vec, size_t param_index, double gen_coeff, size_t gate_index) -> void { - graph_.append(this->template build_evolve_result_

(gen_vec, only_rotate_len_k, coeffs, param), - param_index, - gen_coeff, - gate_index); + graph_.append(build_evolve_result_(gen_vec, only_rotate_len_k, coeffs, param), param_index, gen_coeff, gate_index); } template @@ -964,8 +955,7 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optional(pare_graph(graph_, keep, count, comm_, full_cos_of_layer, P::pare_sweep)); + graph = std::make_shared(pare_graph(graph_, keep, count, comm_, full_cos_of_layer)); } else { graph = std::shared_ptr(std::shared_ptr{}, &graph_); @@ -1079,9 +1069,10 @@ auto MonomialPropagator::contract_partially_(const VecD ¶meters, b // The pictures differ in three values only: the source vector, the (phase, reverse) map_params pair, // and the slot that receives an inplace result. Everything else -- and the order of every flop -- is shared. + // The phase is apply_sign because this replays the very angles the build applied. const VecD &source = P::live_coeffs(mp_op_); const auto mapped_params = - map_params(parameters, parameter_mapping, gen_coeffs, P::contract_phase, P::contract_reverse); + map_params(parameters, parameter_mapping, gen_coeffs, P::apply_sign, P::contract_reverse); VecD evolved = evolve_operator_with_recompute_(VecD(source), graph_.contraction_view(), mapped_params); if (inplace) { diff --git a/cpp/monoprop/detail/pare/PareGraph.cpp b/cpp/monoprop/detail/pare/PareGraph.cpp index 114aca01..030fb591 100644 --- a/cpp/monoprop/detail/pare/PareGraph.cpp +++ b/cpp/monoprop/detail/pare/PareGraph.cpp @@ -104,8 +104,7 @@ auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, mpi::Comm comm, - const std::function &full_cos_of_layer, - PareSweep sweep) -> MPGraph { + const std::function &full_cos_of_layer) -> MPGraph { const size_t num_layers = graph.layers(); const auto my_rank = static_cast(mpi::rank(comm)); @@ -124,9 +123,9 @@ auto pare_graph(const MPGraph &graph, // mark_cross_rank_endpoints_kept), so nodes_to_keep stays consistent across ranks with no exchange. // Cross-rank lists are never pruned; the keep-set only has to be right so cos pruning stays exact. for (size_t iter = 0; iter < num_layers; ++iter) { - // The sweep walks away from the seed end, so reachability reaches every kept node before the - // cosine filter runs on the layer that produced it. - const size_t layer_idx = (sweep == PareSweep::FromInput) ? iter : slot_of_layer(iter, num_layers); + // Unbuild order: away from the seed, so reachability reaches every kept node before the cosine + // filter runs on the layer that produced it. + const size_t layer_idx = graph.layer_of_unbuild_step(iter); const auto &layer = graph.get_layer(layer_idx); const auto lt = layer.traversal(); diff --git a/cpp/monoprop/picture/Picture.h b/cpp/monoprop/picture/Picture.h index 030f7a39..8c056b80 100644 --- a/cpp/monoprop/picture/Picture.h +++ b/cpp/monoprop/picture/Picture.h @@ -42,17 +42,18 @@ namespace monoprop { struct HeisenbergPicture { static constexpr Picture picture = Picture::Heisenberg; - static constexpr bool is_schrodinger = false; + static constexpr bool is_schrodinger = (picture == Picture::Schrodinger); // Simulation step i consumes optimizer slot n-1-i: the observable walks the circuit backwards. - static auto gate_slot(size_t i, size_t n) -> size_t { return n - 1 - i; } + static auto gate_slot(size_t i, size_t n) -> size_t { return slot_of_layer(i, n); } // gate_slot's slope, as MPGraph needs it. static constexpr ArrivalOrder arrival_order = ArrivalOrder::DescendingSlot; static constexpr double apply_sign = 1.0; // the applied angle is the build angle - // map_params() arguments for contract_partially: forward phase, written in reverse. - static constexpr double contract_phase = 1.0; - static constexpr bool contract_reverse = true; + // contract_partially replays the build's steps in build order, so its params are indexed by step while + // parameter_mapping is indexed by slot. map_params must therefore invert gate_slot, which is a reversal + // exactly when gate_slot is not the identity -- that is, exactly when the slots descend. + static constexpr bool contract_reverse = (arrival_order == ArrivalOrder::DescendingSlot); // The live vector the gates mutate, and the slot it lives in. template @@ -82,8 +83,6 @@ struct HeisenbergPicture { static auto pare_seed(const EvalState &state, const VecD & /*op*/, double threshold) -> std::pair { return {state.indices_above(threshold), state.length()}; } - // The state contracts against the evolved operator, so the seed sits at the replay's output end. - static constexpr PareSweep pare_sweep = PareSweep::FromOutput; // A perf hint, never a correctness constraint: overflow spills losslessly. The bound is already in // physical slots (CutoffEvaluator::max_slot_bound), so nothing to scale. NumModes is explicit @@ -102,7 +101,7 @@ struct HeisenbergPicture { struct SchrodingerPicture { static constexpr Picture picture = Picture::Schrodinger; - static constexpr bool is_schrodinger = true; + static constexpr bool is_schrodinger = (picture == Picture::Schrodinger); // Simulation step i consumes optimizer slot i: the state walks the circuit front-to-back. static auto gate_slot(size_t i, size_t /*n*/) -> size_t { return i; } @@ -110,8 +109,7 @@ struct SchrodingerPicture { static constexpr ArrivalOrder arrival_order = ArrivalOrder::AscendingSlot; static constexpr double apply_sign = -1.0; // the applied angle is the negated build angle - static constexpr double contract_phase = -1.0; - static constexpr bool contract_reverse = false; + static constexpr bool contract_reverse = (arrival_order == ArrivalOrder::DescendingSlot); // The dense state IS the live evolved vector here, so it is both the source and the slot. template @@ -137,8 +135,6 @@ struct SchrodingerPicture { static auto pare_seed(const EvalState & /*state*/, const VecD &op, double threshold) -> std::pair { return {indices_above(op, threshold), op.size()}; } - // `op` is what the replay evolves, so the seed sits at the replay's input end. - static constexpr PareSweep pare_sweep = PareSweep::FromInput; // The state's monomials come from generate_paired_op(), not from cutoff_fn_, so the cutoff carries // no structural bound on them. @@ -154,10 +150,8 @@ concept PicturePolicy = requires { { P::picture } -> std::convertible_to; { P::is_schrodinger } -> std::convertible_to; { P::apply_sign } -> std::convertible_to; - { P::contract_phase } -> std::convertible_to; { P::contract_reverse } -> std::convertible_to; { P::arrival_order } -> std::convertible_to; - { P::pare_sweep } -> std::convertible_to; }; static_assert(PicturePolicy); diff --git a/cpp/tests/GraphBuildHarness.h b/cpp/tests/GraphBuildHarness.h index 5602d5b0..db8fae21 100644 --- a/cpp/tests/GraphBuildHarness.h +++ b/cpp/tests/GraphBuildHarness.h @@ -20,8 +20,8 @@ #include "monoprop/MPGraph.h" -// gate_index is only a distinguishable tag for asserting slice/view ordering; the rest of a LayerCore -// is left empty. +// gate_index is only a distinguishable tag for asserting layer ordering; the rest of a LayerCore is left +// empty. namespace test_utils { inline auto core_with_gate(std::size_t gate_index) -> std::shared_ptr { diff --git a/cpp/tests/README.md b/cpp/tests/README.md index f167fe1f..b31c44ba 100644 --- a/cpp/tests/README.md +++ b/cpp/tests/README.md @@ -96,8 +96,9 @@ name and cannot address suite-nested cases, tests use flat coalescer, checked_* overflow guards, packed-phase storage + int8 read, build_layer_exchange_layout, and both arms of the D-from-B derivation). - **Graph / paring**: `pare_graph_tests.cpp`, `mpi_pare.cpp`, - `mp_graph_tests.cpp` (MPGraph slice_graph/slice_view transforms, the - front_offset lazy-compaction arms, MPGraphView reverse mapping + OOB throw). + `mp_graph_tests.cpp` (MPGraph layer indexing under either arrival order, + replay_view/contraction_view, replace_layer, clear, MPGraphView reverse + mapping + OOB throw). - **Transports / distribution**: `shm_comm_tests.cpp`, `hybrid_comm_tests.cpp` (MPI-only), `partition_equivalence_tests.cpp`, `mpi_distributed_layer_equivalence.cpp`, `mpi_fresh_insert_equivalence.cpp` diff --git a/cpp/tests/mp_graph_tests.cpp b/cpp/tests/mp_graph_tests.cpp index 4c5c7093..52b88160 100644 --- a/cpp/tests/mp_graph_tests.cpp +++ b/cpp/tests/mp_graph_tests.cpp @@ -24,7 +24,6 @@ #include "monoprop/MPGraph.h" using namespace monoprop; -using test_utils::core_with_gate; using test_utils::graph_with_gates; using test_utils::layer_with_gate; @@ -84,6 +83,21 @@ BOOST_AUTO_TEST_CASE(mp_graph_replace_layer_addresses_the_same_layer_as_get_laye } } +// The equality across arrival orders is the point, and it is what let pare_graph drop its sweep argument: +// unbuild step 0 is the last gate the build applied under either order, so the sweep needs no picture. +BOOST_AUTO_TEST_CASE(mp_graph_unbuild_step_is_newest_arrival_first_under_either_arrival) { + auto descending = graph_with_gates(ArrivalOrder::DescendingSlot, 4); + auto ascending = graph_with_gates(ArrivalOrder::AscendingSlot, 4); + for (std::size_t step = 0; step < 4; ++step) { + const auto d = descending.get_layer_traversal(descending.layer_of_unbuild_step(step)).gate_index(); + const auto a = ascending.get_layer_traversal(ascending.layer_of_unbuild_step(step)).gate_index(); + BOOST_CHECK_EQUAL(d, 3U - step); + BOOST_CHECK_EQUAL(a, 3U - step); + } + BOOST_CHECK_THROW((void)descending.layer_of_unbuild_step(4), std::out_of_range); + BOOST_CHECK_THROW((void)ascending.layer_of_unbuild_step(4), std::out_of_range); +} + BOOST_AUTO_TEST_CASE(mp_graph_clear_empties_and_leaves_the_graph_usable) { auto graph = graph_with_gates(ArrivalOrder::AscendingSlot, 5); graph.clear(); diff --git a/cpp/tests/pare_graph_tests.cpp b/cpp/tests/pare_graph_tests.cpp index b826cf07..c1654a0b 100644 --- a/cpp/tests/pare_graph_tests.cpp +++ b/cpp/tests/pare_graph_tests.cpp @@ -106,7 +106,7 @@ BOOST_AUTO_TEST_CASE(pare_graph_emits_expected_layer_kinds) { } // Heisenberg's sweep: the simulator built above is a Heisenberg one. - auto pared = pare_graph(graph, seed, local_index_count, MPI_COMM_SELF, provider, PareSweep::FromOutput); + auto pared = pare_graph(graph, seed, local_index_count, MPI_COMM_SELF, provider); BOOST_REQUIRE_EQUAL(pared.layers(), graph.layers()); size_t pruned_count = 0; From 2a636ab1fb5f0de31423dbbd95f980841096be5c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 14:25:56 +0200 Subject: [PATCH 12/13] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Roberto Di Remigio Eikås Signed-off-by: Roberto Di Remigio Eikås --- cpp/include/monoprop/MPFunctions.h | 2 -- cpp/monoprop/core/Picture.h | 2 -- 2 files changed, 4 deletions(-) diff --git a/cpp/include/monoprop/MPFunctions.h b/cpp/include/monoprop/MPFunctions.h index 8940683a..087fd0a8 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -107,8 +107,6 @@ monoprop_EXPORT auto ev_and_grad(const EvalRequest &request, const detail::CosCallbacks &cos = {}) -> std::pair; /// Prune `graph` to the subgraph reaching `nonzero_inds`; `full_cos_of_layer(i)` supplies layer i's full cosine set. -// Takes no sweep direction: the seed is the result of the whole evolution under either picture, so the -// sweep is always the graph's own build reversed. MPGraph::layer_of_unbuild_step is that order. monoprop_EXPORT auto pare_graph(const MPGraph &graph, const VecZ &nonzero_inds, size_t local_index_count, diff --git a/cpp/monoprop/core/Picture.h b/cpp/monoprop/core/Picture.h index 221bda74..a8edc812 100644 --- a/cpp/monoprop/core/Picture.h +++ b/cpp/monoprop/core/Picture.h @@ -14,8 +14,6 @@ #pragma once -// Deliberately dependency-free: MonomialPropagator.h carries a Picture in its signatures, and MPOperator.h -// in a parameter. The behaviour behind each value lives in picture/Picture.h. #include #include From 78f7e36aa5130ff9aae627369a809c504b26c195 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Roberto=20Di=20Remigio=20Eik=C3=A5s?= Date: Tue, 18 Aug 2026 14:27:16 +0200 Subject: [PATCH 13/13] chore: run clang-format --- cpp/monoprop/core/Picture.h | 1 - 1 file changed, 1 deletion(-) diff --git a/cpp/monoprop/core/Picture.h b/cpp/monoprop/core/Picture.h index a8edc812..ceb0bdfb 100644 --- a/cpp/monoprop/core/Picture.h +++ b/cpp/monoprop/core/Picture.h @@ -14,7 +14,6 @@ #pragma once - #include #include