From 951af79b62262ae8f0594b1cd0efbe55bd55aa6e Mon Sep 17 00:00:00 2001 From: Josh Cudby Date: Thu, 3 Sep 2026 14:25:07 +0000 Subject: [PATCH 1/2] =?UTF-8?q?feat(propagator):=20=E2=9C=A8=20read=20indi?= =?UTF-8?q?vidual=20evolved=20coefficients=20by=20term?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `evolved_operator` decodes the whole evolved operator: it enumerates the operator index, materialises a key per surviving term, and builds a term dictionary. `evolved_operator_coefficients` instead probes the index with the terms the caller names, so the decode costs one entry per term *requested* rather than one per term the operator carries. Reading a few amplitudes out of an evolved state (Schrodinger picture) is the motivating case. The graph contraction is unchanged, and still the dominant cost. Contracts: a term the operator does not carry reads back as 0; there is no atol, since the caller named the terms it wants and magnitude filtering would silently zero some of them; the result is positional; and it is rank-local, exactly as `evolved_operator` is. Terms must be canonical monomials, as the constructor's initial_operator keys must be. The engine keys terms by an order-insensitive bitset, so a raw `(1, 0)` would resolve to the row of `(0, 1)` and read back its coefficient without the anticommutation sign -- which a lookup has no coefficient of its own to carry. The Majorana front-end therefore validates raw sequences through `Majorana` rather than normalizing them; use `Majorana.from_unsorted` and apply the sign it returns. Under partitions each partition probes on its own master thread into its own vector, and the merge runs on the facade thread, matching how every other partition fan-out in the file collects its results. Assisted-by: Claude Opus 5 (1M context) --- cpp/include/monoprop/MonomialPropagator.h | 15 + .../MonomialPropagator.inl | 58 +++ .../evolved_operator_coefficients_tests.cpp | 323 ++++++++++++++++ docs/content/docs/features/evaluation.mdx | 54 +++ src/monoprop/bindings/binder.h | 6 + src/monoprop/majorana_propagator.py | 40 +- src/monoprop/monomial_propagator.py | 85 ++++- src/monoprop/pauli_propagator.py | 25 +- tests/test_evolved_operator_coefficients.py | 351 ++++++++++++++++++ tests/test_monoprop_mpi.py | 41 ++ 10 files changed, 991 insertions(+), 7 deletions(-) create mode 100644 cpp/tests/evolved_operator_coefficients_tests.cpp create mode 100644 tests/test_evolved_operator_coefficients.py diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index b684ff58..914216ba 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -283,6 +283,21 @@ class MonomialPropagator { auto evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>>; + /// Coefficients of the requested monomials only, in query order; a term the operator does not + /// carry yields 0. The same contraction and decode as evolved_operator_terms(), but the index is + /// probed with the caller's keys rather than enumerated, so nothing is materialised per surviving + /// term. No atol: the caller named the terms it wants, so magnitude filtering would silently zero + /// some of them. The empty term resolves to core_term() in the Heisenberg picture, where the + /// identity is held apart from the index, and through the probe in the Schrodinger picture, where + /// it is an ordinary index row (the state's identity amplitude) -- so on every key this agrees + /// with what evolved_operator_terms()'s callers publish, the core term they re-add included. + /// Keys must be canonical monomials, as the constructor's initial_operator keys must be: the + /// encode is order-insensitive, so a repeated or unsorted index would resolve to the canonical + /// term's row and lose the sign of the reordering. + /// Contracts non-inplace. Rank-local, as evolved_operator_terms() is. + auto evolved_operator_coefficients(const VecD ¶meters, const std::vector &terms) + -> std::vector>; + virtual auto update_initial_operator(const OperatorDict &op_dict) -> void { apply_initial_operator_(op_dict); } protected: diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 34814906..875bf0b8 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -1133,4 +1133,62 @@ auto MonomialPropagator::evolved_operator_terms(const VecD ¶meters return concat_partitions_(collect); } +template +auto MonomialPropagator::evolved_operator_coefficients(const VecD ¶meters, const std::vector &terms) + -> std::vector> { + std::vector> keys; + keys.reserve(terms.size()); + for (const auto &term : terms) { + // Checked: these are user-supplied indices, and the unchecked encode would write out of bounds. + keys.push_back(indices_to_bitset_checked(term, 2 * logical_num_modes_)); + } + + // `p` is always unpartitioned here (a partition, or *this), so indexing() is available. + const auto probe = [&](MonomialPropagator &p) -> std::vector> { + const VecD evolved = p.contract_partially(parameters, false); + std::vector rows(keys.size()); + p.indexing().find_batch(keys.data(), keys.size(), rows.data()); + std::vector> found(keys.size()); + for (size_t q = 0; q < keys.size(); ++q) { + if (rows[q] >= evolved.size()) { // kNotFound is size_t max, so this covers a miss too + continue; + } + found[q] = algebra_decode_coeff(basis_, evolved[rows[q]], keys[q]); + } + return found; + }; + + std::vector> out(terms.size()); + if (!partition_group_) { + out = probe(*this); + } + else { + // The partitions probe concurrently on their own master threads, so each fills its own + // vector and the merge happens here. They are disjoint, so at most one contributes per key. + for (const auto &partition_found : map_partitions_(probe)) { + for (size_t q = 0; q < out.size(); ++q) { + out[q] += partition_found[q]; + } + } + } + + // Match evolved_operator_terms(): round off the anti-hermitian numerical noise. + for (auto &coeff : out) { + coeff = {std::round(coeff.real() * 1e12) / 1e12, std::round(coeff.imag() * 1e12) / 1e12}; + } + + // Heisenberg only: the empty monomial is diverted to core_term_ instead of indexed, so the probe + // cannot have found it. Assigned after the rounding, so the value is bit-for-bit the core_term() + // that evolved_operator's callers publish under the empty key. + if (!schrodinger_) { + const auto core = core_term(); + for (size_t q = 0; q < terms.size(); ++q) { + if (terms[q].empty()) { + out[q] = {core, 0.0}; + } + } + } + return out; +} + } // namespace monoprop diff --git a/cpp/tests/evolved_operator_coefficients_tests.cpp b/cpp/tests/evolved_operator_coefficients_tests.cpp new file mode 100644 index 00000000..ef5b5af2 --- /dev/null +++ b/cpp/tests/evolved_operator_coefficients_tests.cpp @@ -0,0 +1,323 @@ +// 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. + +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "TestUtilities.h" +#include "monoprop/MonomialPropagator.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/detail/mpi/MPICompat.h" + +// evolved_operator_coefficients() probes the operator index with caller-supplied keys instead of +// enumerating it. Oracle throughout: evolved_operator_terms(params, 0.0), which enumerates. + +namespace { + +using namespace monoprop; +using namespace test_utils; + +constexpr size_t kNumModes = 8; +constexpr unsigned int kCutoff = 8; + +auto majorana_sim(const CaseData &data, size_t partitions = 1) -> MonomialPropagator { + return MonomialPropagator(data.hamiltonian, + kCutoff, + data.initial_state, + std::nullopt, + MPI_COMM_SELF, + std::nullopt, + std::nullopt, + CutoffType::Length, + std::nullopt, + kNumModes, + Basis::Majorana, + partitions); +} + +// A propagator whose graph is built and ready to contract. +auto built_sim(const CaseData &data, size_t partitions = 1) -> MonomialPropagator { + auto sim = majorana_sim(data, partitions); + sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); + return sim; +} + +// A Schrodinger-picture propagator: the state's identity amplitude is an ordinary index row there, +// so the empty key resolves through the probe rather than through core_term(). +auto built_schrodinger_sim(const CaseData &data, size_t partitions = 1) -> MonomialPropagator { + auto sim = MonomialPropagator(data.hamiltonian, + kCutoff, + data.initial_state, + std::optional{kCutoff}, + MPI_COMM_SELF, + std::nullopt, + std::nullopt, + CutoffType::Length, + std::nullopt, + kNumModes, + Basis::Majorana, + partitions); + sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); + return sim; +} + +auto keys_of(const std::vector>> &terms) -> std::vector { + std::vector keys; + keys.reserve(terms.size()); + for (const auto &[indices, coeff] : terms) { + keys.push_back(indices); + } + return keys; +} + +} // namespace + +// Querying every enumerated term must reproduce the enumerated coefficients exactly: same +// contraction, same decode, only the index access differs. +BOOST_AUTO_TEST_CASE(coefficients_match_enumerated_terms) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_sim(data); + + const auto terms = sim.evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(!terms.empty()); + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, keys_of(terms)); + BOOST_REQUIRE_EQUAL(coeffs.size(), terms.size()); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST_CONTEXT("term " << i) { + BOOST_TEST(coeffs[i].real() == terms[i].second.real()); + BOOST_TEST(coeffs[i].imag() == terms[i].second.imag()); + } + } +} + +// A monomial the operator does not carry reads back as exactly 0, not as noise and not as a throw. +BOOST_AUTO_TEST_CASE(absent_term_reads_back_zero) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_sim(data); + + const auto terms = sim.evolved_operator_terms(data.parameters, 0.0); + std::map> present(terms.begin(), terms.end()); + + // Search the weight-2 monomials for one the evolved operator misses; the cutoff makes some + // absent for any realistic case, but assert we actually found one rather than trusting it. + std::optional absent; + for (size_t i = 0; i < 2 * kNumModes && !absent; ++i) { + for (size_t j = i + 1; j < 2 * kNumModes; ++j) { + const VecZ candidate{i, j}; + if (!present.contains(candidate)) { + absent = candidate; + break; + } + } + } + BOOST_REQUIRE_MESSAGE(absent.has_value(), "no absent weight-2 monomial to probe"); + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, {*absent}); + BOOST_REQUIRE_EQUAL(coeffs.size(), 1u); + BOOST_TEST(coeffs[0].real() == 0.0); + BOOST_TEST(coeffs[0].imag() == 0.0); +} + +// The result is positional: out[q] belongs to keys[q] whatever order the keys arrive in. +BOOST_AUTO_TEST_CASE(query_order_is_preserved) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_sim(data); + + const auto terms = sim.evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(terms.size() > 4); + + std::vector order(terms.size()); + std::iota(order.begin(), order.end(), 0u); + std::mt19937 rng(1234); + std::shuffle(order.begin(), order.end(), rng); + + std::vector keys; + keys.reserve(order.size()); + for (const auto &i : order) { + keys.push_back(terms[i].first); + } + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, keys); + BOOST_REQUIRE_EQUAL(coeffs.size(), order.size()); + for (size_t q = 0; q < order.size(); ++q) { + BOOST_TEST_CONTEXT("query " << q << " -> term " << order[q]) { + BOOST_TEST(coeffs[q].real() == terms[order[q]].second.real()); + BOOST_TEST(coeffs[q].imag() == terms[order[q]].second.imag()); + } + } +} + +// Partitions hash-split the operator, so a key resolves in exactly one of them. The probe must +// still see the whole operator, matching the single-partition run term for term. +BOOST_AUTO_TEST_CASE(partitioned_probe_matches_single_partition) { + const auto data = load_case_data("random_exact.msgpack"); + auto ref = built_sim(data); + + const auto terms = ref.evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(!terms.empty()); + const auto keys = keys_of(terms); + const auto expected = ref.evolved_operator_coefficients(data.parameters, keys); + + for (const size_t partitions : {size_t{2}, size_t{4}}) { + auto sim = built_sim(data, partitions); + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, keys); + BOOST_REQUIRE_EQUAL(coeffs.size(), expected.size()); + for (size_t q = 0; q < expected.size(); ++q) { + BOOST_TEST_CONTEXT("partitions=" << partitions << " query " << q) { + BOOST_TEST(near(coeffs[q].real(), expected[q].real())); + BOOST_TEST(near(coeffs[q].imag(), expected[q].imag())); + } + } + } +} + +// The keys are user input, so the checked encode is mandatory: an out-of-range slot must throw +// rather than write past the monomial (indices_to_bitset is noexcept and would corrupt memory). +BOOST_AUTO_TEST_CASE(out_of_range_slot_index_throws) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_sim(data); + + const VecZ too_large{2 * kNumModes}; + BOOST_CHECK_THROW(sim.evolved_operator_coefficients(data.parameters, {too_large}), AlgebraIndexOutOfRange); +} + +// An empty query list is a no-op, not a degenerate probe. +BOOST_AUTO_TEST_CASE(empty_query_list_returns_empty) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_sim(data); + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, {}); + BOOST_TEST(coeffs.empty()); +} + +// evolved_operator_terms() excludes the core term and the binder re-adds it under the empty key, so +// the empty query has to resolve to core_term() for the two APIs to agree term for term. The +// fixture Hamiltonian has no identity term, so inject one -- otherwise core_term() is 0 and the +// check cannot tell the core term apart from the absent-term answer. +BOOST_AUTO_TEST_CASE(empty_term_resolves_to_core_term) { + auto data = load_case_data("random_exact.msgpack"); + data.hamiltonian[VecZ{}] = std::complex{0.75, 0.0}; + auto sim = built_sim(data); + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, {VecZ{}}); + BOOST_REQUIRE_EQUAL(coeffs.size(), 1u); + BOOST_TEST(coeffs[0].real() == 0.75); + BOOST_TEST(coeffs[0].imag() == 0.0); + BOOST_TEST(coeffs[0].real() == sim.core_term()); + + // The identity is not an index row in the Heisenberg picture, so the enumeration never sees it + // and the two APIs would otherwise disagree on exactly this key. + const auto terms = sim.evolved_operator_terms(data.parameters, 0.0); + BOOST_TEST(std::none_of(terms.begin(), terms.end(), [](const auto &term) { return term.first.empty(); })); +} + +// The Schrodinger picture is the motivating case: reading a handful of amplitudes out of an evolved +// state. There is no core term to divert the identity there -- it is an ordinary index row carrying +// the state's identity amplitude -- so the empty query resolves through the probe like any other. +BOOST_AUTO_TEST_CASE(schrodinger_coefficients_match_enumerated_terms) { + const auto data = load_case_data("random_exact.msgpack"); + auto sim = built_schrodinger_sim(data); + + const auto terms = sim.evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(!terms.empty()); + + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, keys_of(terms)); + BOOST_REQUIRE_EQUAL(coeffs.size(), terms.size()); + for (size_t i = 0; i < terms.size(); ++i) { + BOOST_TEST_CONTEXT("term " << i) { + BOOST_TEST(coeffs[i].real() == terms[i].second.real()); + BOOST_TEST(coeffs[i].imag() == terms[i].second.imag()); + } + } + + // The enumeration above covers the identity only because the state carries it; pin that, so a + // regression that started diverting it to the core term would be caught here. + const auto identity = std::find_if(terms.begin(), terms.end(), [](const auto &term) { return term.first.empty(); }); + BOOST_REQUIRE_MESSAGE(identity != terms.end(), "the evolved state has no identity amplitude to pin"); + const auto probed = sim.evolved_operator_coefficients(data.parameters, {VecZ{}}); + BOOST_TEST(probed.at(0).real() == identity->second.real()); + BOOST_TEST(probed.at(0).imag() == identity->second.imag()); +} + +// The two axes above cross: this is the only configuration where the identity is a hash-partitioned +// index row, so it exercises the empty key against the concurrent per-partition probe and merge. +BOOST_AUTO_TEST_CASE(partitioned_schrodinger_probe_matches_single_partition) { + const auto data = load_case_data("random_exact.msgpack"); + auto ref = built_schrodinger_sim(data); + + const auto terms = ref.evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(!terms.empty()); + const auto keys = keys_of(terms); + BOOST_REQUIRE(std::any_of(keys.begin(), keys.end(), [](const auto &key) { return key.empty(); })); + const auto expected = ref.evolved_operator_coefficients(data.parameters, keys); + + for (const size_t partitions : {size_t{2}, size_t{4}}) { + auto sim = built_schrodinger_sim(data, partitions); + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, keys); + BOOST_REQUIRE_EQUAL(coeffs.size(), expected.size()); + for (size_t q = 0; q < expected.size(); ++q) { + BOOST_TEST_CONTEXT("partitions=" << partitions << " query " << q) { + BOOST_TEST(near(coeffs[q].real(), expected[q].real())); + BOOST_TEST(near(coeffs[q].imag(), expected[q].imag())); + } + } + } +} + +// core_term_ is replicated on every partition, so the empty key must be assigned from core_term() +// -- routed through partition 0 -- and not summed over the partitions the way an index row is. +BOOST_AUTO_TEST_CASE(partitioned_empty_term_resolves_to_core_term) { + auto data = load_case_data("random_exact.msgpack"); + data.hamiltonian[VecZ{}] = std::complex{0.75, 0.0}; + + for (const size_t partitions : {size_t{1}, size_t{2}, size_t{4}}) { + auto sim = built_sim(data, partitions); + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, {VecZ{}}); + BOOST_TEST_CONTEXT("partitions=" << partitions) { + BOOST_REQUIRE_EQUAL(coeffs.size(), 1u); + BOOST_TEST(coeffs[0].real() == 0.75); + BOOST_TEST(coeffs[0].imag() == 0.0); + } + } +} + +// The query is a list, not a set: a repeated key is answered once per occurrence. Worth pinning +// because the partitioned merge accumulates, so a key counted twice would show up here first. +BOOST_AUTO_TEST_CASE(repeated_key_is_answered_once_per_occurrence) { + const auto data = load_case_data("random_exact.msgpack"); + + const auto terms = built_sim(data).evolved_operator_terms(data.parameters, 0.0); + BOOST_REQUIRE(!terms.empty()); + const auto &key = terms.front().first; + const auto expected = terms.front().second; + + for (const size_t partitions : {size_t{1}, size_t{2}}) { + auto sim = built_sim(data, partitions); + const auto coeffs = sim.evolved_operator_coefficients(data.parameters, {key, VecZ{}, key}); + BOOST_TEST_CONTEXT("partitions=" << partitions) { + BOOST_REQUIRE_EQUAL(coeffs.size(), 3u); + BOOST_TEST(near(coeffs[0].real(), expected.real())); + BOOST_TEST(near(coeffs[0].imag(), expected.imag())); + BOOST_TEST(coeffs[2] == coeffs[0]); + } + } +} diff --git a/docs/content/docs/features/evaluation.mdx b/docs/content/docs/features/evaluation.mdx index 470854b4..0d3dbd31 100644 --- a/docs/content/docs/features/evaluation.mdx +++ b/docs/content/docs/features/evaluation.mdx @@ -102,3 +102,57 @@ To read the fully evolved operator as a [MajoranaOperator][] (or [PauliOperator] [PauliPropagator][]) — without modifying the simulator — use [MajoranaPropagator.evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] (or [PauliPropagator.evolved_operator][pauli_propagator.PauliPropagator.evolved_operator]). + +## Reading individual terms + +[MajoranaPropagator.evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] +decodes the *whole* evolved operator: it enumerates the operator index, +materialises a key per surviving term, and hands back a term dictionary. When +only a handful of terms are wanted — reading a few amplitudes out of an evolved +state in the Schrödinger picture, say — use +[evolved_operator_coefficients][monoprop.monomial_propagator.MonomialPropagator.evolved_operator_coefficients] +instead. It *probes* the index with the terms you name, so the work scales with +the terms you asked for rather than with the terms the operator carries: + +```python notest +# The coefficients of these three terms alone, in this order. +terms = [Majorana(0, 1), Majorana(2, 3), Majorana(0, 3)] +coeffs = sim.evolved_operator_coefficients(terms, parameters) +``` + +The terms are given in the front-end's own vocabulary — [Majorana][] terms (or +raw index tuples) for [MajoranaPropagator][], [Pauli][] terms for +[PauliPropagator][] — and the returned NumPy array carries one coefficient per +requested term, in the order requested. + +Three contracts are worth knowing: + +- **A term the evolved operator does not carry reads back as `0`.** There is + deliberately no `atol`: you named the terms you want, so filtering by + magnitude would silently zero some of them. Threshold the returned array + instead. +- **Terms must be canonical, and a non-canonical one raises.** The engine keys + its terms by an order-insensitive bitset, so a raw `(1, 0)` would resolve to + the row of `(0, 1)` — and reordering Majorana indices anticommutes them, so + its coefficient carries a sign that a lookup has nowhere to put. Raw tuples + are therefore validated, not normalized: sort with + [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply + the sign it hands back yourself. The `terms` keys of a + [MajoranaOperator][] are already canonical, as are the qubits of a + [Pauli][]. +- **It is rank-local, exactly as + [evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] + is.** Under MPI each rank answers only for the terms in its own hash + partition, and a term another rank owns reads back as `0`; sum across the + communicator to reassemble the operator. + +The empty (identity) term agrees with +[evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] +either way: it is the core term in the Heisenberg picture, where the identity is +held apart from the operator index, and the state's identity amplitude in the +Schrödinger picture, where it is an ordinary term. Repeating a term in one query +simply answers it twice. + +This does **not** make the evolution faster. The graph contraction is identical +and still the dominant cost; what disappears is the decode and dictionary build +over every surviving term. diff --git a/src/monoprop/bindings/binder.h b/src/monoprop/bindings/binder.h index 4076c7ff..cc571d90 100644 --- a/src/monoprop/bindings/binder.h +++ b/src/monoprop/bindings/binder.h @@ -219,6 +219,12 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { "atol"_a, "The evolved operator as a {indices: coefficient} dict, keeping terms with |coeff| >= atol"); + cls.def("evolved_operator_coefficients", + &MonomialPropagator::evolved_operator_coefficients, + "parameters"_a, + "terms"_a, + "Coefficients of the given index terms, in query order; terms the operator does not carry are 0"); + cls.def_prop_ro("num_modes", &MonomialPropagator::logical_num_modes, "Number of modes the operator actually uses"); diff --git a/src/monoprop/majorana_propagator.py b/src/monoprop/majorana_propagator.py index 477b70fd..404e5c57 100644 --- a/src/monoprop/majorana_propagator.py +++ b/src/monoprop/majorana_propagator.py @@ -16,19 +16,18 @@ from __future__ import annotations +from collections.abc import Iterable, Sequence from typing import TYPE_CHECKING -from .majorana import MajoranaOperator +from .majorana import Majorana, MajoranaOperator from .monomial_propagator import MonomialPropagator if TYPE_CHECKING: - from collections.abc import Sequence - import numpy as np from mpi4py import MPI from .circuit import Circuit, ExpGate - from .monomial_propagator import ParameterValues + from .monomial_propagator import OperatorTerm, ParameterValues from .quantum_data import IQuantumOperator @@ -122,6 +121,39 @@ def evolved_operator( terms = self._simulator.evolved_operator(self._bind(parameters), atol) return MajoranaOperator(terms, self.num_modes) + def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: + """Encode a Majorana term into the engine's index tuple. + + Majorana indices *are* the engine's keys, so this only unwraps the term. A raw sequence is + round-tripped through [Majorana][monoprop.majorana.Majorana] to reject a non-canonical + product: the engine's encode is order-insensitive, so ``(1, 0)`` would resolve to the row + of ``(0, 1)`` and read back its coefficient without the anticommutation sign, which this + cannot apply to a coefficient it has yet to look up. Canonicalize first with + [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply the sign it + returns yourself; the ``terms`` keys of a + [MajoranaOperator][monoprop.majorana.MajoranaOperator] are already canonical. + + Args: + term: A [Majorana][monoprop.majorana.Majorana] term, or a raw index sequence + (a tuple, or an index array). + + Returns: + The term's Majorana indices. + + Raises: + TypeError: If the term is neither a Majorana term nor an iterable of indices. + ValueError: If a raw sequence is not sorted, or holds a negative or repeated index. + """ + if isinstance(term, Majorana): + return term.indices + # Iterable rather than Sequence: a NumPy array of indices is not a Sequence. + if not isinstance(term, Iterable): + raise TypeError( + "Majorana terms are Majorana objects or index sequences; got " + f"{type(term).__name__}." + ) + return Majorana(*term).indices + def update_initial_operator(self, new_operator: MajoranaOperator) -> None: """Replace the *initial operator* (existing terms only). diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py index 63be1c0b..9b1677ed 100644 --- a/src/monoprop/monomial_propagator.py +++ b/src/monoprop/monomial_propagator.py @@ -43,12 +43,18 @@ from .utils import validate_basis_change if TYPE_CHECKING: - from collections.abc import Callable, Sequence + from collections.abc import Callable, Iterable, Sequence from typing import Self from mpi4py import MPI + from .majorana import Majorana + from .pauli import Pauli + ParameterValues = Circuit | Sequence[float] | np.ndarray | None + # A single operator term in a front-end's own vocabulary: a Pauli, a Majorana, or the raw + # index sequence either engine keys its terms by. + OperatorTerm = Majorana | Pauli | Sequence[int] | np.ndarray logger = logging.getLogger(__name__) @@ -533,6 +539,83 @@ def evolved_operator( The evolved operator (Heisenberg picture) or evolved state (Schrodinger picture). """ + def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: + """Encode one operator term into the raw index tuple the engine keys terms by. + + The front-end counterpart to the decode ``evolved_operator`` performs. Not abstract, so + that a front-end which has no term encoding still constructs; it only has to provide this + to support [evolved_operator_coefficients][]. + + Terms must be canonical -- the engine's encode is order-insensitive, so a non-canonical + product would resolve to the canonical term's row and lose the sign of the reordering -- so + an implementation validates rather than normalizes: it has no coefficient to put the sign + on yet. + + Args: + term: A single term in this front-end's own vocabulary. + + Returns: + The term's engine index tuple (Majorana indices, or symplectic slots in the Pauli + basis). + + Raises: + NotImplementedError: If this front-end does not encode terms. + """ + raise NotImplementedError( + f"{type(self).__name__} does not implement _term_slots, so it cannot look up " + "individual evolved coefficients." + ) + + def evolved_operator_coefficients( + self, + terms: Iterable[OperatorTerm], + parameters: ParameterValues = None, + ) -> np.ndarray: + """Return the coefficients of ``terms`` alone in the evolved operator, in the order given. + + A cheaper + [evolved_operator][monoprop.monomial_propagator.MonomialPropagator.evolved_operator] + for when only a few terms are wanted. The graph contraction is identical -- and still the + dominant cost -- but the operator index is *probed* with these terms rather than + enumerated, so the decode and the term dictionary cost one entry per term *requested* + instead of one per term the evolved operator carries. Reading a handful of amplitudes out + of an evolved state (Schrodinger picture) is the motivating case. + + A term the evolved operator does not carry reads back as ``0``. There is deliberately no + ``atol``: the caller named the terms it wants, so magnitude filtering would silently zero + some of them -- threshold the returned array instead. + + Terms must be canonical, and a non-canonical one is rejected rather than normalized: the + anticommutation sign of the reordering belongs on the coefficient, which is what this is + about to look up. Repeating a term simply answers it twice. + + The empty (identity) term agrees with ``evolved_operator`` too: it is the core term in the + Heisenberg picture, where the identity is held apart from the operator index, and the + state's identity amplitude in the Schrodinger picture, where it is an ordinary term. + + Rank-local, as ``evolved_operator`` is: under MPI each rank answers for the terms it owns, + and a term another rank owns reads back as ``0``. + + Args: + terms: The terms to look up, in this front-end's own vocabulary. + parameters: Variational parameter values (see [expectation_value][]). + + Returns: + A complex NumPy array, one coefficient per requested term, in the order requested. + + Raises: + TypeError: If a term is not in this front-end's vocabulary. + ValueError: If a term is not a canonical monomial. + RuntimeError: If a term index lies outside the propagator's own system. + """ + slots = [self._term_slots(term) for term in terms] + return np.asarray( + self._simulator.evolved_operator_coefficients( + self._bind(parameters), slots + ), + dtype=complex, + ) + @abstractmethod def update_initial_operator(self, new_operator: T_op) -> None: """Replace the *initial operator* (existing terms only). diff --git a/src/monoprop/pauli_propagator.py b/src/monoprop/pauli_propagator.py index 1b0ec4e8..be345229 100644 --- a/src/monoprop/pauli_propagator.py +++ b/src/monoprop/pauli_propagator.py @@ -18,7 +18,7 @@ from typing import TYPE_CHECKING -from .conversion_utils import _local_slots_to_pauli +from .conversion_utils import _local_slots_to_pauli, _pauli_to_local_slots from .monomial_propagator import MonomialPropagator from .pauli import Pauli, PauliOperator @@ -29,7 +29,7 @@ from mpi4py import MPI from .circuit import Circuit, ExpGate - from .monomial_propagator import ParameterValues + from .monomial_propagator import OperatorTerm, ParameterValues class PauliPropagator(MonomialPropagator[PauliOperator]): @@ -117,6 +117,27 @@ def evolved_operator( } return PauliOperator(terms, self.num_qubits) + def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: + """Encode a qubit Pauli term into the engine's symplectic slots. + + Only [Pauli][monoprop.pauli.Pauli] terms are accepted: the slots are an engine-internal + encoding, so a raw slot sequence is rejected rather than passed through. + + Args: + term: A [Pauli][monoprop.pauli.Pauli] term. + + Returns: + The term's symplectic slot indices. + + Raises: + TypeError: If the term is not a [Pauli][monoprop.pauli.Pauli]. + """ + if not isinstance(term, Pauli): + raise TypeError( + f"Pauli terms are Pauli objects; got {type(term).__name__}." + ) + return _pauli_to_local_slots(term.string, term.qubits) + def _circuit_gates(self, circuit: Circuit) -> Sequence[ExpGate]: """Accept a qubit circuit; its gates are expanded by the shared pipeline. diff --git a/tests/test_evolved_operator_coefficients.py b/tests/test_evolved_operator_coefficients.py new file mode 100644 index 00000000..99d84886 --- /dev/null +++ b/tests/test_evolved_operator_coefficients.py @@ -0,0 +1,351 @@ +# 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. + +"""Coverage for ``evolved_operator_coefficients``, the term-probing companion to ``evolved_operator``. + +Oracle throughout: ``evolved_operator(atol=0.0)``, which enumerates the whole index. Both are +rank-local, so every test takes ``serial_comm`` (see the note in ``test_basis.py``). +""" + +from __future__ import annotations + +from pathlib import Path + +import numpy as np +import pytest + +from monoprop import ( + Circuit, + ExpGate, + MajoranaPropagator, + PauliPropagator, +) +from monoprop.majorana import Majorana, MajoranaOperator +from monoprop.monomial_propagator import MonomialPropagator +from monoprop.pauli import Pauli, PauliOperator +from tests.cases import load_problem + +DATA = Path(__file__).parent / "data" + +N_QUBITS = 6 + + +def _majorana_propagator(problem, serial_comm, schrodinger_cutoff=None): + prop = MajoranaPropagator( + problem.operator, + problem.monomial_circuit.initial_state, + cutoff=2 * problem.n_modes, + schrodinger_cutoff=schrodinger_cutoff, + comm=serial_comm, + ) + prop.build_graph(problem.monomial_circuit.to_circuit()) + return prop + + +def _pauli_propagator(serial_comm, schrodinger_cutoff=None): + """A Pauli propagator with a graph whose evolution spreads the operator over many terms.""" + prop = PauliPropagator( + PauliOperator({Pauli("ZZ", (0, 1)): 1.0, Pauli("XY", (2, 3)): 0.5}, N_QUBITS), + initial_state=[], + cutoff=N_QUBITS, + schrodinger_cutoff=schrodinger_cutoff, + comm=serial_comm, + ) + circuit = Circuit( + tuple( + ExpGate(PauliOperator({Pauli(letters, qubits): 1.0}, N_QUBITS)) + for letters, qubits in ( + ("XX", (0, 1)), + ("YZ", (1, 2)), + ("ZX", (2, 3)), + ("XY", (3, 4)), + ("YY", (4, 5)), + ) + ), + N_QUBITS, + ) + prop.build_graph(circuit) + return prop, np.linspace(0.1, 0.9, prop.n_parameters) + + +def test_pauli_coefficients_match_evolved_operator(serial_comm) -> None: + """Every term the evolved Pauli operator carries reads back with the same coefficient.""" + prop, parameters = _pauli_propagator(serial_comm) + + evolved = prop.evolved_operator(parameters, atol=0.0) + terms = list(evolved.terms) + assert len(terms) > 1 + + coefficients = prop.evolved_operator_coefficients(terms, parameters) + + assert coefficients.shape == (len(terms),) + np.testing.assert_allclose( + coefficients, [evolved.terms[term] for term in terms], atol=1e-12 + ) + + +def test_majorana_coefficients_match_evolved_operator(serial_comm) -> None: + """Same property for the Majorana front-end, whose keys are raw index tuples.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + terms = list(evolved.terms) + assert len(terms) > 1 + + coefficients = prop.evolved_operator_coefficients(terms, parameters) + + np.testing.assert_allclose( + coefficients, [evolved.terms[term] for term in terms], atol=1e-12 + ) + + +def test_majorana_accepts_majorana_terms(serial_comm) -> None: + """A Majorana term object keys the same coefficient its raw index tuple does.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + terms = list(evolved.terms)[:4] + + from_tuples = prop.evolved_operator_coefficients(terms, parameters) + from_objects = prop.evolved_operator_coefficients( + [Majorana(*term) for term in terms], parameters + ) + + np.testing.assert_array_equal(from_tuples, from_objects) + + +def test_absent_term_reads_back_zero(serial_comm) -> None: + """A term the evolved operator does not carry is 0, not a raise and not noise.""" + prop = PauliPropagator( + PauliOperator({Pauli("Z", (0,)): 1.0}, N_QUBITS), + initial_state=[], + cutoff=N_QUBITS, + comm=serial_comm, + ) + + absent = Pauli("X", (5,)) + assert absent not in prop.evolved_operator(atol=0.0).terms + + coefficients = prop.evolved_operator_coefficients([absent]) + + assert coefficients.shape == (1,) + assert coefficients[0] == 0 + + +def test_schrodinger_coefficients_match_evolved_state(serial_comm) -> None: + """The motivating case: reading a few amplitudes out of an evolved state.""" + prop, parameters = _pauli_propagator(serial_comm, schrodinger_cutoff=4) + + evolved = prop.evolved_operator(parameters, atol=0.0) + terms = list(evolved.terms) + assert len(terms) > 1 + + coefficients = prop.evolved_operator_coefficients(terms, parameters) + + np.testing.assert_allclose( + coefficients, [evolved.terms[term] for term in terms], atol=1e-12 + ) + + +def test_query_order_is_preserved(serial_comm) -> None: + """``out[i]`` belongs to ``terms[i]`` whatever order the terms arrive in. + + Majorana rather than Pauli so the shuffled query spans many of the engine's batch-probe + prefetch groups, not just the first. + """ + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + terms = list(evolved.terms) + assert len(terms) > 100 + + shuffled = list(np.random.default_rng(0).permutation(len(terms))) + reordered = [terms[i] for i in shuffled] + + coefficients = prop.evolved_operator_coefficients(reordered, parameters) + + np.testing.assert_allclose( + coefficients, [evolved.terms[term] for term in reordered], atol=1e-12 + ) + + +def test_identity_term_is_the_core_term(serial_comm) -> None: + """In the Heisenberg picture the empty term is the core term ``evolved_operator`` re-adds. + + The fixture Hamiltonian carries no identity term, so one is injected: with a zero core term the + check could not tell the core term apart from the absent-term answer. + """ + problem = load_problem(DATA / "random_exact.msgpack") + with_identity = MajoranaOperator( + {**problem.operator.terms, (): 0.75}, problem.n_modes + ) + prop = MajoranaPropagator( + with_identity, + problem.monomial_circuit.initial_state, + cutoff=2 * problem.n_modes, + comm=serial_comm, + ) + prop.build_graph(problem.monomial_circuit.to_circuit()) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + assert evolved.terms[()] == pytest.approx(0.75) + + coefficients = prop.evolved_operator_coefficients([()], parameters) + + assert coefficients[0] == pytest.approx(evolved.terms[()]) + + +def test_schrodinger_identity_is_a_state_amplitude(serial_comm) -> None: + """In the Schrodinger picture the identity is an ordinary term, not a core term.""" + prop, parameters = _pauli_propagator(serial_comm, schrodinger_cutoff=4) + + evolved = prop.evolved_operator(parameters, atol=0.0) + identity = Pauli("", ()) + assert identity in evolved.terms + + coefficients = prop.evolved_operator_coefficients([identity], parameters) + + assert coefficients[0] == pytest.approx(evolved.terms[identity]) + + +def test_empty_term_list_returns_empty_array(serial_comm) -> None: + """An empty query is an empty answer, not a degenerate probe.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + + coefficients = prop.evolved_operator_coefficients( + [], problem.monomial_circuit.parameters + ) + + assert coefficients.shape == (0,) + + +def test_out_of_range_term_raises(serial_comm) -> None: + """The engine encodes with the checked path, so a term outside the system is rejected.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + + with pytest.raises(RuntimeError, match="out of range"): + prop.evolved_operator_coefficients( + [(2 * problem.n_modes,)], problem.monomial_circuit.parameters + ) + + +def test_base_term_slots_hook_raises(serial_comm) -> None: + """The base hook is a default rather than an abstract method. + + A front-end that leaves it alone still constructs; only the lookup fails. + """ + prop = MajoranaPropagator( + MajoranaOperator({(0,): 1.0}, 2), [], cutoff=2, comm=serial_comm + ) + + with pytest.raises(NotImplementedError, match="_term_slots"): + MonomialPropagator._term_slots(prop, (0,)) + + +@pytest.mark.parametrize( + "term", + [ + pytest.param((14, 0), id="unsorted"), + pytest.param((0, 0), id="repeated"), + pytest.param((-1,), id="negative"), + ], +) +def test_non_canonical_majorana_term_is_rejected(serial_comm, term) -> None: + """A raw tuple that is not a canonical monomial raises instead of answering wrongly. + + The engine keys terms by an order-insensitive bitset, so ``(14, 0)`` would resolve to the row + of ``(0, 14)`` and read back its coefficient *without* the anticommutation sign -- and a lookup + has no coefficient of its own to put that sign on. ``(0, 0)`` would likewise resolve to the + identity row. Both are silent wrong answers, so the front-end validates as + [Majorana][monoprop.majorana.Majorana] does. + """ + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + + with pytest.raises(ValueError, match="Majorana indices must be"): + prop.evolved_operator_coefficients([term], problem.monomial_circuit.parameters) + + +def test_canonicalizing_an_unsorted_majorana_term_recovers_the_sign( + serial_comm, +) -> None: + """``Majorana.from_unsorted`` is the supported way to ask for a non-canonical product.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + canonical = next(term for term in evolved.terms if len(term) == 2) + + term, sign = Majorana.from_unsorted(*reversed(canonical)) + assert term.indices == canonical + assert sign == -1.0 + + coefficient = sign * prop.evolved_operator_coefficients([term], parameters)[0] + assert coefficient == pytest.approx(-evolved.terms[canonical]) + + +def test_majorana_accepts_a_numpy_index_array(serial_comm) -> None: + """An index array is a valid term: the guard is on iterability, and ndarray is not a Sequence.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + parameters = problem.monomial_circuit.parameters + + evolved = prop.evolved_operator(parameters, atol=0.0) + term = next(iter(evolved.terms)) + + coefficients = prop.evolved_operator_coefficients([np.array(term)], parameters) + + assert coefficients[0] == pytest.approx(evolved.terms[term]) + + +def test_majorana_front_end_rejects_a_pauli_term(serial_comm) -> None: + """A term from the other front-end's vocabulary raises TypeError, not an obscure failure.""" + problem = load_problem(DATA / "random_exact.msgpack") + prop = _majorana_propagator(problem, serial_comm) + + with pytest.raises(TypeError, match="Majorana objects or index sequences"): + prop.evolved_operator_coefficients( + [Pauli("X", (0,))], problem.monomial_circuit.parameters + ) + + +def test_pauli_front_end_rejects_a_raw_slot_tuple(serial_comm) -> None: + """Symplectic slots are an engine-internal encoding, so the Pauli front-end takes Pauli only.""" + prop, parameters = _pauli_propagator(serial_comm) + + with pytest.raises(TypeError, match="Pauli objects"): + prop.evolved_operator_coefficients([(0, 1)], parameters) + + +def test_repeated_term_is_answered_once_per_occurrence(serial_comm) -> None: + """The query is a list, not a set: each occurrence gets its own slot with the same value.""" + prop, parameters = _pauli_propagator(serial_comm) + + evolved = prop.evolved_operator(parameters, atol=0.0) + term = next(iter(evolved.terms)) + + coefficients = prop.evolved_operator_coefficients([term, term], parameters) + + assert coefficients.shape == (2,) + assert coefficients[0] == coefficients[1] == pytest.approx(evolved.terms[term]) diff --git a/tests/test_monoprop_mpi.py b/tests/test_monoprop_mpi.py index de4accfe..7a31501a 100644 --- a/tests/test_monoprop_mpi.py +++ b/tests/test_monoprop_mpi.py @@ -134,6 +134,47 @@ def test_expectation_value_and_gradient_functional(self, lih_fermionic_spin_exac combined_grad, lih_fermionic_spin_exact.exact_gradient, atol=1e-9 ) + def test_evolved_operator_coefficients_is_rank_local( + self, lih_fermionic_spin_exact + ): + """Each rank answers for the terms it owns, and the ranks' answers sum to the serial result. + + ``evolved_operator_coefficients`` is rank-local exactly as ``evolved_operator`` is: a term another + rank owns reads back as 0. The identity is excluded from the query -- the core term is + replicated on every rank, so the sum below would count it once per rank. + """ + problem = lih_fermionic_spin_exact + parameters = problem.monomial_circuit.parameters + + serial, serial_circuit = _make_mp(problem, MPI.COMM_SELF) + serial.build_graph(serial_circuit) + # Sorted, so every rank queries the same terms in the same order without relying on the + # enumeration order being reproducible across processes. + terms = sorted( + t for t in serial.evolved_operator(parameters, atol=0.0).terms if t + ) + assert terms + expected = serial.evolved_operator_coefficients(terms, parameters) + + world, world_circuit = _make_mp(problem, MPI.COMM_WORLD) + world.build_graph(world_circuit) + local = world.evolved_operator_coefficients(terms, parameters) + + summed = np.zeros_like(local) + MPI.COMM_WORLD.Allreduce(local, summed, op=MPI.SUM) + + np.testing.assert_allclose(summed, expected, atol=1e-9) + + if MPI.COMM_WORLD.size > 1: + # Pin that the split is real, so the sum above is not just every rank answering in + # full. Counted against the terms each rank *owns* rather than against its nonzero + # coefficients, which would drift if an owned term evolved to exactly 0. Collective, + # so the assertion holds or fails identically on every rank. + local_terms = world.evolved_operator(parameters, atol=0.0).terms + owned = sum(1 for t in local_terms if t) # the identity is on every rank + busiest = MPI.COMM_WORLD.allreduce(owned, op=MPI.MAX) + assert busiest < len(terms) + @pytest.mark.mpi(min_size=2) def test_custom_communicator_split(self, lih_fermionic_spin_exact): rank = MPI.COMM_WORLD.Get_rank() From dee3edd38c1e166f8dfa95f3cb285e8b8dd28f8a Mon Sep 17 00:00:00 2001 From: Josh Cudby Date: Thu, 3 Sep 2026 14:49:07 +0000 Subject: [PATCH 2/2] =?UTF-8?q?docs(propagator):=20=F0=9F=93=9D=20trim=20t?= =?UTF-8?q?he=20evolved=5Foperator=5Fcoefficients=20prose?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The comments and docstrings added with the feature were outliers against their neighbours: the header block was 12 lines where the file's median doc block is 1 and its next-longest are 6 and 7, and `evolved_operator_coefficients` was 36 lines against a module median of 4. Cut to the contracts a reader cannot infer, dropping the worked examples and the reasoning that the docs page and the API reference already carry. The `_term_slots` hooks lose their `Raises:` sections -- they are private, so `gen_api_dump.py` never renders them, and the same raises were documented on the public method too. The evaluation page keeps the overview and defers the contract list to the API reference. No behaviour change. Assisted-by: Claude Opus 5 (1M context) --- cpp/include/monoprop/MonomialPropagator.h | 17 +++----- docs/content/docs/features/evaluation.mdx | 44 ++++---------------- src/monoprop/majorana_propagator.py | 19 +++------ src/monoprop/monomial_propagator.py | 50 ++++++++--------------- src/monoprop/pauli_propagator.py | 7 +--- 5 files changed, 35 insertions(+), 102 deletions(-) diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index 914216ba..4f7711f6 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -283,18 +283,11 @@ class MonomialPropagator { auto evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>>; - /// Coefficients of the requested monomials only, in query order; a term the operator does not - /// carry yields 0. The same contraction and decode as evolved_operator_terms(), but the index is - /// probed with the caller's keys rather than enumerated, so nothing is materialised per surviving - /// term. No atol: the caller named the terms it wants, so magnitude filtering would silently zero - /// some of them. The empty term resolves to core_term() in the Heisenberg picture, where the - /// identity is held apart from the index, and through the probe in the Schrodinger picture, where - /// it is an ordinary index row (the state's identity amplitude) -- so on every key this agrees - /// with what evolved_operator_terms()'s callers publish, the core term they re-add included. - /// Keys must be canonical monomials, as the constructor's initial_operator keys must be: the - /// encode is order-insensitive, so a repeated or unsorted index would resolve to the canonical - /// term's row and lose the sign of the reordering. - /// Contracts non-inplace. Rank-local, as evolved_operator_terms() is. + /// Coefficients of the requested monomials only, in query order; a term the operator does not carry yields 0. + /// As evolved_operator_terms(), but the index is probed with the caller's keys rather than enumerated, and the + /// empty key yields core_term() (Heisenberg) or the indexed identity amplitude (Schrodinger). No atol: the caller + /// named its terms, so filtering by magnitude would silently zero some. Keys must be canonical -- the encode is + /// order-insensitive, so an unsorted or repeated index drops the reordering's sign. Non-inplace. Rank-local. auto evolved_operator_coefficients(const VecD ¶meters, const std::vector &terms) -> std::vector>; diff --git a/docs/content/docs/features/evaluation.mdx b/docs/content/docs/features/evaluation.mdx index 0d3dbd31..ffdb4ec9 100644 --- a/docs/content/docs/features/evaluation.mdx +++ b/docs/content/docs/features/evaluation.mdx @@ -107,12 +107,10 @@ simulator — use [MajoranaPropagator.evolved_operator][majorana_propagator.Majo [MajoranaPropagator.evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] decodes the *whole* evolved operator: it enumerates the operator index, -materialises a key per surviving term, and hands back a term dictionary. When -only a handful of terms are wanted — reading a few amplitudes out of an evolved -state in the Schrödinger picture, say — use -[evolved_operator_coefficients][monoprop.monomial_propagator.MonomialPropagator.evolved_operator_coefficients] -instead. It *probes* the index with the terms you name, so the work scales with -the terms you asked for rather than with the terms the operator carries: +materialises a key per surviving term, and hands back a term dictionary. +Use [evolved_operator_coefficients][monoprop.monomial_propagator.MonomialPropagator.evolved_operator_coefficients] +when only a handful of terms are wanted; reading a few amplitudes out of an evolved +state in the Schrödinger picture, for example. ```python notest # The coefficients of these three terms alone, in this order. @@ -120,39 +118,11 @@ terms = [Majorana(0, 1), Majorana(2, 3), Majorana(0, 3)] coeffs = sim.evolved_operator_coefficients(terms, parameters) ``` -The terms are given in the front-end's own vocabulary — [Majorana][] terms (or -raw index tuples) for [MajoranaPropagator][], [Pauli][] terms for -[PauliPropagator][] — and the returned NumPy array carries one coefficient per +The terms are given in the front-end's own vocabulary: [Majorana][] terms (or +raw index tuples) for [MajoranaPropagator][]; [Pauli][] terms for +[PauliPropagator][]. The returned NumPy array carries one coefficient per requested term, in the order requested. -Three contracts are worth knowing: - -- **A term the evolved operator does not carry reads back as `0`.** There is - deliberately no `atol`: you named the terms you want, so filtering by - magnitude would silently zero some of them. Threshold the returned array - instead. -- **Terms must be canonical, and a non-canonical one raises.** The engine keys - its terms by an order-insensitive bitset, so a raw `(1, 0)` would resolve to - the row of `(0, 1)` — and reordering Majorana indices anticommutes them, so - its coefficient carries a sign that a lookup has nowhere to put. Raw tuples - are therefore validated, not normalized: sort with - [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply - the sign it hands back yourself. The `terms` keys of a - [MajoranaOperator][] are already canonical, as are the qubits of a - [Pauli][]. -- **It is rank-local, exactly as - [evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] - is.** Under MPI each rank answers only for the terms in its own hash - partition, and a term another rank owns reads back as `0`; sum across the - communicator to reassemble the operator. - -The empty (identity) term agrees with -[evolved_operator][majorana_propagator.MajoranaPropagator.evolved_operator] -either way: it is the core term in the Heisenberg picture, where the identity is -held apart from the operator index, and the state's identity amplitude in the -Schrödinger picture, where it is an ordinary term. Repeating a term in one query -simply answers it twice. - This does **not** make the evolution faster. The graph contraction is identical and still the dominant cost; what disappears is the decode and dictionary build over every surviving term. diff --git a/src/monoprop/majorana_propagator.py b/src/monoprop/majorana_propagator.py index 404e5c57..9b20d6ac 100644 --- a/src/monoprop/majorana_propagator.py +++ b/src/monoprop/majorana_propagator.py @@ -124,25 +124,16 @@ def evolved_operator( def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: """Encode a Majorana term into the engine's index tuple. - Majorana indices *are* the engine's keys, so this only unwraps the term. A raw sequence is - round-tripped through [Majorana][monoprop.majorana.Majorana] to reject a non-canonical - product: the engine's encode is order-insensitive, so ``(1, 0)`` would resolve to the row - of ``(0, 1)`` and read back its coefficient without the anticommutation sign, which this - cannot apply to a coefficient it has yet to look up. Canonicalize first with - [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply the sign it - returns yourself; the ``terms`` keys of a - [MajoranaOperator][monoprop.majorana.MajoranaOperator] are already canonical. + Majorana indices *are* the engine's keys, so this only unwraps the term. A raw sequence goes + through [Majorana][monoprop.majorana.Majorana], which rejects a non-canonical product rather + than drop the anticommutation sign a lookup cannot carry; for one of those, use + [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply its sign. Args: - term: A [Majorana][monoprop.majorana.Majorana] term, or a raw index sequence - (a tuple, or an index array). + term: A [Majorana][monoprop.majorana.Majorana] term, or a raw index sequence. Returns: The term's Majorana indices. - - Raises: - TypeError: If the term is neither a Majorana term nor an iterable of indices. - ValueError: If a raw sequence is not sorted, or holds a negative or repeated index. """ if isinstance(term, Majorana): return term.indices diff --git a/src/monoprop/monomial_propagator.py b/src/monoprop/monomial_propagator.py index 9b1677ed..8a0d54fa 100644 --- a/src/monoprop/monomial_propagator.py +++ b/src/monoprop/monomial_propagator.py @@ -542,24 +542,17 @@ def evolved_operator( def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: """Encode one operator term into the raw index tuple the engine keys terms by. - The front-end counterpart to the decode ``evolved_operator`` performs. Not abstract, so - that a front-end which has no term encoding still constructs; it only has to provide this - to support [evolved_operator_coefficients][]. - - Terms must be canonical -- the engine's encode is order-insensitive, so a non-canonical - product would resolve to the canonical term's row and lose the sign of the reordering -- so - an implementation validates rather than normalizes: it has no coefficient to put the sign - on yet. + The front-end counterpart to the decode ``evolved_operator`` performs; a default rather than + an abstract method, so a front-end with no term encoding still constructs and only needs + this for [evolved_operator_coefficients][]. Implementations *validate* canonical terms + rather than normalizing them: the encode is order-insensitive, and a normalizing encode has + no coefficient to put the reordering's sign on. Args: term: A single term in this front-end's own vocabulary. Returns: - The term's engine index tuple (Majorana indices, or symplectic slots in the Pauli - basis). - - Raises: - NotImplementedError: If this front-end does not encode terms. + The term's engine index tuple: Majorana indices, or symplectic slots in the Pauli basis. """ raise NotImplementedError( f"{type(self).__name__} does not implement _term_slots, so it cannot look up " @@ -574,27 +567,16 @@ def evolved_operator_coefficients( """Return the coefficients of ``terms`` alone in the evolved operator, in the order given. A cheaper - [evolved_operator][monoprop.monomial_propagator.MonomialPropagator.evolved_operator] - for when only a few terms are wanted. The graph contraction is identical -- and still the - dominant cost -- but the operator index is *probed* with these terms rather than - enumerated, so the decode and the term dictionary cost one entry per term *requested* - instead of one per term the evolved operator carries. Reading a handful of amplitudes out - of an evolved state (Schrodinger picture) is the motivating case. - - A term the evolved operator does not carry reads back as ``0``. There is deliberately no - ``atol``: the caller named the terms it wants, so magnitude filtering would silently zero - some of them -- threshold the returned array instead. - - Terms must be canonical, and a non-canonical one is rejected rather than normalized: the - anticommutation sign of the reordering belongs on the coefficient, which is what this is - about to look up. Repeating a term simply answers it twice. - - The empty (identity) term agrees with ``evolved_operator`` too: it is the core term in the - Heisenberg picture, where the identity is held apart from the operator index, and the - state's identity amplitude in the Schrodinger picture, where it is an ordinary term. - - Rank-local, as ``evolved_operator`` is: under MPI each rank answers for the terms it owns, - and a term another rank owns reads back as ``0``. + [evolved_operator][monoprop.monomial_propagator.MonomialPropagator.evolved_operator] when + only a few terms are wanted: the contraction is identical and still dominant, but the index + is *probed* with these terms rather than enumerated, so the decode costs one entry per term + *requested* rather than one per term the evolved operator carries. + + A term the operator does not carry reads back as ``0``, and there is deliberately no + ``atol`` -- the caller named its terms, so magnitude filtering would silently zero some of + them. Terms must be canonical; for a Majorana product that is not, use + [Majorana.from_unsorted][monoprop.majorana.Majorana.from_unsorted] and apply the sign it + returns. A repeated term is answered once per occurrence. Args: terms: The terms to look up, in this front-end's own vocabulary. diff --git a/src/monoprop/pauli_propagator.py b/src/monoprop/pauli_propagator.py index be345229..ab3e3c0f 100644 --- a/src/monoprop/pauli_propagator.py +++ b/src/monoprop/pauli_propagator.py @@ -120,17 +120,14 @@ def evolved_operator( def _term_slots(self, term: OperatorTerm) -> tuple[int, ...]: """Encode a qubit Pauli term into the engine's symplectic slots. - Only [Pauli][monoprop.pauli.Pauli] terms are accepted: the slots are an engine-internal - encoding, so a raw slot sequence is rejected rather than passed through. + Slots are an engine-internal encoding, so only [Pauli][monoprop.pauli.Pauli] terms are + accepted; a raw slot sequence is rejected rather than passed through. Args: term: A [Pauli][monoprop.pauli.Pauli] term. Returns: The term's symplectic slot indices. - - Raises: - TypeError: If the term is not a [Pauli][monoprop.pauli.Pauli]. """ if not isinstance(term, Pauli): raise TypeError(