refactor(c++)!: 🎨 de-tangle picture decisions

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<Heisenberg, Schrodinger>`, 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_`

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<CosMask(size_t)> &full_cos_of_layer) -> MPGraph;
 } // namespace monoprop

cpp/include/monoprop/MPGraph.h

@@ -15,8 +15,7 @@
 #pragma once
 
 #include <cstddef>
-#include <format>
-#include <stdexcept>
+#include <cstdint>
 #include <string_view>
 #include <utility>
 #include <vector>
@@ -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<Layer>::iterator;
-    using ConstLayerIterator = std::vector<Layer>::const_iterator;
-
-    bool schrodinger_;
+    ArrivalOrder arrival_;
     std::vector<Layer> 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<std::ptrdiff_t>(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<std::ptrdiff_t>(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<Layer> 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.

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 <size_t NumModes>
 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<unsigned int> schrodinger_cutoff,
+                       const PictureSpec &picture,
                        mpi::Comm comm,
                        std::optional<double> lower_atol = std::nullopt,
                        std::optional<double> 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 &parameters, double atol)
         -> std::vector<std::pair<VecZ, std::complex<double>>>;
 
-    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<MonomialList<NumModes>, 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<NumModes> cutoff_fn_;
     detail::MPOperator<NumModes> mp_op_;
     MPGraph graph_;
@@ -397,40 +398,53 @@ class MonomialPropagator {
     auto validate_cutoff_config_(CutoffType cutoff_type, const std::optional<std::vector<VecZ>> &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 <typename P>
+    auto initialize_operator_caches_() -> void;
 
+    // Grow `coeffs` to the operator's current term count, filling from the picture's live vector.
+    template <typename P>
     auto extend_coeffs_from_current_picture_if_needed_(VecD &coeffs) -> void;
 
+    template <typename P>
     auto evolve_mode_build_graph_(const std::vector<VecZ> &majoranas,
                                   const VecZ &parameter_mapping,
                                   const VecD &gen_coeffs,
                                   const VecZ &gate_indices,
                                   std::optional<size_t> 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<double, double> {
-        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 <typename P>
+    static auto gate_angle_(const VecD &mapped_params, size_t slot) -> std::pair<double, double> {
+        const double build_angle = mapped_params[slot];
+        return {build_angle, P::apply_sign * build_angle};
     }
 
+    template <typename P>
     auto evolve_mode_graph_with_coeffs_(const std::vector<VecZ> &majoranas,
                                         const VecZ &parameter_mapping,
                                         const VecD &gen_coeffs,
                                         const VecZ &gate_indices,
                                         const VecD &parameters,
-                                        const VecD &operator_coeffs,
+                                        VecD operator_coeffs, // by value: the caller's seed is dead after the call
                                         std::optional<size_t> only_rotate_len_k) -> void;
 
+    // build_layer resolves the same policy for its fused sink, so the cosine sweep and the apply agree.
+    template <typename P>
     auto evolve_mode_contract_immediately_(const std::vector<VecZ> &majoranas,
                                            const VecZ &parameter_mapping,
                                            const VecD &gen_coeffs,
                                            const VecD &parameters,
                                            std::optional<size_t> only_rotate_len_k) -> void;
 
-    template <typename EvolutionFunc>
+    // 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 <typename P, typename EvolutionFunc>
     auto run_gate_loop_(const std::vector<VecZ> &majoranas,
                         std::optional<size_t> 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<LayerCore>;
 
-    template <typename Fn,
+    template <typename P>
+    auto contract_partially_(const VecD &parameters, bool inplace) -> VecD;
+
+    template <typename P,
+              typename Fn,
               typename R = std::invoke_result_t<Fn, const EvalRequest &, mpi::Comm, const detail::CosCallbacks &>>
     auto make_functional_(Fn &&func, std::optional<double> pare_threshold) -> std::function<R(const VecD &)>;
 

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

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 &parameters,
                         const VecZ &parameter_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 &parameter_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_,