diff --git a/AGENTS.md b/AGENTS.md index a932c4ab..114efd52 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -67,6 +67,23 @@ 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, 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 + `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 `ArrivalOrder` naming the slope of the optimizer slots a build hands to `append()`, and each + 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 c7419ffe..087fd0a8 100644 --- a/cpp/include/monoprop/MPFunctions.h +++ b/cpp/include/monoprop/MPFunctions.h @@ -110,7 +110,6 @@ 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, 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..4d74e6af 100644 --- a/cpp/include/monoprop/MPGraph.h +++ b/cpp/include/monoprop/MPGraph.h @@ -15,8 +15,7 @@ #pragma once #include -#include -#include +#include #include #include #include @@ -27,49 +26,37 @@ namespace monoprop { +/// 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 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 +}; + /// 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 +// 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: - using LayerIterator = std::vector::iterator; - using ConstLayerIterator = std::vector::const_iterator; - - bool schrodinger_; + ArrivalOrder arrival_; 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(); } + // 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 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 schrodinger_ ? active_begin_iterator() : active_end_iterator(); } - - auto append_layer(Layer layer) -> void { layers_.emplace(append_position(), std::move(layer)); } - - auto checked_layer_offset(size_t layer_idx) const -> size_t { - if (layer_idx >= layers()) { - throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, layers())); - } - return active_begin_index() + layer_idx; + auto stored_offset_(size_t layer_idx) const -> size_t { + return checked_layer_offset(layer_idx, layers_.size(), reverse_indexing_()); } public: - explicit MPGraph(bool schrodinger) : schrodinger_(schrodinger) {} - - explicit MPGraph(bool schrodinger, std::vector layers) - : schrodinger_(schrodinger), - 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. @@ -80,27 +67,45 @@ 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_.emplace_back(std::move(storage)); } - /// Slice the graph at `key` (the number of earliest operations to include); `contract` also removes - /// the sliced part from this graph. - auto slice_graph(size_t key, bool contract = false) -> MPGraph; + /// 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); } - auto slice_view(size_t key) const -> MPGraphView; + /// 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 active_end_index() - active_begin_index(); } + auto layers() const -> size_t { return layers_.size(); } - auto get_layer(size_t layer_idx) -> Layer& { return layers_[checked_layer_offset(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& { return layers_[checked_layer_offset(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 active layers, in build order. - auto replay_view() const -> MPGraphView { return {layers_, active_begin_index(), layers(), false}; } + /// 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_()}; } - auto is_schrodinger() const -> bool { return schrodinger_; } + /// The layers in the order this graph's own build walked them, which is the order a contraction must + /// replay them in. + // 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/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index f80f3e87..a32e0c3b 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,7 @@ class MonomialPropagator { }); } - auto schrodinger() const -> bool { return schrodinger_; } + auto schrodinger() const -> bool { return picture_ == Picture::Schrodinger; } auto basis() const -> Basis { return basis_; } @@ -274,7 +275,7 @@ 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 { apply_initial_operator_(op_dict); } protected: static inline const auto ev_fn = [](const EvalRequest &request, @@ -290,8 +291,8 @@ class MonomialPropagator { /// so caches can refresh. 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,40 +398,53 @@ class MonomialPropagator { auto validate_cutoff_config_(CutoffType cutoff_type, const std::optional> &basis_change) const -> void; - auto initialize_operator_caches_() -> void; + // 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. - 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, 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. + 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; @@ -454,7 +468,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/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/MPGraph.cpp b/cpp/monoprop/MPGraph.cpp index a425582d..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,48 +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 (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]); - } - - if (contract && k != 0) { - layers_.resize(active_end - 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(schrodinger_, std::move(sliced_layers)); -} - -auto MPGraph::slice_view(size_t key) const -> MPGraphView { - const auto k = std::min(key, layers()); - if (schrodinger_) { - return {layers_, active_end_index() - k, k, true}; - } - return {layers_, active_begin_index(), k, false}; -} - 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; } @@ -108,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/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..ceb0bdfb --- /dev/null +++ b/cpp/monoprop/core/Picture.h @@ -0,0 +1,39 @@ +// 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 + +#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..4eb78ddf 100644 --- a/cpp/monoprop/detail/evolution/layer_build/FusedApply.h +++ b/cpp/monoprop/detail/evolution/layer_build/FusedApply.h @@ -29,15 +29,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/graph/MPGraphViews.h b/cpp/monoprop/detail/graph/MPGraphViews.h index 0af7d7a2..4fcbd16e 100644 --- a/cpp/monoprop/detail/graph/MPGraphViews.h +++ b/cpp/monoprop/detail/graph/MPGraphViews.h @@ -15,9 +15,9 @@ #pragma once #include +#include #include #include -#include #include #include @@ -32,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; @@ -56,34 +78,32 @@ 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, 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 >= count_) { - throw LayerIndexOutOfRange(std::format("Layer {} is out of range (layers={})", layer_idx, count_)); - } - - return base_ + (reverse_ ? count_ - 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/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 117590ec..da647dce 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_(with_picture(picture_, []() { return P::arrival_order; })), 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,17 +343,8 @@ 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 @@ -373,7 +372,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, picture_); } template @@ -496,33 +495,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,87 +523,91 @@ 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->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, 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; - 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); - }); + auto coeffs = std::move(operator_coeffs); + + 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->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->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); + detail::apply_fused_contract

(fc, *op_coeffs, cos, apply_angle, fused_scale); + }); } template @@ -650,10 +643,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 +672,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, + std::move(seed), + only_rotate_len_k); + }); } template @@ -714,22 +712,27 @@ 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 @@ -758,7 +761,7 @@ auto MonomialPropagator::build_evolve_result_(const VecZ &gen_vec, comm_, out_cos, fused_contract, - schrodinger_, + picture_, fused_scale_coeffs, fused_scale, basis_); @@ -833,7 +836,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) { @@ -859,10 +862,9 @@ auto MonomialPropagator::graph_gate_arrays_() const -> std::pair &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 +932,7 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optionalcore_term(); const auto comm = comm_; @@ -958,10 +954,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, comm_, full_cos_of_layer)); } else { graph = std::shared_ptr(std::shared_ptr{}, &graph_); @@ -1013,7 +1007,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 +1023,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 +1051,37 @@ 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(); - // 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; - } + // 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::apply_sign, P::contract_reverse); - 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 = evolve_operator_with_recompute_(VecD(source), graph_.contraction_view(), mapped_params); 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; - } - else { - evolved_op = evolve_operator_with_recompute_(VecD(op), graph_.slice_view(num_majoranas), 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; + graph_.clear(); } - 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..030fb591 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, - bool schrodinger, mpi::Comm comm, const std::function &full_cos_of_layer) -> MPGraph { const size_t num_layers = graph.layers(); @@ -116,13 +115,17 @@ 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. // 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); + // 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(); @@ -132,10 +135,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.is_schrodinger(), std::move(layers)); + return pared; } } // 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..8c056b80 --- /dev/null +++ b/cpp/monoprop/picture/Picture.h @@ -0,0 +1,171 @@ +// 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/MPGraph.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 = (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 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 + + // 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 + 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 = (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; } + // 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 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 + 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_reverse } -> std::convertible_to; + { P::arrival_order } -> 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..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 { @@ -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(bool schrodinger, std::size_t n) -> monoprop::MPGraph { - monoprop::MPGraph graph(schrodinger); +// 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/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/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..52b88160 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 @@ -24,96 +24,87 @@ #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; -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 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); -} - -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 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(/*schrodinger=*/false, 3); - auto sliced = graph.slice_graph(100, /*contract=*/false); - BOOST_CHECK_EQUAL(sliced.layers(), 3U); +// 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); } -// 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); - (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); +// 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(descending_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(ascending_view.get_layer_traversal(i).gate_index(), 3U - i); + } } -BOOST_AUTO_TEST_CASE(mp_graph_contract_noop_arm_keeps_dead_prefix_lazy) { - auto graph = graph_with_gates(/*schrodinger=*/false, 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); +// 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_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(descending_view.get_layer_traversal(i).gate_index(), i); + BOOST_CHECK_EQUAL(ascending_view.get_layer_traversal(i).gate_index(), i); + } } -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 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); +// 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_slice_view_heisenberg_forward_window) { - auto graph = graph_with_gates(/*schrodinger=*/false, 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); +// 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_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 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_clear_empties_and_leaves_the_graph_usable) { + 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); + BOOST_REQUIRE_EQUAL(graph.layers(), 1U); + BOOST_CHECK_EQUAL(graph.get_layer_traversal(0).gate_index(), 7U); } BOOST_AUTO_TEST_CASE(mp_graph_view_reverse_flag_flips_index_mapping) { @@ -122,8 +113,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); @@ -133,7 +124,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(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. 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..c1654a0b 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, /*schrodinger=*/false, 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); 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,