diff --git a/AGENTS.md b/AGENTS.md index 35bcd09e..35a555e6 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -71,6 +71,20 @@ Key files: - `cpp/include/monoprop/MonomialPropagator.h`: the single templated C++ engine `MonomialPropagator` (the Majorana/Pauli choice is a runtime `Basis`, not a separate class). Its `only_rotate_len_k` arguments use `std::optional`; `std::nullopt` means no gate-application length cap. +- `cpp/include/monoprop/Functional.h`: the two functional objects, + `ExpectationValueFunctional` and `ExpectationValueAndGradientFunctional`. Both derive + from `detail::FunctionalHandle` — the shared handle half — and each holds a + `detail::FunctionalPlan`: the propagator snapshot a call replays, plus the checks that say the + snapshot is still that propagator's. The plan's `std::variant` carries the single-partition shape and the + facade shape, so both paths have one public type, and the plan holds the snapshot rather than the + choice of what to compute, so one plan type backs either kind (each factory call builds its own). A + functional borrows from its propagator (the inverted index always, the graph unless pared), so it + must not outlive it; the bindings pin that with + `nb::keep_alive<0, 1>`. It does **not** snapshot the initial-operator weights: it reads the + `detail::OperatorWeights` set the propagator has published, so it follows an + `update_initial_operator` instead of going stale — the one exception being a Schrödinger plan with a + `pare_threshold`, whose keep-set came from the coefficients the re-weight replaced, which throws + (`follows_weights` reports which case an object is). - `src/monoprop/bindings/binder.h`: hand-written binding template; `tools/generate-*.py` generate the per-mode-width `bindings.cpp` and `_dispatch.py` from it (do not hand-edit the generated files). Both generators take the 32-mode storage-block rule from `tools/_binding_layout.py` — they must @@ -98,6 +112,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`. +- **`detail::FunctionalControl`** (`cpp/monoprop/detail/functional/Control.h`): the validity block a + propagator shares with every functional plan it makes — a structure revision, an alive flag, and the + name of the last structural change. A plan borrows from its propagator, so this is how it answers "is + the propagator still there, and does it still hold what I replay?" without dereferencing it. **Every + new mutating method must call `bump_structure_("its_name()")`** once the mutation has committed (a + rejected mutation must not bump); the settings that only gate the next build — the atols, the cutoff, + the cutoff type, the basis change — deliberately do not. Two halves go with that success-path bump: + guard the mutation itself with `bump_structure_on_unwind_` (`detail::BumpOnUnwind`), constructed after + the last rejection check and before the first write, so a mutator that throws part-way still + invalidates; and do not bump a fan-out whose children all no-op — decide with the facade-transparent + readers *before* fanning out, or `monoprop_PARTITIONS=auto` invalidates where `off` does not. A plan + additionally re-derives the operator's store pointer and inverted-index row count as a backstop, so a + missing bump reports staleness instead of folding a rebuilt index. The block also carries the + published `OperatorWeights`: a re-weight publishes a new set rather than bumping, which is what lets a + live functional follow it, and bumps only if it fails part-way. Publishing runs on the propagator's own thread (a facade publishes through + `for_each_partition_`), and a plan reads the set once per call so `op` and `core_term` cannot come from + two publications. - **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_` @@ -137,6 +168,13 @@ mp = MajoranaPropagator(operator, initial_state, cutoff=4) ### Testing Structure +- `cpp/tests/functional_validity.cpp` is the **mutation table**: one row per public mutating method of + `MonomialPropagator`, recording what a functional built *before* that mutator ran does when called + *after* it — throw, or answer from its own snapshot. `MonomialPropagator::num_mutating_methods` pins + the roster and the table `static_assert`s against it, so adding a mutator means bumping that + constant and adding a row (the build fails until you do). + `tests/test_parameter_validation.py::TestFunctionalValidityTable` mirrors the same rows through the + Python front end, over `monoprop_PARTITIONS=off` and `=auto`. - `tests/cases.py`: Parametrized test cases using `pytest-cases`; `load_problem()` loads a `tests/data/*.msgpack` fixture directly into the public API (`MonomialCircuit` + `MonomialOperator`). C++ tests use the equivalent `test_utils::load_case()` in `cpp/tests/TestData.h` - Fixture msgpack schema is documented in `tests/data/README.md` - Tests validate against exact solutions for small systems @@ -167,6 +205,10 @@ mp = MajoranaPropagator(operator, initial_state, cutoff=4) 7. Add Python bindings in `src/monoprop/bindings/binder.h` 8. Regenerate bindings with `tools/generate-binders.py` 9. Test with both C++ and Python tests +10. If the new method mutates a `MonomialPropagator`: call `bump_structure_` from it and guard it with + `bump_structure_on_unwind_` (see `detail::FunctionalControl`), bump + `MonomialPropagator::num_mutating_methods`, and add its row to the mutation table (see "Testing + Structure"). A mutator with no row leaves its effect on a live functional unrecorded. ## Documentation Maintenance Policy diff --git a/cpp/include/monoprop/CMakeLists.txt b/cpp/include/monoprop/CMakeLists.txt index 98159824..9fb002ea 100644 --- a/cpp/include/monoprop/CMakeLists.txt +++ b/cpp/include/monoprop/CMakeLists.txt @@ -16,6 +16,7 @@ target_sources( "${PROJECT_BINARY_DIR}/include/${PROJECT_NAME}/Info.h" "${PROJECT_BINARY_DIR}/include/${PROJECT_NAME}/Variants.h" "${PROJECT_SOURCE_DIR}/cpp/include/${PROJECT_NAME}/Evolution.h" + "${PROJECT_SOURCE_DIR}/cpp/include/${PROJECT_NAME}/Functional.h" "${PROJECT_SOURCE_DIR}/cpp/include/${PROJECT_NAME}/MPFunctions.h" "${PROJECT_SOURCE_DIR}/cpp/include/${PROJECT_NAME}/MPGraph.h" "${PROJECT_SOURCE_DIR}/cpp/include/${PROJECT_NAME}/MonomialPropagator.h" diff --git a/cpp/include/monoprop/Functional.h b/cpp/include/monoprop/Functional.h new file mode 100644 index 00000000..bef5ceca --- /dev/null +++ b/cpp/include/monoprop/Functional.h @@ -0,0 +1,217 @@ +// 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 +#include +#include +#include +#include + +#include "monoprop/MPFunctions.h" +#include "monoprop/MPGraph.h" +#include "monoprop/TypeAliases.h" +#include "monoprop/Validation.h" +#include "monoprop/detail/evolution/CosineRecomputeCallbacks.h" +#include "monoprop/detail/functional/Control.h" +#include "monoprop/detail/mpi/Comm.h" +#include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/detail/partition/PartitionGroup.h" + +namespace monoprop { + +template +class MonomialPropagator; + +namespace detail { + +/// Immutable, shared functional replay plan. +/// Borrowed fields require the functional to outlive neither its propagator nor its partitions. +template +class FunctionalPlan { +public: + /// A single-partition propagator's snapshot: one replay of its graph against its operator. + struct Local { + // Build-time weights; also used until a re-weight publishes new weights. + std::shared_ptr weights; + // Owned snapshot: operator rows can grow, but `op` cannot. + EvalState state; ///< the contraction partner, sparse (Heisenberg) or dense (Schrodinger) + VecZ parameter_mapping; ///< optimizer order: which parameter drives graph layer i + VecD gen_coeffs; ///< optimizer order, parallel to parameter_mapping + // Owned because `cos` holds raw pointers into graph layers. + std::shared_ptr graph; + // `cos` borrows columns from the propagator's inverted index. + CosCallbacks cos; + mpi::Comm comm{}; ///< real MPI across nodes, or the in-process comm across partitions + + // Borrowed operator and inverted-index identity check. + // This catches a missing revision bump before the stale index is used. + const MPOperator *mp_op{nullptr}; + const OperatorIndex *op_store{nullptr}; + size_t inverted_index_rows{0}; + + // A coefficient-pared Schrodinger graph cannot follow a re-weight. + bool pared_from_operator{false}; + }; + + /// A partition facade's snapshot: one child plan per partition, replayed together. + struct Fanout { + // Borrowed from the facade propagator. + partition::PartitionGroup *group{nullptr}; + std::vector> partitions; ///< in partition order + }; + + /// Pins the propagator control block and its current revision. + FunctionalPlan(size_t num_params, std::shared_ptr control, Local local) + : num_params_(num_params), + control_(std::move(control)), + expected_revision_(control_->structure_revision.load()), + shape_(std::move(local)) {} + + FunctionalPlan(size_t num_params, std::shared_ptr control, Fanout fanout) + : num_params_(num_params), + control_(std::move(control)), + expected_revision_(control_->structure_revision.load()), + shape_(std::move(fanout)) {} + + /// Required parameter-axis length. + auto num_params() const -> size_t { return num_params_; } + + /// Whether calls may follow re-weighted coefficients. + auto follows_weights() const -> bool { + if (const auto *fanout = std::get_if(&shape_)) { + // Child plans share picture and threshold. + return fanout->partitions.front()->follows_weights(); + } + return !std::get(shape_).pared_from_operator; + } + + /// Throw unless `params` and the propagator still match this plan. + // The facade validates its group; child plans validate on their partition masters. + auto validate(const VecD ¶ms) const -> void { + // Check liveness before reading the borrowed operator for the layout check. + const bool alive = control_->propagator_alive.load(); + const auto *local = alive ? std::get_if(&shape_) : nullptr; + validate_functional_state({.propagator_alive = alive, + .current_revision = control_->structure_revision.load(), + .expected_revision = expected_revision_, + .operator_layout_unchanged = local == nullptr || operator_layout_unchanged(*local), + .last_structural_change = control_->last_structural_change.load()}); + validate_functional_call(params, num_params_); + } + + /// Replay locally, or return partition 0's facade result. + template > + auto evaluate(Fn &&fn, const VecD ¶ms) const -> R { + validate(params); + if (const auto *fanout = std::get_if(&shape_)) { + // Every partition must join its synchronized collective; partition 0 has the result. + return std::move(partition::collect_on_all(*fanout->group, [&](int r) -> R { + return fanout->partitions[static_cast(r)]->evaluate(fn, params); + })[0]); + } + const auto &local = std::get(shape_); + // Keeps `request.op` alive for the call. + const auto weights = resolve_weights(local); + return fn(EvalRequest{.e_core = weights->core_term, + .state = local.state, + .op = weights->op, + .parameter_mapping = local.parameter_mapping, + .gen_coeffs = local.gen_coeffs, + .graph = local.graph->replay_view(), + .params = params}, + local.comm, + local.cos); + } + +private: + // Load matching `op` and `core_term` from the current weight publication. + auto resolve_weights(const Local &local) const -> std::shared_ptr { + auto published = control_->weights.load(); + // No publication or no re-weight: build-time weights are current. + if (published == nullptr || published == local.weights) { + return local.weights; + } + validate_weight_refresh({.weights_revision = published->structure_revision, + .expected_revision = expected_revision_, + .may_follow_weights = follows_weights()}); + return published; + } + + // Do not call inverted_index(): it could rebuild the borrowed index. + static auto operator_layout_unchanged(const Local &local) -> bool { + return local.mp_op->store.get() == local.op_store && local.mp_op->inverted_index_.has_value() + && local.mp_op->inverted_index_->rows() == local.inverted_index_rows; + } + + size_t num_params_{0}; + std::shared_ptr control_; + size_t expected_revision_{0}; + std::variant shape_; +}; + +/// Shared functional handle; derived types differ only in `operator()`. +template +class FunctionalHandle { +public: + /// Required parameter-axis length. + auto num_params() const -> size_t { return plan_->num_params(); } + + /// Whether calls may follow a re-weight. + auto follows_weights() const -> bool { return plan_->follows_weights(); } + +protected: + explicit FunctionalHandle(std::shared_ptr> plan) : plan_(std::move(plan)) {} + + std::shared_ptr> plan_; +}; + +} // namespace detail + +/// Reusable expectation value: `fn(parameters) -> double`. +/// +/// Borrows its propagator and throws after structural mutation. It follows re-weighted initial-operator +/// coefficients unless its graph was coefficient-pared. +template +class ExpectationValueFunctional : public detail::FunctionalHandle { +public: + auto operator()(const VecD ¶meters) const -> double { return this->plan_->evaluate(ev, parameters); } + +private: + friend class MonomialPropagator; + + explicit ExpectationValueFunctional(std::shared_ptr> plan) + : detail::FunctionalHandle(std::move(plan)) {} +}; + +/// Expectation value and gradient from one backward pass. +/// `fn(parameters) -> (value, gradient)`, with the gradient in parameter-axis order. +template +class ExpectationValueAndGradientFunctional : public detail::FunctionalHandle { +public: + auto operator()(const VecD ¶meters) const -> std::pair { + return this->plan_->evaluate(ev_and_grad, parameters); + } + +private: + friend class MonomialPropagator; + + explicit ExpectationValueAndGradientFunctional(std::shared_ptr> plan) + : detail::FunctionalHandle(std::move(plan)) {} +}; + +} // namespace monoprop diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index f80f3e87..7adf5fe0 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -34,6 +34,7 @@ #include #include "monoprop/Evolution.h" +#include "monoprop/Functional.h" #include "monoprop/MPFunctions.h" #include "monoprop/MPGraph.h" #include "monoprop/TypeAliases.h" @@ -94,6 +95,12 @@ class MonomialPropagator { static constexpr auto num_modes{NumModes}; static constexpr auto storage_num_modes{NumModes}; + /// How many public methods mutate this propagator, and so may invalidate a live functional. + // cpp/tests/functional_validity.cpp static_asserts that its table has a row for each, so bumping + // this when adding a mutator breaks that build until the new method's effect on a functional is + // recorded. That table is the roster; keeping the names here too would be a copy nothing enforces. + static constexpr size_t num_mutating_methods{10}; + auto logical_num_modes() const -> size_t { return logical_num_modes_; } /// Term count on this rank (allreduce for global). @@ -255,11 +262,12 @@ class MonomialPropagator { auto expectation_value_and_gradient(const VecD ¶meters) -> std::pair; /// `pare_threshold` is the edge-retention cutoff for a masked plan; nullopt keeps the exact graph. + /// The result borrows from this propagator, so it must not outlive it. auto expectation_value_functional(std::optional pare_threshold = std::nullopt) - -> std::function; + -> ExpectationValueFunctional; auto expectation_value_and_gradient_functional(std::optional pare_threshold = std::nullopt) - -> std::function(const VecD &)>; + -> ExpectationValueAndGradientFunctional; /// Contract the graph into the operator (Heisenberg) or state (Schrodinger). `inplace` consumes the /// graph and updates internal state; otherwise nothing is mutated. Core term excluded either way. @@ -274,18 +282,11 @@ class MonomialPropagator { auto evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>>; + /// Re-weight the initial operator in place; every existing term op_dict omits is zeroed, the + /// identity term included. virtual auto update_initial_operator(const OperatorDict &op_dict) -> void { apply_initial_operator_(op_dict); } protected: - static inline const auto ev_fn = [](const EvalRequest &request, - mpi::Comm comm, - const detail::CosCallbacks &cos) -> double { return ev(request, comm, cos); }; - - static inline const auto ev_and_grad_fn = - [](const EvalRequest &request, mpi::Comm comm, const detail::CosCallbacks &cos) -> std::pair { - return ev_and_grad(request, comm, cos); - }; - /// Distribute op_dict across ranks and apply this rank's share; returns its new (terms, coeffs) /// so caches can refresh. auto apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD>; @@ -308,10 +309,6 @@ class MonomialPropagator { std::optional lower_atol_, upper_atol_; double core_term_{0.0}; - // Bumped by every initial-operator re-weight. A functional snapshots the operator coefficients, so - // it captures this and rejects a later call once it moves, as it does for a rebuilt graph. - size_t initial_operator_epoch_{0}; - size_t logical_num_modes_{NumModes}; CutoffType cutoff_type_; @@ -326,6 +323,35 @@ class MonomialPropagator { // PartitionGroup rebinds a cloned partition's comm_ to its own transport during a deep copy. friend class detail::partition::PartitionGroup; + // Shared with every functional plan this propagator makes: they read it, only this propagator writes + // it. A copy gets its own block, because a copy carries no functionals. + std::shared_ptr functional_control_{std::make_shared()}; + + // Record that what a plan replays has moved. `site` must be a string literal: plans read it after + // this propagator is gone. A mutation that is rejected outright must not bump -- it changed nothing, + // so it must not invalidate anything. A mutation that can fail part-way bumps on the failure path + // instead: invalidating a functional that did not need it costs a rebuild, answering from a snapshot + // of a half-written operator costs a wrong number. + auto bump_structure_(const char *site) const -> void { functional_control_->bump(site); } + + // The failure-path half of the rule above: guard the mutation itself, so a mutator that throws + // after committing part of its work still invalidates. `armed` is false where the guarded call + // is a known no-op. See detail::BumpOnUnwind. + auto bump_structure_on_unwind_(const char *site, bool armed = true) const -> detail::BumpOnUnwind { + return detail::BumpOnUnwind(*functional_control_, site, armed); + } + + // Hand the current initial-operator weights to every functional over this propagator, as one + // immutable set stamped with the revision it belongs to. Called on a re-weight, which is what makes a + // functional follow it, and never on a facade -- a facade holds no terms, and its partitions publish + // their own on their own masters. + auto publish_weights_() -> std::shared_ptr; + + // The weight set a new plan is built over: the published one when it still belongs to this revision, + // a fresh publication otherwise. Reusing matters -- a plan detects a re-weight by comparing pointers, + // so republishing an identical set would read as one. + auto weights_for_plan_() -> std::shared_ptr; + // A facade's own graph_/mp_op_ are never populated, so handing them out would return plausible-looking // empty state; there is no meaningful merge either, since the callers want one partition's raw layout. auto require_single_partition_(const char *what) const -> void { @@ -454,9 +480,9 @@ class MonomialPropagator { VecD *fused_scale_coeffs = nullptr, bool *fused_scale = nullptr) -> std::shared_ptr; - template > - auto make_functional_(Fn &&func, std::optional pare_threshold) -> std::function; + // One plan serves both functional kinds: it holds the snapshot and the validity checks, not the + // choice of what to compute. On a facade it holds one child plan per partition. + auto make_plan_(std::optional pare_threshold) -> std::shared_ptr>; // Reconstruct the optimizer-order (parameter_mapping, gen_coeffs) arrays from the layers' gate info. auto graph_gate_arrays_() const -> std::pair; diff --git a/cpp/monoprop/Validation.cpp b/cpp/monoprop/Validation.cpp index 46e5252d..aa7cde97 100644 --- a/cpp/monoprop/Validation.cpp +++ b/cpp/monoprop/Validation.cpp @@ -20,21 +20,19 @@ namespace monoprop { -// Every validate_* precondition on a caller-supplied argument reports through this one type. +// Invalid caller input reports through this type. class ValidationError : public std::runtime_error { public: using std::runtime_error::runtime_error; }; -// The propagator was mutated after a functional captured what it replays: a rebuilt graph leaves the -// functional's parameter mapping describing a graph that is gone, a re-weight leaves its snapshotted -// operator coefficients stale. +// A functional no longer matches its propagator. class StaleFunctionalGraph : public std::runtime_error { public: using std::runtime_error::runtime_error; }; -// Declared in Validation.h and used across translation units, so internal linkage does not apply. +// Exported validators below have external linkage. // NOLINTBEGIN(misc-use-internal-linkage) namespace { @@ -99,20 +97,37 @@ auto validate_functional_call(const VecD ¶meters, size_t expected_num_params } } -auto validate_expected_graph_layers(size_t current_layers, size_t expected_layers) -> void { - if (current_layers != expected_layers) { - throw StaleFunctionalGraph(std::format("MP object has been modified since the functional was created. " - "Previous number of graph layers was {} and now is {}.", - expected_layers, - current_layers)); +auto validate_functional_state(const FunctionalState &state) -> void { + if (!state.propagator_alive) { + throw StaleFunctionalGraph("The propagator this functional was built from has been destroyed. " + "A functional reads the propagator's operator index directly, so it " + "cannot outlive it; keep the propagator alive, or build the functional " + "again from a live one."); + } + if (state.current_revision != state.expected_revision) { + throw StaleFunctionalGraph(std::format( + "MP object has been modified since the functional was created: {} changed the " + "graph or operator the functional replays. Create a new functional.", + state.last_structural_change != nullptr ? state.last_structural_change : "a structural mutation")); + } + if (!state.operator_layout_unchanged) { + throw StaleFunctionalGraph("MP object has been modified since the functional was created: the " + "operator index it reads was rebuilt or grown. Create a new functional."); } } -auto validate_expected_initial_operator(size_t current_epoch, size_t expected_epoch) -> void { - if (current_epoch != expected_epoch) { - throw StaleFunctionalGraph("MP object has been modified since the functional was created. " - "The initial operator was re-weighted, so the coefficients the functional " - "snapshotted are stale; create a new functional."); +auto validate_weight_refresh(const WeightRefresh &refresh) -> void { + if (!refresh.may_follow_weights) { + throw StaleFunctionalGraph("MP object has been modified since the functional was created: the " + "initial operator was re-weighted, and this functional pares its graph " + "against the operator coefficients (a Schrodinger picture functional " + "with a pare_threshold), so it cannot follow the new weights. Create a " + "new functional."); + } + if (refresh.weights_revision != refresh.expected_revision) { + throw StaleFunctionalGraph("MP object has been modified since the functional was created: the " + "initial-operator weights it reads belong to a different graph. Create " + "a new functional."); } } diff --git a/cpp/monoprop/Validation.h b/cpp/monoprop/Validation.h index db93e622..7e7cf941 100644 --- a/cpp/monoprop/Validation.h +++ b/cpp/monoprop/Validation.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include @@ -22,15 +23,11 @@ namespace monoprop { -// Each of these throws when the stated condition does not hold: ValidationError for inconsistent -// arguments, StaleFunctionalGraph when the propagator was mutated after a functional captured the -// graph and operator it replays. Both derive from std::runtime_error, so catching that still catches -// either. +// Validators throw ValidationError for invalid arguments or StaleFunctionalGraph for stale plans. monoprop_EXPORT auto validate_coefficient_lengths(const VecZ ¶meter_mapping, const VecD &gen_coeffs) -> void; -// gate_indices records which ingested gate each monomial came from: one entry per monomial, forming -// contiguous runs from 0. +// gate_indices has one contiguous-run index per ingested monomial. monoprop_EXPORT auto validate_gate_indices(const VecZ &gate_indices, size_t num_monomials) -> void; // params must have max(parameter_mapping)+1 entries. @@ -38,11 +35,25 @@ monoprop_EXPORT auto validate_parameters_length(const VecD ¶ms, const VecZ & monoprop_EXPORT auto validate_functional_call(const VecD ¶meters, size_t expected_num_params) -> void; -// The graph must still have the layer count the functional was built against. -monoprop_EXPORT auto validate_expected_graph_layers(size_t current_layers, size_t expected_layers) -> void; - -// The initial operator must not have been re-weighted since the functional snapshotted its coefficients. -monoprop_EXPORT auto validate_expected_initial_operator(size_t current_epoch, size_t expected_epoch) -> void; +/// Propagator state required before reading a functional's borrowed data. +struct FunctionalState { + bool propagator_alive; ///< False after destruction begins. + size_t current_revision; ///< Current structure revision. + size_t expected_revision; ///< Build-time structure revision. + bool operator_layout_unchanged; ///< Borrowed index still matches its store and rows. + const char *last_structural_change; ///< Last mutation, or nullptr. +}; + +monoprop_EXPORT auto validate_functional_state(const FunctionalState &state) -> void; + +/// Inputs for checking whether a functional may follow new weights. +struct WeightRefresh { + size_t weights_revision; ///< Revision at publication. + size_t expected_revision; ///< Build-time revision. + bool may_follow_weights; ///< Functional weight-following policy. +}; + +monoprop_EXPORT auto validate_weight_refresh(const WeightRefresh &refresh) -> void; // only_rotate_len_k is optional; when set it must satisfy 0 < k <= max_k. monoprop_EXPORT auto validate_only_rotate_len_k_(std::optional only_rotate_len_k, size_t max_k) -> void; diff --git a/cpp/monoprop/detail/CMakeLists.txt b/cpp/monoprop/detail/CMakeLists.txt index f75f51c5..d075e550 100644 --- a/cpp/monoprop/detail/CMakeLists.txt +++ b/cpp/monoprop/detail/CMakeLists.txt @@ -8,6 +8,7 @@ target_sources( ) add_subdirectory(evolution) +add_subdirectory(functional) add_subdirectory(graph) add_subdirectory(graph_encoding) add_subdirectory(monomial_propagator) diff --git a/cpp/monoprop/detail/functional/CMakeLists.txt b/cpp/monoprop/detail/functional/CMakeLists.txt new file mode 100644 index 00000000..a8b518aa --- /dev/null +++ b/cpp/monoprop/detail/functional/CMakeLists.txt @@ -0,0 +1,8 @@ +target_sources( + monoprop + PUBLIC + FILE_SET headers + TYPE HEADERS + FILES + "Control.h" +) diff --git a/cpp/monoprop/detail/functional/Control.h b/cpp/monoprop/detail/functional/Control.h new file mode 100644 index 00000000..5c288a2c --- /dev/null +++ b/cpp/monoprop/detail/functional/Control.h @@ -0,0 +1,100 @@ +// 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 +#include +#include +#include + +#include "monoprop/TypeAliases.h" + +namespace monoprop::detail { + +// Initial-operator weights published together. +struct OperatorWeights { + VecD op; // One coefficient per store row. + double core_term{0.0}; // Identity contribution. + size_t structure_revision{0}; // Publication revision. +}; + +// Thread-safe slot for the current weights. +class WeightsSlot { +public: + auto load() const -> std::shared_ptr { + const std::lock_guard lock(mutex_); + return weights_; + } + + auto store(std::shared_ptr weights) -> void { + const std::lock_guard lock(mutex_); + weights_ = std::move(weights); + } + +private: + mutable std::mutex mutex_; // Allows reads through const FunctionalControl. + std::shared_ptr weights_; +}; + +// State shared by a propagator and its functionals. +struct FunctionalControl { + // Changes when replayed structure changes. + std::atomic structure_revision{0}; + + // Cleared before destruction. + std::atomic propagator_alive{true}; + + // Last structural change. + std::atomic last_structural_change{nullptr}; + + // Current weights, null until the first plan. + WeightsSlot weights; + + // Record a structure change. + auto bump(const char *site) -> void { + last_structural_change.store(site); + structure_revision.fetch_add(1); + } +}; + +// Bumps the control block if a mutation throws after a possible state change. +class [[nodiscard]] BumpOnUnwind { +public: + BumpOnUnwind(FunctionalControl &control, const char *site, bool armed = true) + : control_(control), + site_(site), + armed_(armed), + uncaught_(std::uncaught_exceptions()) {} + + BumpOnUnwind(const BumpOnUnwind &) = delete; + BumpOnUnwind(BumpOnUnwind &&) = delete; + auto operator=(const BumpOnUnwind &) -> BumpOnUnwind & = delete; + auto operator=(BumpOnUnwind &&) -> BumpOnUnwind & = delete; + + ~BumpOnUnwind() { + if (armed_ && std::uncaught_exceptions() > uncaught_) { + control_.bump(site_); + } + } + +private: + FunctionalControl &control_; + const char *site_; + bool armed_; + int uncaught_; +}; + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 117590ec..aaf10779 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -110,8 +110,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope } const size_t n_partitions = resolve_partition_count_(partitions, comm); - // The R ranks x S partitions form one flat P = R*S SPMD world, so a mismatch across ranks would - // deadlock at the first hybrid collective. + // All ranks need the same partition count for hybrid collectives. if (comm.kind == mpi::Comm::Kind::Mpi && mpi::size(comm) > 1 && mpi::allreduce_sum(n_partitions, comm) != n_partitions * static_cast(mpi::size(comm))) { throw PartitionCountMismatch( @@ -164,15 +163,15 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope auto op = schrodinger_ ? generate_paired_op(sc / 2 + sc % 2, logical_num_modes_) : 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_. + // packed_inline_width_() depends on cutoff_fn_. regenerate_cutoff_fn_(); mp_op_.store = std::make_unique>(packed_inline_width_()); mp_op_.store->reserve(expected_local_terms); - // Store replaced: drop the stale lazy inverted index so it rebuilds against the new store. + // Rebuild the lazy index for the new store. mp_op_.inverted_index_.reset(); size_t i = 0; - // The initial monomials are distinct, so emplace (insert-if-absent) is an assigning insert here. + // Initial monomials are distinct, so emplace assigns here. for (size_t r = 0; r < op.size(); ++r) { const auto &mono = materialize_row(op, r); if (my_rank == find_rank(mono, num_ranks)) { @@ -188,7 +187,10 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope } template -MonomialPropagator::~MonomialPropagator() = default; +MonomialPropagator::~MonomialPropagator() { + // Mark borrowed state unavailable before members are destroyed. + functional_control_->propagator_alive.store(false); +} template MonomialPropagator::MonomialPropagator(const MonomialPropagator &other) @@ -202,7 +204,6 @@ MonomialPropagator::MonomialPropagator(const MonomialPropagator &other lower_atol_(other.lower_atol_), upper_atol_(other.upper_atol_), core_term_(other.core_term_), - initial_operator_epoch_(other.initial_operator_epoch_), logical_num_modes_(other.logical_num_modes_), cutoff_type_(other.cutoff_type_), basis_change_(other.basis_change_), @@ -351,9 +352,13 @@ auto MonomialPropagator::packed_inline_width_() const -> size_t { template auto MonomialPropagator::apply_initial_operator_(const OperatorDict &op_dict) -> std::pair, VecD> { - ++initial_operator_epoch_; + // A failed re-weight cannot be followed: the facade fan-out can commit some partitions before + // another rejects a term only it holds, and a functional cannot tell that apart from a dict + // rejected before anything committed. Guarding the whole body keeps today's verdict for both. + auto guard = bump_structure_on_unwind_("update_initial_operator()"); if (partition_group_) { - // The facade holds no local terms of its own, so the return is empty. + // The facade holds no local terms of its own, so the return is empty. Each partition publishes + // its own weights on its own master, which is what its children's functionals read. for_each_partition_([&](MonomialPropagator &s) { s.update_initial_operator(op_dict); }); return {}; } @@ -361,10 +366,14 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o const size_t my_rank = static_cast(mpi::rank(comm_)); OperatorDict new_op; + // A dict that omits the identity term means zero, which is what MPOperator::update_initial_operator + // does with every row this dict leaves out; accumulated locally so a rejected dict leaves the + // propagator's own expectation value where it was. + double core_term = 0.0; for (const auto &[ind, coeff] : op_dict) { const auto mono = indices_to_bitset_checked(ind, 2 * logical_num_modes_); if (ind.empty()) { // Core term, store in all - core_term_ = algebra_encode_coeff(basis_, coeff, mono); + core_term = algebra_encode_coeff(basis_, coeff, mono); continue; } if (my_rank == find_rank(mono, num_ranks)) { @@ -373,7 +382,40 @@ auto MonomialPropagator::apply_initial_operator_(const OperatorDict &o } } - return mp_op_.update_initial_operator(new_op, schrodinger_); + // No revision bump: a re-weight leaves the store, the inverted index and the graph where they are, + // so a live functional follows the new coefficients instead of going stale. Publication comes after + // the commit, so a functional never sees half of one. + auto applied = mp_op_.update_initial_operator(new_op, schrodinger_); + core_term_ = core_term; + // Only refresh a set somebody has already been given: a null one means no plan was ever built, + // and weights_for_plan_() publishes on demand, so there is nothing left stale by skipping this. + // Publishing copies the whole coefficient vector, which a re-weight loop would otherwise pay per + // call on a propagator with no functionals. + if (functional_control_->weights.load() != nullptr) { + publish_weights_(); + } + return applied; +} + +template +auto MonomialPropagator::publish_weights_() -> std::shared_ptr { + // get_operator() merges pending terms and must run on the propagator thread. + auto weights = std::make_shared( + detail::OperatorWeights{.op = mp_op_.get_operator(), + .core_term = core_term_, + .structure_revision = functional_control_->structure_revision.load()}); + functional_control_->weights.store(weights); + return weights; +} + +template +auto MonomialPropagator::weights_for_plan_() -> std::shared_ptr { + // Weights from the current revision are still current. + if (auto published = functional_control_->weights.load(); + published != nullptr && published->structure_revision == functional_control_->structure_revision.load()) { + return published; + } + return publish_weights_(); } template @@ -625,16 +667,21 @@ auto MonomialPropagator::build_graph(const std::vector &majorana std::optional parameters, std::optional only_rotate_len_k) -> void { validate_only_rotate_len_k_(only_rotate_len_k, 2 * logical_num_modes_); + // Hoisted above the fan-out so a call that appends nothing bumps on neither shape: the children + // re-run both identically, so deciding here moves nothing observable. + if (majoranas.empty()) { + return; + } + validate_coefficient_lengths(parameter_mapping, gen_coeffs); if (partition_group_) { + // A fan-out that throws part-way has already mutated the partitions it reached. + auto guard = bump_structure_on_unwind_("build_graph()"); for_each_partition_([&](MonomialPropagator &s) { s.build_graph(majoranas, parameter_mapping, gen_coeffs, gate_indices, parameters, only_rotate_len_k); }); + bump_structure_("build_graph()"); return; } - if (majoranas.empty()) { - return; - } - validate_coefficient_lengths(parameter_mapping, gen_coeffs); VecZ local_gates; if (gate_indices.has_value()) { @@ -650,15 +697,14 @@ 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 { + // The seed is computed before the guard below: it validates and contracts out of place, so a + // rejection here has changed nothing yet. + VecD seed; + if (parameters.has_value()) { // 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 // realistic coefficients. That graph covers the parameter prefix [0, m). - VecD seed; if (graph_layers() > 0) { const auto existing = graph_gate_arrays_(); const size_t m = expected_num_params(existing.first); @@ -679,6 +725,15 @@ auto MonomialPropagator::build_graph(const std::vector &majorana else { seed = current_picture_coeffs_(); } + } + + // The gate loop bounds-checks each generator as it reaches it, so a bad generator in a multi-gate + // call throws with the earlier layers already appended. + auto guard = bump_structure_on_unwind_("build_graph()"); + if (!parameters.has_value()) { + evolve_mode_build_graph_(majoranas, parameter_mapping, gen_coeffs, local_gates, only_rotate_len_k); + } + else { evolve_mode_graph_with_coeffs_(majoranas, parameter_mapping, gen_coeffs, @@ -687,6 +742,7 @@ auto MonomialPropagator::build_graph(const std::vector &majorana seed, only_rotate_len_k); } + bump_structure_("build_graph()"); } template @@ -696,12 +752,8 @@ auto MonomialPropagator::propagate(const std::vector &majoranas, const VecD ¶meters, std::optional only_rotate_len_k) -> void { validate_only_rotate_len_k_(only_rotate_len_k, 2 * logical_num_modes_); - if (partition_group_) { - for_each_partition_([&](MonomialPropagator &s) { - s.propagate(majoranas, parameter_mapping, gen_coeffs, parameters, only_rotate_len_k); - }); - return; - } + // Hoisted above the fan-out so a call that evolves nothing bumps on neither shape: the children + // re-run all three identically (graph_layers() is facade-transparent), so nothing observable moves. if (majoranas.empty()) { return; } @@ -714,7 +766,20 @@ auto MonomialPropagator::propagate(const std::vector &majoranas, "build_graph() to extend it.", graph_layers())); } + if (partition_group_) { + // A fan-out that throws part-way has already mutated the partitions it reached. + auto guard = bump_structure_on_unwind_("propagate()"); + for_each_partition_([&](MonomialPropagator &s) { + s.propagate(majoranas, parameter_mapping, gen_coeffs, parameters, only_rotate_len_k); + }); + bump_structure_("propagate()"); + return; + } + // The gate loop bounds-checks each generator as it reaches it, so a bad generator in a multi-gate + // call throws with the earlier gates already folded into the operator. + auto guard = bump_structure_on_unwind_("propagate()"); evolve_mode_contract_immediately_(majoranas, parameter_mapping, gen_coeffs, parameters, only_rotate_len_k); + bump_structure_("propagate()"); } template @@ -810,6 +875,7 @@ template auto MonomialPropagator::set_parameter_mapping(const VecZ ¶meter_mapping) -> void { if (partition_group_) { for_each_partition_([&](MonomialPropagator &s) { s.set_parameter_mapping(parameter_mapping); }); + bump_structure_("set_parameter_mapping()"); return; } const size_t count = graph_.layers(); @@ -849,6 +915,7 @@ auto MonomialPropagator::set_parameter_mapping(const VecZ ¶meter_m count, gates)); } + bump_structure_("set_parameter_mapping()"); } template @@ -916,40 +983,48 @@ auto build_cos_callbacks(const detail::InvertedIndex &inverted_index, } template -template -auto MonomialPropagator::make_functional_(Fn &&func, std::optional pare_threshold) - -> std::function { +auto MonomialPropagator::make_plan_(std::optional pare_threshold) + -> std::shared_ptr> { + using Plan = detail::FunctionalPlan; + + if (partition_group_) { + // Both functional types share this fan-out plan. + typename Plan::Fanout fanout; + fanout.group = partition_group_.get(); + fanout.partitions = map_partitions_([&](MonomialPropagator &s) { return s.make_plan_(pare_threshold); }); + // Child plans share graph structure and parameter-axis length. + const auto num_params = fanout.partitions.front()->num_params(); + return std::make_shared(num_params, functional_control_, std::move(fanout)); + } + + typename Plan::Local local; + auto gate_arrays = graph_gate_arrays_(); - auto parameter_mapping = std::move(gate_arrays.first); - auto gen_coeffs = std::move(gate_arrays.second); - const auto num_params = expected_num_params(parameter_mapping); - - // Nothing here needs a dense state: energy only dots it against the evolved operator, and the gradient - // scatters it into its own thread-local scratch before back-evolving. So Heisenberg hands over just the - // sparse scores; Schrödinger's state is the live evolved vector, snapshotted whole. - // EvalState owns its rows and snapshots the term count -- a later append push_backs onto the operator's - // sparse rows, which would both dangle a view and outrun the `op` captured below. + local.parameter_mapping = std::move(gate_arrays.first); + local.gen_coeffs = std::move(gate_arrays.second); + const auto num_params = expected_num_params(local.parameter_mapping); + + // Heisenberg snapshots sparse scores; Schrodinger snapshots the dense state. const auto num_terms = mp_op_.size(); - auto state = [&] { + local.state = [&] { if (schrodinger_) { return EvalState::dense(mp_op_.dense_state()); } const auto sparse = mp_op_.sparse_state(); return EvalState::sparse(num_terms, sparse.rows, sparse.values); }(); - VecD op = mp_op_.get_operator(); - const auto core_term = this->core_term(); - const auto comm = comm_; - - const auto expected_layers = graph_layers(); - // Aliased rather than copied: the check below needs the live counter, like graph->layers(). - const auto *epoch = &initial_operator_epoch_; - const auto expected_epoch = initial_operator_epoch_; + local.weights = weights_for_plan_(); + local.comm = comm_; + const auto &inverted_index = mp_op_.inverted_index(); + local.mp_op = &mp_op_; + local.op_store = mp_op_.store.get(); + local.inverted_index_rows = inverted_index.rows(); - // One owning handle either way: pare hands back a heap-owned MPGraph the functional must keep alive - // (build_cos_callbacks holds pointers into its layers' stored cos); non-pare aliases graph_. - std::shared_ptr graph; + // Only coefficient-pared Schrodinger plans reject re-weighting. + local.pared_from_operator = schrodinger_ && pare_threshold.has_value(); + + // Own graph layers because `cos` stores raw pointers into them. if (pare_threshold.has_value()) { auto full_cos_of_layer = [this, &inverted_index](size_t i) -> CosMask { const auto layer = graph_.get_layer_traversal(i); @@ -957,84 +1032,44 @@ auto MonomialPropagator::make_functional_(Fn &&func, std::optional(inverted_index, gen, layer.scaled_count(), basis_); return detail::fold_to_cos_mask(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 = + // Threshold the Hamiltonian in Schrodinger, otherwise the state. + const auto keep = schrodinger_ ? indices_above(local.weights->op, *pare_threshold) + : local.state.indices_above(*pare_threshold); + const auto count = schrodinger_ ? local.weights->op.size() : local.state.length(); + local.graph = std::make_shared(pare_graph(graph_, keep, count, schrodinger_, comm_, full_cos_of_layer)); } else { - graph = std::shared_ptr(std::shared_ptr{}, &graph_); + // Snapshot active layers only; graph_ retains layers retired by slicing. + std::vector owned; + owned.reserve(graph_.layers()); + for (size_t i = 0; i < graph_.layers(); ++i) { + owned.push_back(graph_.get_layer(i)); + } + local.graph = std::make_shared(graph_.is_schrodinger(), std::move(owned)); } - // The folds keep raw column pointers into this propagator's inverted index, so the returned callable - // must not outlive the propagator. - auto cos = build_cos_callbacks(inverted_index, graph->replay_view(), basis_); - - return [func = std::move(func), - core_term, - state = std::move(state), - op = std::move(op), - graph = std::move(graph), - parameter_mapping, - gen_coeffs, - num_params, - epoch, - expected_epoch, - expected_layers, - cos = std::move(cos), - comm](const VecD ¶ms) -> R { - validate_expected_initial_operator(*epoch, expected_epoch); - validate_functional_call(params, num_params); - validate_expected_graph_layers(graph->layers(), expected_layers); - return func(EvalRequest{.e_core = core_term, - .state = state, - .op = op, - .parameter_mapping = parameter_mapping, - .gen_coeffs = gen_coeffs, - .graph = graph->replay_view(), - .params = params}, - comm, - cos); - }; + local.cos = build_cos_callbacks(inverted_index, local.graph->replay_view(), basis_); + + return std::make_shared(num_params, functional_control_, std::move(local)); } template auto MonomialPropagator::expectation_value_functional(std::optional pare_threshold) - -> std::function { - if (partition_group_) { - // Each partition allreduces internally, so partition 0 is the global value. The group is captured by - // raw pointer, so the returned callable must not outlive this propagator. - auto fns = std::make_shared>>( - map_partitions_([&](MonomialPropagator &s) { return s.expectation_value_functional(pare_threshold); })); - auto *grp = partition_group_.get(); - return [grp, fns](const VecD ¶ms) -> double { - return detail::partition::collect_on_all(*grp, - [&](int r) { return (*fns)[static_cast(r)](params); })[0]; - }; - } - return make_functional_(ev_fn, pare_threshold); + -> ExpectationValueFunctional { + return ExpectationValueFunctional(make_plan_(pare_threshold)); } template auto MonomialPropagator::expectation_value_and_gradient_functional(std::optional pare_threshold) - -> std::function(const VecD &)> { - if (partition_group_) { - auto fns = std::make_shared(const VecD &)>>>(map_partitions_( - [&](MonomialPropagator &s) { return s.expectation_value_and_gradient_functional(pare_threshold); })); - auto *grp = partition_group_.get(); - return [grp, fns](const VecD ¶ms) -> std::pair { - return detail::partition::collect_on_all(*grp, - [&](int r) { return (*fns)[static_cast(r)](params); })[0]; - }; - } - return make_functional_(ev_and_grad_fn, pare_threshold); + -> ExpectationValueAndGradientFunctional { + return ExpectationValueAndGradientFunctional(make_plan_(pare_threshold)); } template auto MonomialPropagator::expectation_value(const VecD ¶meters) -> double { if (partition_group_) { - // Each partition allreduces internally, so every partition returns the global value; take partition 0. + // Each partition returns the same allreduced value. return map_partitions_([&](MonomialPropagator &s) { return s.expectation_value(parameters); })[0]; } return expectation_value_functional(std::nullopt)(parameters); @@ -1043,7 +1078,7 @@ auto MonomialPropagator::expectation_value(const VecD ¶meters) -> template auto MonomialPropagator::expectation_value_and_gradient(const VecD ¶meters) -> std::pair { if (partition_group_) { - // As in expectation_value(): the gradient is allreduced inside each partition. + // Each partition returns the same allreduced gradient. return map_partitions_([&](MonomialPropagator &s) { return s.expectation_value_and_gradient(parameters); })[0]; } return expectation_value_and_gradient_functional(std::nullopt)(parameters); @@ -1051,8 +1086,21 @@ auto MonomialPropagator::expectation_value_and_gradient(const VecD &pa template auto MonomialPropagator::contract_partially(const VecD ¶meters, bool inplace) -> VecD { + // The one site name all three bump paths below report; bump_structure_ stores the pointer, so it has + // to be a literal that outlives every propagator. + static constexpr const char *kInplaceSite = "contract_partially(inplace=true)"; if (partition_group_) { - return concat_partitions_([&](MonomialPropagator &s) { return s.contract_partially(parameters, inplace); }); + // Decided before the fan-out, so an empty parameter vector -- which every child early-returns on + // without folding anything -- leaves a functional valid here exactly as it does on one partition. + validate_parameters_length(parameters, parameter_mapping()); + const bool folds = !parameters.empty(); + auto guard = bump_structure_on_unwind_(kInplaceSite, inplace && folds); + auto merged = + concat_partitions_([&](MonomialPropagator &s) { return s.contract_partially(parameters, inplace); }); + if (inplace && folds) { + bump_structure_(kInplaceSite); + } + return merged; } const auto gate_arrays = graph_gate_arrays_(); const auto ¶meter_mapping = gate_arrays.first; @@ -1060,10 +1108,14 @@ auto MonomialPropagator::contract_partially(const VecD ¶meters, bo validate_parameters_length(parameters, parameter_mapping); if (parameters.empty()) { + // Nothing is folded and the graph is left whole, so this is not a mutation even when inplace. return current_picture_coeffs_(); } const size_t num_majoranas = parameter_mapping.size(); + // Retiring the folded layers commits before the evolve that can throw, so the failure path has to + // bump as well; out of place nothing is written, so the guard stays disarmed there. + auto guard = bump_structure_on_unwind_(kInplaceSite, inplace); // 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_) { @@ -1074,6 +1126,9 @@ auto MonomialPropagator::contract_partially(const VecD ¶meters, bo 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; + // The graph lost the layers it folded and the state moved, so every plan built against + // either is stale. Out of place nothing changes, so there is no bump on that branch. + bump_structure_(kInplaceSite); } else { evolved_state = @@ -1089,6 +1144,7 @@ auto MonomialPropagator::contract_partially(const VecD ¶meters, bo 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; + bump_structure_(kInplaceSite); } else { evolved_op = evolve_operator_with_recompute_(VecD(op), graph_.slice_view(num_majoranas), mapped_params); diff --git a/cpp/tests/functional_validity.cpp b/cpp/tests/functional_validity.cpp new file mode 100644 index 00000000..69766744 --- /dev/null +++ b/cpp/tests/functional_validity.cpp @@ -0,0 +1,705 @@ +// 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. + +// What each public mutator does to a functional that was built before it ran. The table is the +// contract: one row per mutating method, asserted for the value and the gradient functional, with +// and without a pare threshold, in both pictures. + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/MonomialPropagator.h" +#include "monoprop/detail/functional/Control.h" +#include "monoprop/detail/mpi/MPICompat.h" + +using namespace monoprop; +namespace tt = boost::test_tools; + +namespace { + +constexpr size_t kNumModes = 2; +using Prop = MonomialPropagator; + +// Small enough to reason about by hand, large enough that both gates rotate something. +constexpr double kPareThreshold = 1e-12; +const VecD kBaseParams{0.3, 0.7}; +// One gate per Hamiltonian term, and the terms carry different weights, so the answer is not +// symmetric under swapping the two angles -- which is what makes the set_parameter_mapping row bite. +const std::vector kBaseGates{VecZ{0}, VecZ{2}}; + +// The weight mutate_update_initial_operator() writes onto term (0, 1). A re-weighted propagator must +// answer exactly like one built with it from the start, so both sides read it from here. +constexpr double kReweightedFirstWeight = 2.75; + +// The identity weight the core-term cases start from; every other case leaves the row out entirely, +// which is what makes a re-weight that also leaves it out a no-op there. +constexpr double kCoreTerm = 0.25; + +// `partitions` is passed explicitly, so it wins over the suite-wide monoprop_PARTITIONS=off. +// `core_term` adds the identity row, which only the core-term cases below need. +auto make_propagator(bool schrodinger, + size_t partitions = 1, + double first_weight = 1.0, + std::optional core_term = std::nullopt) -> Prop { + OperatorDict initial_ham; + initial_ham[VecZ{0, 1}] = std::complex{0.0, first_weight}; + initial_ham[VecZ{2, 3}] = std::complex{0.0, 0.5}; + if (core_term.has_value()) { + initial_ham[VecZ{}] = std::complex{*core_term, 0.0}; + } + const auto cutoff = static_cast(2 * kNumModes); + return Prop(initial_ham, + cutoff, + VecZ{0, 1}, + schrodinger ? std::optional{cutoff} : std::nullopt, + MPI_COMM_SELF, + std::nullopt, + std::nullopt, + CutoffType::Support, + std::nullopt, + kNumModes, + Basis::Majorana, + partitions); +} + +auto build_base_graph(Prop &prop) -> void { + prop.build_graph(kBaseGates, VecZ{0, 1}, VecD{1.0, 1.0}); +} + +// The mutators, one per public mutating method. Each runs against a propagator that already carries +// the two-layer base graph, except propagate(), which refuses a non-empty graph (see +// propagate_on_non_empty_graph_leaves_functional_valid). + +auto mutate_build_graph(Prop &prop) -> void { + prop.build_graph({VecZ{1}}, VecZ{0}, VecD{1.0}); +} + +auto mutate_propagate(Prop &prop) -> void { + prop.propagate({VecZ{0}}, VecZ{0}, VecD{1.0}, VecD{0.4}); +} + +auto mutate_contract_partially(Prop &prop) -> void { + prop.contract_partially(kBaseParams, /*inplace=*/true); +} + +auto mutate_update_initial_operator(Prop &prop) -> void { + OperatorDict updated; + updated[VecZ{0, 1}] = std::complex{0.0, kReweightedFirstWeight}; + updated[VecZ{2, 3}] = std::complex{0.0, 0.5}; + prop.update_initial_operator(updated); +} + +// Swaps which parameter drives which layer, so the answer moves whenever the two angles differ. +auto mutate_set_parameter_mapping(Prop &prop) -> void { + prop.set_parameter_mapping(VecZ{1, 0}); +} + +auto mutate_update_cutoff(Prop &prop) -> void { + prop.update_cutoff(2); +} + +auto mutate_update_cutoff_type(Prop &prop) -> void { + prop.update_cutoff_type(CutoffType::Length); +} + +auto mutate_update_basis_change(Prop &prop) -> void { + prop.update_basis_change(std::vector{VecZ{0}, VecZ{1}, VecZ{0, 1, 2}, VecZ{0, 1, 3}}); +} + +auto mutate_update_lower_atol(Prop &prop) -> void { + prop.update_lower_atol(1e-12); +} + +auto mutate_update_upper_atol(Prop &prop) -> void { + prop.update_upper_atol(1e-3); +} + +// What calling the pre-built functional does after the row's mutator ran. +enum class Outcome : std::uint8_t { + Stale, // throws, reporting the propagator moved under the functional + Answers, // returns, and returns exactly what it returned before the mutation + Refreshes, // returns, and now returns what a functional built after the mutation returns + RefusesRefresh, // throws, reporting weights it cannot follow rather than a moved structure +}; + +struct MutatorRow { + std::string_view method; // the public method this row covers + void (*apply)(Prop &); + bool needs_empty_graph; // build the functional with no graph, so the mutator is accepted + // Paring only changes the verdict where the keep-set came from the operator coefficients, so the + // Schrodinger-pared column is the only one that can differ from the general one. + Outcome outcome; // no pare_threshold in either picture, or pared Heisenberg + Outcome pared_schrodinger; // pare_threshold == kPareThreshold, Schrodinger: pares from `op` + std::string_view rationale; +}; + +constexpr std::array kMutatorTable{ + MutatorRow{.method = "build_graph", + .apply = &mutate_build_graph, + .needs_empty_graph = false, + .outcome = Outcome::Stale, + .pared_schrodinger = Outcome::Stale, + .rationale = "Appending a layer moves the structure revision, which a pared plan reads as " + "readily as an exact one."}, + MutatorRow{.method = "propagate", + .apply = &mutate_propagate, + .needs_empty_graph = true, + .outcome = Outcome::Stale, + .pared_schrodinger = Outcome::Stale, + .rationale = "Re-evolves the operator in place. It leaves the layer count at zero, so the " + "revision is the only thing that sees it."}, + MutatorRow{.method = "contract_partially", + .apply = &mutate_contract_partially, + .needs_empty_graph = false, + .outcome = Outcome::Stale, + .pared_schrodinger = Outcome::Stale, + .rationale = "Consumes the folded layers and rewrites the coefficients. Only inplace=true " + "bumps; see contract_partially_out_of_place_keeps_functional_valid."}, + MutatorRow{.method = "update_initial_operator", + .apply = &mutate_update_initial_operator, + .needs_empty_graph = false, + .outcome = Outcome::Refreshes, + .pared_schrodinger = Outcome::RefusesRefresh, + .rationale = "A re-weight moves no structure, so the functional follows the new " + "coefficients -- unless its keep-set was thresholded from those very " + "coefficients, which is Schrodinger with a pare threshold."}, + MutatorRow{.method = "set_parameter_mapping", + .apply = &mutate_set_parameter_mapping, + .needs_empty_graph = false, + .outcome = Outcome::Stale, + .pared_schrodinger = Outcome::Stale, + .rationale = "Relabels the layers in place, which changes neither the layer count nor the " + "operator -- the revision is the only thing that sees it."}, + MutatorRow{.method = "update_cutoff", + .apply = &mutate_update_cutoff, + .needs_empty_graph = false, + .outcome = Outcome::Answers, + .pared_schrodinger = Outcome::Answers, + .rationale = "Intended: a cutoff gates the next build and changes nothing the plan holds."}, + MutatorRow{.method = "update_cutoff_type", + .apply = &mutate_update_cutoff_type, + .needs_empty_graph = false, + .outcome = Outcome::Answers, + .pared_schrodinger = Outcome::Answers, + .rationale = "Intended: as update_cutoff."}, + MutatorRow{.method = "update_basis_change", + .apply = &mutate_update_basis_change, + .needs_empty_graph = false, + .outcome = Outcome::Answers, + .pared_schrodinger = Outcome::Answers, + .rationale = "Intended: as update_cutoff."}, + MutatorRow{.method = "update_lower_atol", + .apply = &mutate_update_lower_atol, + .needs_empty_graph = false, + .outcome = Outcome::Answers, + .pared_schrodinger = Outcome::Answers, + .rationale = "Intended: as update_cutoff."}, + MutatorRow{.method = "update_upper_atol", + .apply = &mutate_update_upper_atol, + .needs_empty_graph = false, + .outcome = Outcome::Answers, + .pared_schrodinger = Outcome::Answers, + .rationale = "Intended: as update_cutoff."}, +}; + +// A mutating method with no row would leave its effect on a live functional unrecorded, so the count +// is pinned to the roster on the class. Adding a mutator means bumping that constant, which breaks +// this build until a row lands here. +constexpr auto distinct_methods() -> size_t { + size_t distinct = 0; + for (size_t i = 0; i < kMutatorTable.size(); ++i) { + bool seen_earlier = false; + for (size_t j = 0; j < i; ++j) { + seen_earlier = seen_earlier || kMutatorTable[j].method == kMutatorTable[i].method; + } + distinct += seen_earlier ? 0U : 1U; + } + return distinct; +} + +static_assert(kMutatorTable.size() == Prop::num_mutating_methods && distinct_methods() == kMutatorTable.size(), + "cpp/tests/functional_validity.cpp must carry one row per public mutating method of " + "MonomialPropagator; see num_mutating_methods in MonomialPropagator.h."); + +// Both functional kinds behind one signature; the gradient kind is judged on its value component. +auto make_call(Prop &prop, bool gradient, std::optional pare_threshold) -> std::function { + if (gradient) { + auto fn = prop.expectation_value_and_gradient_functional(pare_threshold); + return [fn = std::move(fn)](const VecD ¶ms) { return fn(params).first; }; + } + auto fn = prop.expectation_value_functional(pare_threshold); + return [fn = std::move(fn)](const VecD ¶ms) { return fn(params); }; +} + +auto expected_outcome(const MutatorRow &row, bool schrodinger, std::optional pare_threshold) -> Outcome { + return pare_threshold.has_value() && schrodinger ? row.pared_schrodinger : row.outcome; +} + +auto run_row(const MutatorRow &row, bool schrodinger, bool gradient, std::optional pare_threshold) -> void { + auto prop = make_propagator(schrodinger); + if (!row.needs_empty_graph) { + build_base_graph(prop); + } + const VecD params = row.needs_empty_graph ? VecD{} : kBaseParams; + + auto call = make_call(prop, gradient, pare_threshold); + const double before = call(params); + + row.apply(prop); + + switch (expected_outcome(row, schrodinger, pare_threshold)) { + case Outcome::Stale: + BOOST_CHECK_EXCEPTION(call(params), std::runtime_error, [](const std::runtime_error &e) { + BOOST_TEST_INFO("message: " << e.what()); + return std::string_view(e.what()).find("MP object has been modified") != std::string_view::npos; + }); + return; + case Outcome::RefusesRefresh: + BOOST_CHECK_EXCEPTION(call(params), std::runtime_error, [](const std::runtime_error &e) { + BOOST_TEST_INFO("message: " << e.what()); + return std::string_view(e.what()).find("cannot follow the new weights") != std::string_view::npos; + }); + return; + case Outcome::Answers: { + // The plan replays its own snapshot, so the number cannot have moved. + double after = 0.0; + BOOST_REQUIRE_NO_THROW(after = call(params)); + BOOST_TEST(after == before, tt::tolerance(1e-12)); + return; + } + case Outcome::Refreshes: { + // The plan reads the propagator's live weights, so it must now agree exactly with a functional + // built after the mutation -- and disagree with what it answered before it. + double after = 0.0; + BOOST_REQUIRE_NO_THROW(after = call(params)); + BOOST_TEST(after == make_call(prop, gradient, pare_threshold)(params)); + BOOST_TEST(after != before); + return; + } + } +} + +auto run_table(bool schrodinger, std::optional pare_threshold) -> void { + for (const auto &row : kMutatorTable) { + for (const bool gradient : {false, true}) { + BOOST_TEST_CONTEXT("method=" << row.method << " gradient=" << gradient << " rationale=" << row.rationale) { + run_row(row, schrodinger, gradient, pare_threshold); + } + } + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(functional_validity_table_heisenberg_exact) { + run_table(/*schrodinger=*/false, std::nullopt); +} + +BOOST_AUTO_TEST_CASE(functional_validity_table_heisenberg_pared) { + run_table(/*schrodinger=*/false, kPareThreshold); +} + +BOOST_AUTO_TEST_CASE(functional_validity_table_schrodinger_exact) { + run_table(/*schrodinger=*/true, std::nullopt); +} + +BOOST_AUTO_TEST_CASE(functional_validity_table_schrodinger_pared) { + run_table(/*schrodinger=*/true, kPareThreshold); +} + +// contract_partially(inplace=false) mutates nothing, so it belongs outside the table: it must never +// invalidate a functional. +BOOST_AUTO_TEST_CASE(contract_partially_out_of_place_keeps_functional_valid) { + auto prop = make_propagator(/*schrodinger=*/false); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + const double before = call(kBaseParams); + + prop.contract_partially(kBaseParams, /*inplace=*/false); + + BOOST_TEST(call(kBaseParams) == before, tt::tolerance(1e-12)); +} + +// propagate() refuses to run on top of a stored graph, so it cannot invalidate a functional built +// against one: the rejection leaves the propagator untouched. +BOOST_AUTO_TEST_CASE(propagate_on_non_empty_graph_leaves_functional_valid) { + auto prop = make_propagator(/*schrodinger=*/false); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + const double before = call(kBaseParams); + + BOOST_CHECK_THROW(mutate_propagate(prop), std::runtime_error); + + BOOST_TEST(call(kBaseParams) == before, tt::tolerance(1e-12)); +} + +// Two of the table's Stale rows, spelled out: the propagator's own answer moves, and the functional +// refuses the call rather than following it half way. + +BOOST_AUTO_TEST_CASE(set_parameter_mapping_invalidates_functional_it_desynchronises) { + auto prop = make_propagator(/*schrodinger=*/false); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + const double before = call(kBaseParams); + + mutate_set_parameter_mapping(prop); + + BOOST_TEST(prop.expectation_value(kBaseParams) != before); + BOOST_CHECK_THROW(call(kBaseParams), std::runtime_error); +} + +BOOST_AUTO_TEST_CASE(pared_functional_is_invalidated_by_build_graph) { + auto prop = make_propagator(/*schrodinger=*/false); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, kPareThreshold); + const double before = call(kBaseParams); + + mutate_build_graph(prop); + + BOOST_TEST(prop.expectation_value(kBaseParams) != before); + BOOST_CHECK_THROW(call(kBaseParams), std::runtime_error); +} + +// A functional reads its propagator's operator index directly, so it cannot outlive it. The control +// block records the destruction, which is the only thing left to read once the handles are dangling. +BOOST_AUTO_TEST_CASE(functional_reports_a_destroyed_propagator) { + auto prop = std::make_unique(make_propagator(/*schrodinger=*/false)); + build_base_graph(*prop); + auto call = make_call(*prop, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + + prop.reset(); + + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, [](const std::runtime_error &e) { + return std::string_view(e.what()).find("has been destroyed") != std::string_view::npos; + }); +} + +// The functional objects report the axis they were built against, so a caller can size its parameter +// vector without going back to the propagator. +BOOST_AUTO_TEST_CASE(functional_reports_its_parameter_axis) { + auto prop = make_propagator(/*schrodinger=*/false); + BOOST_TEST(prop.expectation_value_functional().num_params() == 0U); + + build_base_graph(prop); + BOOST_TEST(prop.expectation_value_functional().num_params() == kBaseParams.size()); + BOOST_TEST(prop.expectation_value_and_gradient_functional().num_params() == kBaseParams.size()); + BOOST_TEST(prop.expectation_value_functional(kPareThreshold).num_params() == kBaseParams.size()); +} + +// The contract read off the object, for callers that hold a functional and not the propagator: only a +// Schrodinger plan pared against the operator's own coefficients refuses to follow a re-weight. +BOOST_AUTO_TEST_CASE(functional_reports_whether_it_follows_the_weights) { + for (const bool schrodinger : {false, true}) { + auto prop = make_propagator(schrodinger); + build_base_graph(prop); + BOOST_TEST(prop.expectation_value_functional().follows_weights()); + BOOST_TEST(prop.expectation_value_and_gradient_functional().follows_weights()); + BOOST_TEST(prop.expectation_value_functional(kPareThreshold).follows_weights() == !schrodinger); + } + // A facade answers for its partitions, which were all built the same way. + auto facade = make_propagator(/*schrodinger=*/true, /*partitions=*/2); + build_base_graph(facade); + BOOST_TEST(facade.expectation_value_functional().follows_weights()); + BOOST_TEST(!facade.expectation_value_functional(kPareThreshold).follows_weights()); +} + +// A facade's plan holds the partition group by raw pointer, so it needs the facade's own control block: +// the children cannot report a destruction that takes their group with it. +BOOST_AUTO_TEST_CASE(fanned_out_functional_reports_a_destroyed_facade) { + auto facade = std::make_unique(make_propagator(/*schrodinger=*/false, /*partitions=*/2)); + build_base_graph(*facade); + auto call = make_call(*facade, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + + facade.reset(); + + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, [](const std::runtime_error &e) { + return std::string_view(e.what()).find("has been destroyed") != std::string_view::npos; + }); +} + +// The facade bumps its own revision as it fans a mutation out, so the error is reported on the calling +// thread rather than surfacing out of a partition's collective. +BOOST_AUTO_TEST_CASE(fanned_out_functional_is_invalidated_by_build_graph) { + auto facade = make_propagator(/*schrodinger=*/false, /*partitions=*/2); + build_base_graph(facade); + auto call = make_call(facade, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + + mutate_build_graph(facade); + + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, [](const std::runtime_error &e) { + return std::string_view(e.what()).find("build_graph()") != std::string_view::npos; + }); +} + +// The refresh, end to end and at full precision: a re-weighted propagator's functional must answer what +// a propagator built with those coefficients answers, to the last bit. The two run the same arithmetic +// over the same store order, so anything less than equality means the refresh reached a different +// vector -- there is no rounding to hide behind here. +namespace { + +auto check_refresh_matches_fresh_propagator(bool gradient, std::optional pare_threshold, size_t partitions) + -> void { + auto reweighted = make_propagator(/*schrodinger=*/false, partitions); + build_base_graph(reweighted); + auto call = make_call(reweighted, gradient, pare_threshold); + const double before = call(kBaseParams); + mutate_update_initial_operator(reweighted); + + auto fresh = make_propagator(/*schrodinger=*/false, partitions, kReweightedFirstWeight); + build_base_graph(fresh); + + BOOST_TEST(call(kBaseParams) == make_call(fresh, gradient, pare_threshold)(kBaseParams)); + BOOST_TEST(call(kBaseParams) != before); +} + +} // namespace + +BOOST_AUTO_TEST_CASE(reweighted_functional_matches_a_fresh_propagator) { + for (const bool gradient : {false, true}) { + for (const auto pare_threshold : {std::optional{}, std::optional{kPareThreshold}}) { + for (const size_t partitions : {1U, 2U}) { + BOOST_TEST_CONTEXT("gradient=" << gradient << " pared=" << pare_threshold.has_value() + << " partitions=" << partitions) { + check_refresh_matches_fresh_propagator(gradient, pare_threshold, partitions); + } + } + } + } +} + +// The gradient follows the weights too, component by component: a re-weight scales what the backward +// pass carries, so a value-only check would pass on a gradient that stayed behind. +BOOST_AUTO_TEST_CASE(reweighted_gradient_matches_a_fresh_propagator) { + auto reweighted = make_propagator(/*schrodinger=*/false); + build_base_graph(reweighted); + auto fn = reweighted.expectation_value_and_gradient_functional(std::nullopt); + const auto before = fn(kBaseParams); + mutate_update_initial_operator(reweighted); + + auto fresh = make_propagator(/*schrodinger=*/false, /*partitions=*/1, kReweightedFirstWeight); + build_base_graph(fresh); + const auto expected = fresh.expectation_value_and_gradient_functional(std::nullopt)(kBaseParams); + + const auto after = fn(kBaseParams); + BOOST_TEST(after.second == expected.second, tt::per_element()); + BOOST_CHECK(after.second != before.second); +} + +// Invariant 4: Schrodinger thresholds its keep-set from the operator coefficients, so new coefficients +// select a different keep-set and the plan cannot replay this one. It says so rather than answering for a +// paring nobody asked for. Heisenberg thresholds the state, which a re-weight leaves alone. +BOOST_AUTO_TEST_CASE(pared_schrodinger_functional_refuses_to_follow_a_reweight) { + auto prop = make_propagator(/*schrodinger=*/true); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, kPareThreshold); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + + mutate_update_initial_operator(prop); + + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, [](const std::runtime_error &e) { + const std::string_view what(e.what()); + BOOST_TEST_INFO("message: " << what); + return what.find("pares its graph") != std::string_view::npos + && what.find("cannot follow the new weights") != std::string_view::npos; + }); + // The unpared plan over the same propagator follows the re-weight, so the refusal is the paring's and + // not the picture's. + BOOST_CHECK_NO_THROW(make_call(prop, /*gradient=*/false, std::nullopt)(kBaseParams)); +} + +// A rejected re-weight is not a refresh: MPOperator::update_initial_operator throws for a term the +// operator does not hold, and core_term_ may already carry the new value by then. A functional must +// report that rather than answer from weights the propagator disagrees with. +BOOST_AUTO_TEST_CASE(a_failed_reweight_invalidates_the_functional) { + auto prop = make_propagator(/*schrodinger=*/false); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + + OperatorDict unknown_term; + unknown_term[VecZ{0, 2}] = std::complex{0.0, 1.0}; + BOOST_CHECK_THROW(prop.update_initial_operator(unknown_term), std::runtime_error); + + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, [](const std::runtime_error &e) { + return std::string_view(e.what()).find("update_initial_operator()") != std::string_view::npos; + }); +} + +// Building a second functional must not look like a re-weight to the first: both are built over the one +// published weight set, so the pared Schrodinger plan -- the only one that refuses to follow a new set -- +// keeps answering. +BOOST_AUTO_TEST_CASE(building_another_functional_is_not_a_reweight) { + auto prop = make_propagator(/*schrodinger=*/true); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, kPareThreshold); + const double before = call(kBaseParams); + + auto other = make_call(prop, /*gradient=*/true, std::nullopt); + BOOST_CHECK_NO_THROW(other(kBaseParams)); + + BOOST_TEST(call(kBaseParams) == before); +} + +// The scope guard behind the failure-path rule, over a bare control block: a scope that returns +// normally leaves the revision alone, and one that throws bumps it once, recording its site. +BOOST_AUTO_TEST_CASE(bump_on_unwind_bumps_only_when_the_scope_throws) { + detail::FunctionalControl control; + + { + auto guard = detail::BumpOnUnwind(control, "site()"); + } + BOOST_TEST(control.structure_revision.load() == 0U); + BOOST_TEST((control.last_structural_change.load() == nullptr)); + + BOOST_CHECK_THROW( + [&] { + auto guard = detail::BumpOnUnwind(control, "site()"); + throw std::runtime_error("part-way"); + }(), + std::runtime_error); + BOOST_TEST(control.structure_revision.load() == 1U); + BOOST_TEST(std::string_view(control.last_structural_change.load()) == std::string_view("site()")); + + // Disarmed is how a known no-op opts out, so it must stay silent on both paths. + BOOST_CHECK_THROW( + [&] { + auto guard = detail::BumpOnUnwind(control, "other()", /*armed=*/false); + throw std::runtime_error("part-way"); + }(), + std::runtime_error); + BOOST_TEST(control.structure_revision.load() == 1U); +} + +// A gate generator is bounds-checked only as the gate loop reaches it, so a bad one in a multi-gate +// call throws with the earlier gates already committed: valid, out of range, valid appends before it +// throws in both pictures (Heisenberg walks the gates in reverse). +namespace { + +const std::vector kPartWayGates{VecZ{0}, VecZ{2 * kNumModes + 1}, VecZ{2}}; +const VecZ kPartWayMapping{0, 1, 2}; +const VecD kPartWayCoeffs{1.0, 1.0, 1.0}; + +auto reports_site(std::string_view site) { + return [site](const std::runtime_error &e) { + const std::string_view what(e.what()); + BOOST_TEST_INFO("message: " << what); + return what.find(site) != std::string_view::npos; + }; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(a_part_way_build_graph_failure_invalidates_the_functional) { + for (const bool schrodinger : {false, true}) { + BOOST_TEST_CONTEXT("schrodinger=" << schrodinger) { + auto prop = make_propagator(schrodinger); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(kBaseParams)); + const size_t layers_before = prop.graph_layers(); + + BOOST_CHECK_THROW(prop.build_graph(kPartWayGates, kPartWayMapping, kPartWayCoeffs), std::runtime_error); + + // The graph grew, so the call did mutate on its way to the throw. + BOOST_TEST(prop.graph_layers() > layers_before); + BOOST_CHECK_EXCEPTION(call(kBaseParams), std::runtime_error, reports_site("build_graph()")); + } + } +} + +// propagate() folds into the operator instead of appending, so nothing counts the mutation: without +// the failure-path bump the functional would keep answering for coefficients that had moved. +BOOST_AUTO_TEST_CASE(a_part_way_propagate_failure_invalidates_the_functional) { + for (const bool schrodinger : {false, true}) { + BOOST_TEST_CONTEXT("schrodinger=" << schrodinger) { + auto prop = make_propagator(schrodinger); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + BOOST_CHECK_NO_THROW(call(VecD{})); + + BOOST_CHECK_THROW(prop.propagate(kPartWayGates, kPartWayMapping, kPartWayCoeffs, VecD{0.3, 0.7, 0.5}), + std::runtime_error); + + BOOST_CHECK_EXCEPTION(call(VecD{}), std::runtime_error, reports_site("propagate()")); + } + } +} + +// A call that appends or folds nothing is not a mutation. The facade decides that before it fans out, +// so both partition counts must agree -- this is the off/auto divergence the table exists to rule out. +BOOST_AUTO_TEST_CASE(no_op_mutators_keep_a_functional_valid) { + for (const size_t partitions : {1U, 2U}) { + BOOST_TEST_CONTEXT("partitions=" << partitions) { + auto prop = make_propagator(/*schrodinger=*/false, partitions); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + const double before = call(VecD{}); + + // No graph, so there is nothing to fold and no layer to retire; the return is the current + // coefficients either way, which is what it was before the call. + const VecD folded = prop.contract_partially(VecD{}, /*inplace=*/true); + BOOST_TEST(folded == prop.contract_partially(VecD{}, /*inplace=*/false), tt::per_element()); + prop.build_graph({}, VecZ{}, VecD{}); + prop.propagate({}, VecZ{}, VecD{}, VecD{}); + BOOST_TEST(prop.graph_layers() == 0U); + + BOOST_TEST(call(VecD{}) == before, tt::tolerance(1e-12)); + } + } +} + +// The core term describes the dict that committed, like every other row: a re-weight that leaves the +// identity out zeroes it, which is what MPOperator::update_initial_operator does with the rows the dict +// omits. The live functional and the propagator's own answer must both see that. +BOOST_AUTO_TEST_CASE(a_reweight_that_drops_the_core_term_zeroes_it) { + auto prop = make_propagator(/*schrodinger=*/false, /*partitions=*/1, /*first_weight=*/1.0, kCoreTerm); + build_base_graph(prop); + auto call = make_call(prop, /*gradient=*/false, std::nullopt); + BOOST_TEST(prop.core_term() == kCoreTerm); + + mutate_update_initial_operator(prop); // carries no identity row + + BOOST_TEST(prop.core_term() == 0.0); + auto fresh = make_propagator(/*schrodinger=*/false, /*partitions=*/1, kReweightedFirstWeight); + build_base_graph(fresh); + BOOST_TEST(call(kBaseParams) == make_call(fresh, /*gradient=*/false, std::nullopt)(kBaseParams)); + BOOST_TEST(call(kBaseParams) == prop.expectation_value(kBaseParams)); +} + +// The other half of the same rule: a dict the store rejects commits nothing, so the core term stays +// where it was and the propagator's own expectation value is untouched. +BOOST_AUTO_TEST_CASE(a_rejected_reweight_leaves_the_core_term_alone) { + auto prop = make_propagator(/*schrodinger=*/false, /*partitions=*/1, /*first_weight=*/1.0, kCoreTerm); + build_base_graph(prop); + const double before = prop.expectation_value(kBaseParams); + + OperatorDict rejected; + rejected[VecZ{}] = std::complex{2.0, 0.0}; + rejected[VecZ{0, 2}] = std::complex{0.0, 1.0}; // a term the operator does not hold + BOOST_CHECK_THROW(prop.update_initial_operator(rejected), std::runtime_error); + + BOOST_TEST(prop.core_term() == kCoreTerm); + BOOST_TEST(prop.expectation_value(kBaseParams) == before, tt::tolerance(1e-12)); +} diff --git a/cpp/tests/update_initial_operator.cpp b/cpp/tests/update_initial_operator.cpp index 76dc74ac..850ebe3d 100644 --- a/cpp/tests/update_initial_operator.cpp +++ b/cpp/tests/update_initial_operator.cpp @@ -53,12 +53,12 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_updates_core_expval) { auto updated_fn = simulator.expectation_value_functional(std::nullopt); BOOST_TEST(updated_fn(empty_params) == 2.75, tt::tolerance(1e-12)); - // The functional built before the re-weight snapshotted the old coefficients, so it must reject - // the call rather than answer for an operator the propagator no longer holds. - BOOST_CHECK_THROW(expval_fn(empty_params), std::runtime_error); + // A re-weight moves no structure, so the functional built before it follows the new coefficients: + // the same call now gives the new core term, not the old one it was built over. + BOOST_TEST(expval_fn(empty_params) == 2.75, tt::tolerance(1e-12)); } -BOOST_AUTO_TEST_CASE(update_initial_operator_invalidates_gradient_functional) { +BOOST_AUTO_TEST_CASE(update_initial_operator_refreshes_gradient_functional) { constexpr size_t n_modes = 2; OperatorDict initial_ham; initial_ham[VecZ{}] = std::complex{1.0, 0.0}; @@ -82,7 +82,8 @@ BOOST_AUTO_TEST_CASE(update_initial_operator_invalidates_gradient_functional) { updated[VecZ{}] = std::complex{2.75, 0.0}; simulator.update_initial_operator(updated); - BOOST_CHECK_THROW(grad_fn(empty_params), std::runtime_error); + // As for the value functional: the gradient functional follows the weights rather than going stale. + BOOST_TEST(grad_fn(empty_params).first == 2.75, tt::tolerance(1e-12)); BOOST_TEST(simulator.expectation_value_and_gradient_functional(std::nullopt)(empty_params).first == 2.75, tt::tolerance(1e-12)); } diff --git a/cpp/tests/validation_tests.cpp b/cpp/tests/validation_tests.cpp index 173c3a2d..4de8a96d 100644 --- a/cpp/tests/validation_tests.cpp +++ b/cpp/tests/validation_tests.cpp @@ -17,6 +17,7 @@ #include #include +#include #include "monoprop/TypeAliases.h" #include "monoprop/Validation.h" @@ -56,9 +57,48 @@ BOOST_AUTO_TEST_CASE(validation_functional_call) { BOOST_CHECK_THROW(validate_functional_call(VecD{0.1}, 2), std::runtime_error); } -BOOST_AUTO_TEST_CASE(validation_expected_graph_layers) { - BOOST_CHECK_NO_THROW(validate_expected_graph_layers(3, 3)); - BOOST_CHECK_THROW(validate_expected_graph_layers(4, 3), std::runtime_error); +BOOST_AUTO_TEST_CASE(validation_functional_state) { + const FunctionalState healthy{.propagator_alive = true, + .current_revision = 3, + .expected_revision = 3, + .operator_layout_unchanged = true, + .last_structural_change = nullptr}; + BOOST_CHECK_NO_THROW(validate_functional_state(healthy)); + + auto destroyed = healthy; + destroyed.propagator_alive = false; + BOOST_CHECK_THROW(validate_functional_state(destroyed), std::runtime_error); + + // The revision names the mutation that moved the structure, so the message can point at it. + auto mutated = healthy; + mutated.current_revision = 4; + mutated.last_structural_change = "build_graph()"; + BOOST_CHECK_EXCEPTION(validate_functional_state(mutated), std::runtime_error, [](const auto &e) { + return std::string_view(e.what()).find("build_graph()") != std::string_view::npos; + }); + + // The backstop: the operator moved without a revision bump. + auto rebuilt = healthy; + rebuilt.operator_layout_unchanged = false; + BOOST_CHECK_THROW(validate_functional_state(rebuilt), std::runtime_error); +} + +BOOST_AUTO_TEST_CASE(validation_weight_refresh) { + const WeightRefresh followable{.weights_revision = 3, .expected_revision = 3, .may_follow_weights = true}; + BOOST_CHECK_NO_THROW(validate_weight_refresh(followable)); + + // What a pared Schrodinger functional reports: its keep-set came from the coefficients the re-weight + // replaced, so replaying it would answer for a paring nobody asked for. + auto refuses = followable; + refuses.may_follow_weights = false; + BOOST_CHECK_EXCEPTION(validate_weight_refresh(refuses), std::runtime_error, [](const auto &e) { + return std::string_view(e.what()).find("cannot follow the new weights") != std::string_view::npos; + }); + + // The backstop: weights from another revision reached a functional whose own revision still matches. + auto other_revision = followable; + other_revision.weights_revision = 2; + BOOST_CHECK_THROW(validate_weight_refresh(other_revision), std::runtime_error); } BOOST_AUTO_TEST_CASE(validation_only_rotate_len_k) { diff --git a/cspell.json b/cspell.json index b1f6a54a..e622cfed 100644 --- a/cspell.json +++ b/cspell.json @@ -56,6 +56,8 @@ "qubit", "qubits", "Remigio", + "reweight", + "reweighted", "tracemalloc", "uct", "ucx", diff --git a/docs/content/docs/features/evaluation.mdx b/docs/content/docs/features/evaluation.mdx index 470854b4..38278373 100644 --- a/docs/content/docs/features/evaluation.mdx +++ b/docs/content/docs/features/evaluation.mdx @@ -47,15 +47,38 @@ expval, grad = expval_grad_fn(parameters) Shorter aliases are also available: `sim.expval_functional()` and `sim.expval_and_grad_functional()`. -A functional is built against the graph and initial-operator coefficients present when it is -created, so mutating either — -[build_graph][monoprop.monomial_propagator.MonomialPropagator.build_graph], -[contract_partially][monoprop.monomial_propagator.MonomialPropagator.contract_partially], -[update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator] — -invalidates it: calling it afterwards raises `RuntimeError` rather than returning a value for state -the propagator no longer holds. Build a new functional after such a call. The direct -[expectation_value][monoprop.monomial_propagator.MonomialPropagator.expectation_value] path always -reflects the current operator. +A functional is built against the graph present when it is created, and it reads through to the +propagator it came from, which it keeps alive for as long as you hold it. A structural change to that +graph invalidates the functional: calling it afterwards raises `RuntimeError` rather than returning a +value for a circuit the propagator no longer holds. Build a new functional after such a call. The +direct [expectation_value][monoprop.monomial_propagator.MonomialPropagator.expectation_value] path +always reflects the current propagator. + +A re-weight is the one mutation a functional survives: +[update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator] +leaves the graph, the operator's terms and their index where they are, so a live functional follows +the new coefficients. The price is that a functional is a *live view* of those coefficients — two +calls with the same parameters give two answers across a re-weight. Build the functional again after +the last re-weight if you need a frozen value. Each functional carries the rule as a +`follows_weights` attribute, so a caller holding one need not re-derive it from the picture and the +threshold. + +This is what every public mutating method does to a functional built before it ran: + +| Method | A call afterwards | +| --- | --- | +| [build_graph][monoprop.monomial_propagator.MonomialPropagator.build_graph] | raises | +| [propagate][monoprop.monomial_propagator.MonomialPropagator.propagate] | raises | +| [contract_partially][monoprop.monomial_propagator.MonomialPropagator.contract_partially] with `inplace=True` | raises | +| [parameter_mapping][monoprop.monomial_propagator.MonomialPropagator.parameter_mapping] | raises | +| [update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator] | answers for the new coefficients — except in the Schrödinger picture with a `pare_threshold`, where it raises, because the pared graph was selected from the coefficients the re-weight replaced | +| a rejected `update_initial_operator` | raises: a partitioned propagator applies the re-weight one partition at a time, so a term only one of them holds is rejected with its siblings already re-weighted | +| a `build_graph` or `propagate` that fails part-way | raises: a gate generator is bounds-checked only as the gate loop reaches it, so the gates before it are already committed | +| `cutoff`, `cutoff_type`, `basis_change`, `lower_atol`, `upper_atol` | answers, unchanged: these gate the next build and touch nothing a functional replays | +| `contract_partially` with `inplace=False` | answers, unchanged: it mutates nothing | +| a call that folds or appends nothing — `contract_partially([], inplace=True)`, `build_graph` or `propagate` with no gates | answers, unchanged: nothing moved, on any partition count | + +The table is executable — see [the mutation table](/testing#the-mutation-table). Both functionals accept an optional `pare_threshold` — see *Paring* below. diff --git a/docs/content/docs/features/initialisation.mdx b/docs/content/docs/features/initialisation.mdx index 78d49cd9..36920229 100644 --- a/docs/content/docs/features/initialisation.mdx +++ b/docs/content/docs/features/initialisation.mdx @@ -44,6 +44,10 @@ Schrödinger picture the state is evolved independently, so you can change the observable arbitrarily and read it off against the same evolved state; in the Heisenberg picture the observable can be re-weighted. +Functionals built before the update stay valid and answer for the new coefficients, so a re-weight +loop does not have to rebuild them — see [Reusable functionals](/features/evaluation#reusable-functionals) +for the one exception and for what a live view of the weights costs. + Only terms that already exist in the initial operator can be updated; supplying a monomial that is not present (or a non-Hermitian coefficient) raises diff --git a/docs/content/docs/testing.mdx b/docs/content/docs/testing.mdx index 51326cff..118a3faf 100644 --- a/docs/content/docs/testing.mdx +++ b/docs/content/docs/testing.mdx @@ -194,6 +194,18 @@ BOOST_AUTO_TEST_CASE(my_basic_check) { } ``` +#### The mutation table + +`cpp/tests/functional_validity.cpp` records, for every public method that mutates a +`MonomialPropagator`, what a functional built before that call does when called after it: throw +("stale"), or answer from its own snapshot ("answers"). Each row is asserted for the value and the +gradient functional, with and without a pare threshold, in both pictures. + +The roster is pinned by `MonomialPropagator::num_mutating_methods`, and the table `static_assert`s +its row count against it — so adding a mutating method means bumping that constant, which breaks the +build until the new method gets a row. `tests/test_parameter_validation.py::TestFunctionalValidityTable` +mirrors the same rows through the Python front end, over both `monoprop_PARTITIONS=off` and `=auto`. + #### Data-driven test using reference fixtures Use the `ExampleDataFix` fixture class from `TestUtilities.h` to load the same msgpack data as Python tests, and `BOOST_DATA_TEST_CASE_F` to parametrize over it: diff --git a/src/monoprop/bindings/binder.h b/src/monoprop/bindings/binder.h index 9aad5003..d9a18b8f 100644 --- a/src/monoprop/bindings/binder.h +++ b/src/monoprop/bindings/binder.h @@ -50,6 +50,21 @@ auto cutoff_type_enum_2_str(CutoffType cutoff_type) -> std::string; auto basis_str_2_enum(const std::string &basis) -> Basis; auto basis_enum_2_str(Basis basis) -> std::string; +// Opaque Python functional types, created only by their propagator factories. +template +auto bind_functional(nb::module_ &mod, const std::string &name, const char *call_doc) -> void { + nb::class_(mod, name.c_str()) + .def("__call__", &Functional::operator(), "parameters"_a, call_doc) + .def_prop_ro( + "num_params", + [](const Functional &f) { return f.num_params(); }, + "Parameter-axis length this functional was built against") + .def_prop_ro( + "follows_weights", + [](const Functional &f) { return f.follows_weights(); }, + "Whether a call after update_initial_operator() answers for the new weights"); +} + template auto bind_monomial_propagator(nb::module_ &mod) -> void { using namespace monoprop; @@ -57,6 +72,14 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { auto name = std::format("MonomialPropagator{:03d}", NumModes); auto cls = nb::class_>(mod, name.c_str()); + bind_functional>(mod, + std::format("ExpectationValueFunctional{:03d}", NumModes), + "Expectation value at the given variational parameters"); + bind_functional>( + mod, + std::format("ExpectationValueAndGradientFunctional{:03d}", NumModes), + "(expectation value, gradient) at the given variational parameters"); + cls.def( "__init__", [](MonomialPropagator *t, @@ -134,14 +157,17 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { "parameters"_a, "Expectation value and its gradient at the given variational parameters"); + // Functionals borrow this propagator's index and graph. cls.def("expectation_value_functional", &MonomialPropagator::expectation_value_functional, "pare_threshold"_a = std::nullopt, + nb::keep_alive<0, 1>(), "Reusable callable giving the expectation value from parameters; None keeps the exact graph"); cls.def("expectation_value_and_gradient_functional", &MonomialPropagator::expectation_value_and_gradient_functional, "pare_threshold"_a = std::nullopt, + nb::keep_alive<0, 1>(), "Reusable callable giving (expectation value, gradient) from parameters; None keeps the exact graph"); cls.def("contract_partially", @@ -153,7 +179,8 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { cls.def("update_initial_operator", &MonomialPropagator::update_initial_operator, "op_dict"_a, - "Rewrite the initial operator from an {indices: coefficient} dict"); + "Rewrite the initial operator from an {indices: coefficient} dict; every existing term " + "the dict omits is zeroed, the identity term included"); cls.def_prop_rw("lower_atol", &MonomialPropagator::lower_atol, diff --git a/src/monoprop/majorana_propagator.py b/src/monoprop/majorana_propagator.py index 477b70fd..0c4b491f 100644 --- a/src/monoprop/majorana_propagator.py +++ b/src/monoprop/majorana_propagator.py @@ -128,8 +128,9 @@ def update_initial_operator(self, new_operator: MajoranaOperator) -> None: Re-weights the initial operator the graph is evaluated against, without touching the evolution graph or rebuilding the simulator. Only the initial operator is affected -- the gates and their generator coefficients are unchanged. Functionals created - earlier are invalidated; see - [update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator]. + earlier follow the new coefficients rather than being invalidated; see + [update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator] + for the exception and the price. Args: new_operator: A [MajoranaOperator][monoprop.majorana.MajoranaOperator] whose terms diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py index 63be1c0b..73bb151f 100644 --- a/src/monoprop/monomial_propagator.py +++ b/src/monoprop/monomial_propagator.py @@ -26,7 +26,7 @@ import logging from abc import ABC, abstractmethod -from typing import TYPE_CHECKING, Generic, TypeVar +from typing import TYPE_CHECKING, Generic, Protocol, TypeVar import numpy as np @@ -43,7 +43,7 @@ from .utils import validate_basis_change if TYPE_CHECKING: - from collections.abc import Callable, Sequence + from collections.abc import Sequence from typing import Self from mpi4py import MPI @@ -55,6 +55,56 @@ T_op = TypeVar("T_op", MajoranaOperator, PauliOperator) +T_ret_co = TypeVar("T_ret_co", covariant=True) + + +class _EngineFunctional(Protocol[T_ret_co]): + """Interface exposed by engine functionals.""" + + def __call__(self, parameters: list[float], /) -> T_ret_co: ... + + @property + def num_params(self) -> int: ... + + @property + def follows_weights(self) -> bool: ... + + +class _BoundFunctionalBase(Generic[T_ret_co]): + """Shared parameter binding and forwarded functional attributes.""" + + def __init__( + self, propagator: MonomialPropagator, functional: _EngineFunctional[T_ret_co] + ) -> None: + self._propagator = propagator + self._functional = functional + + @property + def num_params(self) -> int: + """Parameter-axis length this functional was built against.""" + return self._functional.num_params + + @property + def follows_weights(self) -> bool: + """Whether a call after [update_initial_operator][] answers for the new coefficients.""" + return self._functional.follows_weights + + +class _BoundFunctional(_BoundFunctionalBase[float]): + """One engine functional returning the expectation value.""" + + def __call__(self, parameters: ParameterValues = None) -> float: + return self._functional(self._propagator._bind(parameters)) + + +class _BoundGradientFunctional(_BoundFunctionalBase[tuple[float, np.ndarray]]): + """As the value functional, but returning ``(value, gradient)`` with the gradient as ``float64``.""" + + def __call__(self, parameters: ParameterValues = None) -> tuple[float, np.ndarray]: + value, grad = self._functional(self._propagator._bind(parameters)) + return value, np.asarray(grad, dtype=np.float64) + + class MonomialPropagator(ABC, Generic[T_op]): """Abstract base for the classical monomial-propagation simulators. @@ -379,12 +429,18 @@ def gradient( def expectation_value_functional( self, pare_threshold: float | None = None - ) -> Callable[..., float]: + ) -> _BoundFunctional: """Return a reusable callable computing the expectation value from parameters. - The callable is built against the graph and initial-operator coefficients present now, so - mutating either -- [build_graph][], [contract_partially][], [update_initial_operator][] -- - invalidates it; build a new one after such a call. + The callable is built against the graph present now, so a structural change to it -- + [build_graph][], [propagate][], [contract_partially][] with ``inplace=True``, or a new + [parameter_mapping][] -- invalidates the callable; build a new one after such a call. + + A re-weight is not a structural change: the callable follows the propagator's current + initial-operator coefficients, so an [update_initial_operator][] leaves it valid and it + answers for the new weights. The price is that it is a live view of those weights rather + than a frozen number -- two calls with the same parameters give two answers across a + re-weight. Build the callable again after the last re-weight to freeze a value. Args: pare_threshold: Edge-retention cutoff for this functional's masked plan: edges @@ -392,40 +448,44 @@ def expectation_value_functional( accuracy for speed. ``None`` (default) disables paring. Returns: - A callable ``fn(parameters=None) -> float``. + A callable ``fn(parameters=None) -> float``, exposing the rule above as + ``follows_weights`` and its parameter-axis length as ``num_params``. Raises: - RuntimeError: From the returned callable, if the propagator was mutated after this - functional was created. + RuntimeError: From the returned callable, if the propagator was structurally mutated + after this functional was created, or if it was re-weighted and this functional + cannot follow the new weights -- which is the Schrodinger picture with a + ``pare_threshold``, whose pared graph was selected from the very coefficients the + re-weight replaced. """ - fn = self._simulator.expectation_value_functional(pare_threshold) - return lambda parameters=None: fn(self._bind(parameters)) + return _BoundFunctional( + self, self._simulator.expectation_value_functional(pare_threshold) + ) def expectation_value_and_gradient_functional( self, pare_threshold: float | None = None - ) -> Callable[..., tuple]: + ) -> _BoundGradientFunctional: """Return a reusable callable computing (expectation value, gradient). Like [expectation_value_functional][], but one backward pass also yields the gradient. It is - invalidated by the same mutations. + invalidated by the same mutations, and follows the initial-operator weights on the same + terms: value and gradient both answer for the current coefficients. Args: pare_threshold: See [expectation_value_functional][]. Returns: - A callable ``fn(parameters=None) -> (float, np.ndarray)``, gradient in parameter order. + A callable ``fn(parameters=None) -> (float, np.ndarray)``, gradient in parameter order, + with the same ``follows_weights`` and ``num_params`` as [expectation_value_functional][]. Raises: - RuntimeError: From the returned callable, if the propagator was mutated after this - functional was created. + RuntimeError: From the returned callable, on the same conditions as + [expectation_value_functional][]. """ - fn = self._simulator.expectation_value_and_gradient_functional(pare_threshold) - - def _call(parameters=None): # noqa: ANN001, ANN202 - value, grad = fn(self._bind(parameters)) - return value, np.asarray(grad, dtype=np.float64) - - return _call + return _BoundGradientFunctional( + self, + self._simulator.expectation_value_and_gradient_functional(pare_threshold), + ) def expval( self, @@ -459,7 +519,7 @@ def expval_and_grad( def expval_functional( self, pare_threshold: float | None = None - ) -> Callable[..., float]: + ) -> _BoundFunctional: """Shorthand for [expectation_value_functional][]. See [expectation_value_functional][] for full documentation. @@ -468,7 +528,7 @@ def expval_functional( def expval_and_grad_functional( self, pare_threshold: float | None = None - ) -> Callable[..., tuple]: + ) -> _BoundGradientFunctional: """Shorthand for [expectation_value_and_gradient_functional][]. See [expectation_value_and_gradient_functional][] for full documentation. @@ -541,13 +601,23 @@ def update_initial_operator(self, new_operator: T_op) -> None: Each concrete front-end implements this over its own operator type, encoding the terms into the engine's raw index tuples. - Functionals hold the coefficients they were built with, so any created earlier are - invalidated -- they raise instead of answering for the replaced operator. + Functionals created earlier stay valid and follow the new coefficients: a re-weight moves no + structure, so there is nothing for them to be stale about. They are therefore a live view of + the weights -- two calls with the same parameters give two answers across this call. The one + exception is a functional built in the Schrodinger picture with a ``pare_threshold``: its + pared graph was selected from the coefficients this call replaces, so it raises instead of + following them. + + A rejected re-weight invalidates them: a partitioned propagator applies the new weights one + partition at a time, so a term only one of them holds is rejected with its siblings already + re-weighted, and a functional cannot tell that apart from a rejection that committed nothing. + They raise rather than answer for weights the propagator disagrees with. Args: new_operator: A [MajoranaOperator][monoprop.majorana.MajoranaOperator] or [PauliOperator][monoprop.pauli.PauliOperator], per the front-end, whose terms replace - the matching initial-operator. + the matching initial-operator. Every existing term it leaves out is zeroed, the + identity term included. Raises: RuntimeError: In the Heisenberg picture, if a term is absent from the current operator. diff --git a/src/monoprop/pauli_propagator.py b/src/monoprop/pauli_propagator.py index 1b0ec4e8..de026c11 100644 --- a/src/monoprop/pauli_propagator.py +++ b/src/monoprop/pauli_propagator.py @@ -221,8 +221,9 @@ def update_initial_operator(self, new_operator: PauliOperator) -> None: base method, which takes the engine's raw symplectic-slot keys, this accepts qubit Pauli terms and encodes them via [get_local_operator][monoprop.pauli.PauliOperator.get_local_operator]. Functionals created - earlier are invalidated; see - [update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator]. + earlier follow the new coefficients rather than being invalidated; see + [update_initial_operator][monoprop.monomial_propagator.MonomialPropagator.update_initial_operator] + for the exception and the price. Args: new_operator: A [PauliOperator][monoprop.pauli.PauliOperator] whose terms replace the diff --git a/tests/test_parameter_validation.py b/tests/test_parameter_validation.py index a38bbbea..423cc1d8 100644 --- a/tests/test_parameter_validation.py +++ b/tests/test_parameter_validation.py @@ -16,13 +16,62 @@ from __future__ import annotations +import numpy as np import pytest -from monoprop import Circuit, ExpGate, MajoranaPropagator, PauliPropagator +from monoprop import ( + Circuit, + ExpGate, + MajoranaPropagator, + PauliPropagator, + jordan_wigner_basis_change, +) from monoprop.majorana import MajoranaOperator from monoprop.pauli import PauliOperator +def _components(result): + """A functional's answer as ``(value, gradient)``; ``None`` gradient for the value-only kind.""" + return result if isinstance(result, tuple) else (result, None) + + +def _assert_answers_match(actual, expected, *, exact, context=""): + """Assert both components of two answers agree, bit-exactly or to ``pytest.approx``. + + The two functional kinds return different shapes, so a test parametrized over both has to + compare whatever the kind under test returned -- a gradient dropped here is a gradient nobody + checks. + """ + value, gradient = _components(actual) + expected_value, expected_gradient = _components(expected) + assert (gradient is None) == (expected_gradient is None), context + if exact: + assert value == expected_value, context + if gradient is not None: + assert np.array_equal(gradient, expected_gradient), context + else: + assert value == pytest.approx(expected_value), context + if gradient is not None: + assert gradient == pytest.approx(expected_gradient), context + + +def _assert_answers_differ(actual, other, context=""): + """Assert both components of two answers moved (a gradient counts as moved if any entry did).""" + value, gradient = _components(actual) + other_value, other_gradient = _components(other) + assert value != pytest.approx(other_value), context + if gradient is not None: + assert gradient != pytest.approx(other_gradient), context + + +# The propagator method each functional factory is the reusable form of, so a test can compare a +# functional's answer against the direct call of the same shape. +_DIRECT_CALL = { + "expectation_value_functional": "expectation_value", + "expectation_value_and_gradient_functional": "expectation_value_and_gradient", +} + + def _two_gate_graph(serial_comm): """A propagator with a two-layer, two-parameter graph already built.""" operator = MajoranaOperator({(0, 1): 1.0j, (2, 3): 0.5j}, num_modes=2) @@ -163,7 +212,7 @@ def test_functional_invalidated_after_graph_mutation(self, serial_comm): "expectation_value_and_gradient_functional", ], ) - def test_functional_invalidated_after_initial_operator_update( + def test_functional_follows_initial_operator_update( self, serial_comm, propagator_cls, @@ -186,18 +235,529 @@ def test_functional_invalidated_after_initial_operator_update( ) functional = getattr(mp, functional_name)() parameters = [0.3, 0.7] - functional(parameters) + before = functional(parameters) mp.update_initial_operator(updated_operator) - # The functional snapshotted the old coefficients, so it must reject the call rather than - # keep answering for the operator the propagator no longer holds. - with pytest.raises(RuntimeError, match=r"MP object has been modified"): + # A re-weight moves no structure, so the functional built before it follows the new + # coefficients instead of refusing the call: it now answers what the propagator answers. + after = functional(parameters) + _assert_answers_differ(after, before) + direct = getattr(mp, _DIRECT_CALL[functional_name])(parameters) + _assert_answers_match(after, direct, exact=False) + _assert_answers_match( + after, getattr(mp, functional_name)()(parameters), exact=True + ) + + +class TestFunctionalValidityTable: + """What each public mutator does to a functional built before it ran. + + The Python mirror of ``cpp/tests/functional_validity.cpp``: same rows, same expectations, run + over both partition settings. The gate that a new mutator gets a row is the C++ static_assert + against ``MonomialPropagator::num_mutating_methods``; these rows mirror that table. + """ + + _MODES = 2 + _CUTOFF = 4 + _PARAMS = (0.3, 0.7) + _PARE_THRESHOLD = 1e-12 + # What _mutate_update_initial_operator() writes onto term (0, 1); a re-weighted propagator must + # answer exactly like one built with it from the start, so both sides read it from here. + _REWEIGHTED_FIRST_WEIGHT = 2.75 + + # One gate per Hamiltonian term, and the terms carry different weights, so the answer is not + # symmetric under swapping the two angles -- which is what makes the parameter_mapping row bite. + @staticmethod + def _generator(index): + return MajoranaOperator({(index,): 1.0}, num_modes=2) + + @classmethod + def _base_circuit(cls): + return Circuit( + (ExpGate(cls._generator(0)), ExpGate(cls._generator(2))), + cls._MODES, + cls._PARAMS, + ) + + @classmethod + def _propagator( + cls, comm, *, schrodinger, with_graph, first_weight=1.0, core_term=None + ): + # The identity row is only carried by the core-term case; every other case leaves it out, + # which is what makes a re-weight that also leaves it out a no-op there. + terms = {(0, 1): first_weight * 1j, (2, 3): 0.5j} + if core_term is not None: + terms[()] = core_term + mp = MajoranaPropagator( + MajoranaOperator(terms, num_modes=cls._MODES), + [0, 1], + cutoff=cls._CUTOFF, + schrodinger_cutoff=cls._CUTOFF if schrodinger else None, + comm=comm, + ) + if with_graph: + mp.build_graph(cls._base_circuit()) + return mp + + # The mutators, one per public mutating method. + @classmethod + def _mutate_build_graph(cls, mp): + mp.build_graph(Circuit((ExpGate(cls._generator(1)),), cls._MODES, (0.4,))) + + @classmethod + def _mutate_propagate(cls, mp): + mp.propagate(Circuit((ExpGate(cls._generator(0)),), cls._MODES, (0.4,))) + + @classmethod + def _mutate_contract_partially(cls, mp): + mp.contract_partially(list(cls._PARAMS), inplace=True) + + @classmethod + def _mutate_update_initial_operator(cls, mp): + mp.update_initial_operator( + MajoranaOperator( + {(0, 1): cls._REWEIGHTED_FIRST_WEIGHT * 1j, (2, 3): 0.5j}, + num_modes=cls._MODES, + ) + ) + + @staticmethod + def _mutate_parameter_mapping(mp): + mp.parameter_mapping = [1, 0] + + @staticmethod + def _mutate_cutoff(mp): + mp.cutoff = 2 + + @staticmethod + def _mutate_cutoff_type(mp): + mp.cutoff_type = "length" + + @classmethod + def _mutate_basis_change(cls, mp): + # No front-end setter; the engine property is the only way in (see tests/test_basis.py). + mp._simulator.basis_change = jordan_wigner_basis_change(cls._MODES) + + @staticmethod + def _mutate_lower_atol(mp): + mp.lower_atol = 1e-12 + + @staticmethod + def _mutate_upper_atol(mp): + mp.upper_atol = 1e-3 + + # (method, mutator, needs_empty_graph, outcome, pared_schrodinger, rationale). Paring only changes + # the verdict where the keep-set came from the operator coefficients, so the Schrodinger-pared + # column is the only one that can differ from the general one -- as in the C++ table. + # "stale" = the call must throw that the propagator moved, "answers" = it must return exactly what + # it returned before the mutation, "refreshes" = it must return what a functional built after the + # mutation returns, "refuses-refresh" = it must throw that it cannot follow the new weights. + ROWS = ( + ( + "build_graph", + "_mutate_build_graph", + False, + "stale", + "stale", + "Appending a layer moves the structure revision, which a pared plan reads as readily as " + "an exact one.", + ), + ( + "propagate", + "_mutate_propagate", + True, + "stale", + "stale", + "Re-evolves the operator in place. It leaves the layer count at zero, so the revision " + "is the only thing that sees it.", + ), + ( + "contract_partially", + "_mutate_contract_partially", + False, + "stale", + "stale", + "Consumes the folded layers and rewrites the coefficients. Only inplace=True bumps.", + ), + ( + "update_initial_operator", + "_mutate_update_initial_operator", + False, + "refreshes", + "refuses-refresh", + "A re-weight moves no structure, so the functional follows the new coefficients -- " + "unless its keep-set was thresholded from those very coefficients, which is " + "Schrodinger with a pare threshold.", + ), + ( + "parameter_mapping", + "_mutate_parameter_mapping", + False, + "stale", + "stale", + "Relabels the layers in place, which changes neither the layer count nor the operator -- " + "the revision is the only thing that sees it.", + ), + ( + "cutoff", + "_mutate_cutoff", + False, + "answers", + "answers", + "Intended: a cutoff gates the next build and changes nothing the plan holds.", + ), + ( + "cutoff_type", + "_mutate_cutoff_type", + False, + "answers", + "answers", + "Intended: as cutoff.", + ), + ( + "basis_change", + "_mutate_basis_change", + False, + "answers", + "answers", + "Intended: as cutoff.", + ), + ( + "lower_atol", + "_mutate_lower_atol", + False, + "answers", + "answers", + "Intended: as cutoff.", + ), + ( + "upper_atol", + "_mutate_upper_atol", + False, + "answers", + "answers", + "Intended: as cutoff.", + ), + ) + + @pytest.mark.parametrize("row", ROWS, ids=[row[0] for row in ROWS]) + @pytest.mark.parametrize("partitions", ["off", "auto"]) + @pytest.mark.parametrize("pared", [False, True], ids=["exact", "pared"]) + @pytest.mark.parametrize( + "schrodinger", [False, True], ids=["heisenberg", "schrodinger"] + ) + @pytest.mark.parametrize( + "functional_name", + [ + "expectation_value_functional", + "expectation_value_and_gradient_functional", + ], + ) + def test_mutator_effect_on_live_functional( + self, + monkeypatch, + serial_comm, + functional_name, + schrodinger, + pared, + partitions, + row, + ): + ( + method, + mutator, + needs_empty_graph, + outcome, + pared_schrodinger, + rationale, + ) = row + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + + mp = self._propagator( + serial_comm, schrodinger=schrodinger, with_graph=not needs_empty_graph + ) + parameters = [] if needs_empty_graph else list(self._PARAMS) + threshold = self._PARE_THRESHOLD if pared else None + functional = getattr(mp, functional_name)(threshold) + + before = functional(parameters) + getattr(self, mutator)(mp) + + expected = pared_schrodinger if pared and schrodinger else outcome + context = f"{method}: {rationale}" + if expected == "stale": + with pytest.raises(RuntimeError, match=r"MP object has been modified"): + functional(parameters) + elif expected == "refuses-refresh": + with pytest.raises(RuntimeError, match=r"cannot follow the new weights"): + functional(parameters) + elif expected == "refreshes": + after = functional(parameters) + fresh = getattr(mp, functional_name)(threshold) + _assert_answers_match(after, fresh(parameters), exact=True, context=context) + _assert_answers_differ(after, before, context=context) + else: + _assert_answers_match( + functional(parameters), before, exact=False, context=context + ) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + @pytest.mark.parametrize( + "factory", + [ + "expectation_value_functional", + "expectation_value_and_gradient_functional", + ], + ) + def test_bound_functional_reports_its_parameter_axis( + self, monkeypatch, serial_comm, factory, partitions + ): + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=True) + assert getattr(mp._simulator, factory)(None).num_params == len(self._PARAMS) + # The front end hands back its own callable, which forwards what the engine object exposes. + assert getattr(mp, factory)(None).num_params == len(self._PARAMS) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_parameter_mapping_invalidates_functional_it_desynchronises( + self, monkeypatch, serial_comm, partitions + ): + """The relabel moves the propagator's own answer, and the functional refuses the call rather + than following it half way. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=True) + functional = mp.expectation_value_functional() + parameters = list(self._PARAMS) + before = functional(parameters) + + self._mutate_parameter_mapping(mp) + + assert mp.expval(parameters) != pytest.approx(before) + with pytest.raises(RuntimeError, match=r"set_parameter_mapping"): + functional(parameters) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + @pytest.mark.parametrize( + "factory", + [ + "expectation_value_functional", + "expectation_value_and_gradient_functional", + ], + ) + def test_bound_functional_reports_whether_it_follows_weights( + self, monkeypatch, serial_comm, factory, partitions + ): + """The contract read off the object: only a Schrodinger plan pared against the operator's own + coefficients refuses to follow a re-weight. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + for schrodinger in (False, True): + mp = self._propagator(serial_comm, schrodinger=schrodinger, with_graph=True) + for threshold in (None, self._PARE_THRESHOLD): + follows = not (schrodinger and threshold is not None) + engine = getattr(mp._simulator, factory)(threshold) + assert engine.follows_weights is follows + # The engine object is not part of the public surface, so the rule has to be readable + # off what the front end returns. + assert getattr(mp, factory)(threshold).follows_weights is follows + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + @pytest.mark.parametrize("pared", [False, True], ids=["exact", "pared"]) + @pytest.mark.parametrize( + "functional_name", + [ + "expectation_value_functional", + "expectation_value_and_gradient_functional", + ], + ) + def test_reweighted_functional_matches_a_fresh_propagator( + self, monkeypatch, serial_comm, functional_name, pared, partitions + ): + """The refresh at full precision: a re-weighted propagator's functional must answer what a + propagator built with those coefficients answers, to the last bit. Both run the same + arithmetic over the same store order, so there is no rounding to hide behind. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + threshold = self._PARE_THRESHOLD if pared else None + parameters = list(self._PARAMS) + + mp = self._propagator(serial_comm, schrodinger=False, with_graph=True) + functional = getattr(mp, functional_name)(threshold) + before = functional(parameters) + self._mutate_update_initial_operator(mp) + + fresh = self._propagator( + serial_comm, + schrodinger=False, + with_graph=True, + first_weight=self._REWEIGHTED_FIRST_WEIGHT, + ) + expected = getattr(fresh, functional_name)(threshold)(parameters) + + after = functional(parameters) + _assert_answers_match(after, expected, exact=True) + _assert_answers_differ(after, before) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_pared_schrodinger_functional_refuses_to_follow_a_reweight( + self, monkeypatch, serial_comm, partitions + ): + """Schrodinger thresholds its keep-set from the operator coefficients, so new coefficients + select a different keep-set: the plan says so rather than replaying a paring nobody asked + for. The unpared functional over the same propagator follows the re-weight, which is what + makes the refusal the paring's and not the picture's. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=True, with_graph=True) + parameters = list(self._PARAMS) + pared = mp.expectation_value_functional(self._PARE_THRESHOLD) + exact = mp.expectation_value_functional() + pared(parameters) + + self._mutate_update_initial_operator(mp) + + with pytest.raises(RuntimeError, match=r"cannot follow the new weights"): + pared(parameters) + assert exact(parameters) == pytest.approx(mp.expval(parameters)) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_no_op_mutators_keep_a_functional_valid( + self, monkeypatch, serial_comm, partitions + ): + """A call that appends or folds nothing is not a mutation, on either partition setting. + + The facade decides that before it fans out, so ``off`` and ``auto`` must agree -- an + unconditional bump on the fan-out is exactly the divergence this table exists to rule out. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=False) + functional = mp.expectation_value_functional() + before = functional([]) + + # No graph, so there is nothing to fold and no layer to retire; the return is the current + # coefficients either way. + folded = mp.contract_partially([], inplace=True) + assert folded == pytest.approx(mp.contract_partially([], inplace=False)) + empty = Circuit((), self._MODES, ()) + mp.build_graph(empty) + mp.propagate(empty) + assert mp.n_parameters == 0 + + assert functional([]) == pytest.approx(before) + + # A gate generator is bounds-checked only as the gate loop reaches it, so a bad one in a + # multi-gate call throws with the earlier gates already committed. The front end validates + # generators against num_modes, so the engine is the only way to express one -- as in + # _mutate_basis_change. + _PART_WAY_GATES = ((0,), (2 * _MODES + 1,), (2,)) + _PART_WAY_MAPPING = (0, 1, 2) + _PART_WAY_COEFFS = (1.0, 1.0, 1.0) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_part_way_build_graph_failure_invalidates_the_functional( + self, monkeypatch, serial_comm, partitions + ): + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=True) + functional = mp.expectation_value_functional() + parameters = list(self._PARAMS) + functional(parameters) + layers_before = mp.graph_layers + + with pytest.raises(RuntimeError): + mp._simulator.build_graph( + self._PART_WAY_GATES, self._PART_WAY_MAPPING, self._PART_WAY_COEFFS + ) + + # The graph grew, so the call did mutate on its way to the throw. + assert mp.graph_layers > layers_before + with pytest.raises(RuntimeError, match=r"build_graph"): functional(parameters) - rebuilt = getattr(mp, functional_name)()(parameters) - rebuilt_expval = rebuilt[0] if isinstance(rebuilt, tuple) else rebuilt - assert rebuilt_expval == pytest.approx(mp.expval(parameters)) + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_part_way_propagate_failure_invalidates_the_functional( + self, monkeypatch, serial_comm, partitions + ): + """propagate() folds into the operator instead of appending, so nothing counts the + mutation: without the failure-path bump the functional would keep answering for + coefficients that had moved. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=False) + functional = mp.expectation_value_functional() + functional([]) + + with pytest.raises(RuntimeError): + mp._simulator.propagate( + self._PART_WAY_GATES, + self._PART_WAY_MAPPING, + self._PART_WAY_COEFFS, + [0.3, 0.7, 0.5], + ) + + with pytest.raises(RuntimeError, match=r"propagate"): + functional([]) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + @pytest.mark.parametrize( + "functional_name", + [ + "expectation_value_functional", + "expectation_value_and_gradient_functional", + ], + ) + def test_reweight_that_drops_the_core_term_zeroes_it( + self, monkeypatch, serial_comm, functional_name, partitions + ): + """The identity row describes the dict that committed, like every other row: a re-weight + that leaves it out zeroes it, which is what the store does with every row a dict omits. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + parameters = list(self._PARAMS) + mp = self._propagator( + serial_comm, schrodinger=False, with_graph=True, core_term=0.25 + ) + functional = getattr(mp, functional_name)(None) + before = functional(parameters) + + self._mutate_update_initial_operator(mp) # carries no identity row + + fresh = self._propagator( + serial_comm, + schrodinger=False, + with_graph=True, + first_weight=self._REWEIGHTED_FIRST_WEIGHT, + ) + after = functional(parameters) + _assert_answers_match( + after, getattr(fresh, functional_name)(None)(parameters), exact=True + ) + _assert_answers_match( + after, getattr(mp, _DIRECT_CALL[functional_name])(parameters), exact=False + ) + _assert_answers_differ(after, before) + + @pytest.mark.parametrize("partitions", ["off", "auto"]) + def test_pared_functional_is_invalidated_by_build_graph( + self, monkeypatch, serial_comm, partitions + ): + """A pared plan owns its layers, so its layer count cannot move; the structure revision is + what sees the appended layer. + """ + monkeypatch.setenv("monoprop_PARTITIONS", partitions) + mp = self._propagator(serial_comm, schrodinger=False, with_graph=True) + functional = mp.expectation_value_functional(self._PARE_THRESHOLD) + parameters = list(self._PARAMS) + before = functional(parameters) + + self._mutate_build_graph(mp) + + # The appended gate claims a fresh angle, so the propagator's own axis is now three long. + assert mp.expval([*parameters, 0.4]) != pytest.approx(before) + with pytest.raises(RuntimeError, match=r"build_graph"): + functional(parameters) class TestEvolvedOperatorBothPictures: