fix(propagator): 🐛 invalidate functionals after update_initial_operator

cpp/include/monoprop/MonomialPropagator.h

@@ -308,6 +308,10 @@ class MonomialPropagator {
     std::optional<double> 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_;

cpp/monoprop/Validation.cpp

@@ -26,8 +26,9 @@ class ValidationError : public std::runtime_error {
     using std::runtime_error::runtime_error;
 };
 
-// The graph was rebuilt after a functional captured its layer count, so the functional's parameter
-// mapping no longer describes the graph.
+// 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.
 class StaleFunctionalGraph : public std::runtime_error {
 public:
     using std::runtime_error::runtime_error;
@@ -107,6 +108,14 @@ auto validate_expected_graph_layers(size_t current_layers, size_t expected_layer
     }
 }
 
+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_only_rotate_len_k_(std::optional<size_t> only_rotate_len_k, size_t max_k) -> void {
     if (!only_rotate_len_k.has_value()) {
         return;

cpp/monoprop/Validation.h

@@ -23,8 +23,9 @@
 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 its
-// layer count. Both derive from std::runtime_error, so catching that still catches either.
+// 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.
 
 monoprop_EXPORT auto validate_coefficient_lengths(const VecZ &parameter_mapping, const VecD &gen_coeffs) -> void;
 
@@ -40,6 +41,9 @@ monoprop_EXPORT auto validate_functional_call(const VecD &parameters, size_t exp
 // 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;
+
 // 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<size_t> only_rotate_len_k, size_t max_k) -> void;
 

cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl

@@ -202,6 +202,7 @@ MonomialPropagator<NumModes>::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_),
@@ -350,6 +351,7 @@ auto MonomialPropagator<NumModes>::packed_inline_width_() const -> size_t {
 template <size_t NumModes>
 auto MonomialPropagator<NumModes>::apply_initial_operator_(const OperatorDict &op_dict)
     -> std::pair<MonomialList<NumModes>, VecD> {
+    ++initial_operator_epoch_;
     if (partition_group_) {
         // The facade holds no local terms of its own, so the return is empty.
         for_each_partition_([&](MonomialPropagator &s) { s.update_initial_operator(op_dict); });
@@ -940,6 +942,9 @@ auto MonomialPropagator<NumModes>::make_functional_(Fn &&func, std::optional<dou
     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_;
     const auto &inverted_index = mp_op_.inverted_index();
 
     // One owning handle either way: pare hands back a heap-owned MPGraph the functional must keep alive
@@ -974,9 +979,12 @@ auto MonomialPropagator<NumModes>::make_functional_(Fn &&func, std::optional<dou
             parameter_mapping,
             gen_coeffs,
             num_params,
+            epoch,
+            expected_epoch,
             expected_layers,
             cos = std::move(cos),
             comm](const VecD &params) -> 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,

cpp/tests/update_initial_operator.cpp

@@ -52,6 +52,39 @@ 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);
+}
+
+BOOST_AUTO_TEST_CASE(update_initial_operator_invalidates_gradient_functional) {
+    constexpr size_t n_modes = 2;
+    OperatorDict initial_ham;
+    initial_ham[VecZ{}] = std::complex<double>{1.0, 0.0};
+
+    VecZ initial_state{0, 1};
+    MonomialPropagator<n_modes> simulator(initial_ham,
+                                          2 * n_modes,
+                                          initial_state,
+                                          std::nullopt,
+                                          MPI_COMM_SELF,
+                                          std::nullopt,
+                                          std::nullopt,
+                                          CutoffType::Support,
+                                          std::nullopt);
+
+    const VecD empty_params;
+    auto grad_fn = simulator.expectation_value_and_gradient_functional(std::nullopt);
+    BOOST_TEST(grad_fn(empty_params).first == 1.0, tt::tolerance(1e-12));
+
+    OperatorDict updated;
+    updated[VecZ{}] = std::complex<double>{2.75, 0.0};
+    simulator.update_initial_operator(updated);
+
+    BOOST_CHECK_THROW(grad_fn(empty_params), std::runtime_error);
+    BOOST_TEST(simulator.expectation_value_and_gradient_functional(std::nullopt)(empty_params).first == 2.75,
+               tt::tolerance(1e-12));
 }
 
 BOOST_AUTO_TEST_CASE(update_initial_operator_throws_for_unknown_term_in_heisenberg) {

docs/content/docs/features/evaluation.mdx

@@ -47,6 +47,16 @@ 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.
+
 Both functionals accept an optional `pare_threshold` — see *Paring* below.
 
 ## Paring

src/monoprop/majorana_propagator.py

@@ -127,7 +127,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.
+        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].
 
         Args:
             new_operator: A [MajoranaOperator][monoprop.majorana.MajoranaOperator] whose terms

src/monoprop/monomial_propagator.py

@@ -382,13 +382,21 @@ def expectation_value_functional(
     ) -> Callable[..., float]:
         """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.
+
         Args:
             pare_threshold: Edge-retention cutoff for this functional's masked plan: edges
                 contributing below it are pared away and skipped during replay, trading memory and
                 accuracy for speed. ``None`` (default) disables paring.
 
         Returns:
             A callable ``fn(parameters=None) -> float``.
+
+        Raises:
+            RuntimeError: From the returned callable, if the propagator was mutated after this
+                functional was created.
         """
         fn = self._simulator.expectation_value_functional(pare_threshold)
         return lambda parameters=None: fn(self._bind(parameters))
@@ -398,13 +406,18 @@ def expectation_value_and_gradient_functional(
     ) -> Callable[..., tuple]:
         """Return a reusable callable computing (expectation value, gradient).
 
-        Like [expectation_value_functional][], but one backward pass also yields the gradient.
+        Like [expectation_value_functional][], but one backward pass also yields the gradient. It is
+        invalidated by the same mutations.
 
         Args:
             pare_threshold: See [expectation_value_functional][].
 
         Returns:
             A callable ``fn(parameters=None) -> (float, np.ndarray)``, gradient in parameter order.
+
+        Raises:
+            RuntimeError: From the returned callable, if the propagator was mutated after this
+                functional was created.
         """
         fn = self._simulator.expectation_value_and_gradient_functional(pare_threshold)
 
@@ -528,6 +541,9 @@ 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.
+
         Args:
             new_operator: A [MajoranaOperator][monoprop.majorana.MajoranaOperator] or
                 [PauliOperator][monoprop.pauli.PauliOperator], per the front-end, whose terms replace

src/monoprop/pauli_propagator.py

@@ -220,7 +220,9 @@ def update_initial_operator(self, new_operator: PauliOperator) -> None:
         affected -- the gates and their generator coefficients are unchanged. Unlike the
         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].
+        [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].
 
         Args:
             new_operator: A [PauliOperator][monoprop.pauli.PauliOperator] whose terms replace the

tests/test_parameter_validation.py

@@ -18,8 +18,9 @@
 
 import pytest
 
-from monoprop import Circuit, ExpGate, MajoranaPropagator
+from monoprop import Circuit, ExpGate, MajoranaPropagator, PauliPropagator
 from monoprop.majorana import MajoranaOperator
+from monoprop.pauli import PauliOperator
 
 
 def _two_gate_graph(serial_comm):
@@ -119,6 +120,85 @@ def test_functional_invalidated_after_graph_mutation(self, serial_comm):
         with pytest.raises(RuntimeError, match=r"MP object has been modified"):
             functional([1.0, 2.0])
 
+    @pytest.mark.parametrize(
+        (
+            "propagator_cls",
+            "initial_operator",
+            "updated_operator",
+            "gate_generators",
+            "cutoff",
+        ),
+        [
+            pytest.param(
+                MajoranaPropagator,
+                MajoranaOperator({(0, 1): 1.0j, (2, 3): 0.5j}, num_modes=2),
+                MajoranaOperator({(0, 1): 2.0j, (2, 3): 0.5j}, num_modes=2),
+                (
+                    MajoranaOperator({(0,): 1.0}, num_modes=2),
+                    MajoranaOperator({(1,): 1.0}, num_modes=2),
+                ),
+                4,
+                id="majorana",
+            ),
+            pytest.param(
+                PauliPropagator,
+                PauliOperator({"ZZ": 1.0, "XX": 0.5}, num_qubits=2),
+                PauliOperator({"ZZ": 2.0, "XX": 0.5}, num_qubits=2),
+                (
+                    PauliOperator({"XI": 1.0}, num_qubits=2),
+                    PauliOperator({"IY": 1.0}, num_qubits=2),
+                ),
+                2,
+                id="pauli",
+            ),
+        ],
+    )
+    @pytest.mark.parametrize(
+        "schrodinger_cutoff", [None, 2], ids=["heisenberg", "schrodinger"]
+    )
+    @pytest.mark.parametrize(
+        "functional_name",
+        [
+            "expectation_value_functional",
+            "expectation_value_and_gradient_functional",
+        ],
+    )
+    def test_functional_invalidated_after_initial_operator_update(
+        self,
+        serial_comm,
+        propagator_cls,
+        initial_operator,
+        updated_operator,
+        gate_generators,
+        cutoff,
+        schrodinger_cutoff,
+        functional_name,
+    ):
+        mp = propagator_cls(
+            initial_operator,
+            [0, 1],
+            cutoff=cutoff,
+            schrodinger_cutoff=schrodinger_cutoff,
+            comm=serial_comm,
+        )
+        mp.build_graph(
+            Circuit(tuple(ExpGate(generator) for generator in gate_generators), 2)
+        )
+        functional = getattr(mp, functional_name)()
+        parameters = [0.3, 0.7]
+        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"):
+            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))
+
 
 class TestEvolvedOperatorBothPictures:
     def test_schrodinger_returns_state_dict(self, serial_comm):