diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index b684ff58..4f7711f6 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -283,6 +283,14 @@ 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. + /// 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>; + 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 e97555b6..054c9d5e 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -1169,4 +1169,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..ffdb4ec9 100644 --- a/docs/content/docs/features/evaluation.mdx +++ b/docs/content/docs/features/evaluation.mdx @@ -102,3 +102,27 @@ 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. +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. +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][]. The returned NumPy array carries one coefficient per +requested term, in the order requested. + +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 e9c4d77b..10185fcc 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 @@ -126,6 +125,30 @@ 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 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. + + Returns: + The term's Majorana indices. + """ + 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..8a0d54fa 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,65 @@ 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; 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. + """ + 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] 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. + 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 dd23a1de..f026f555 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]): @@ -120,6 +120,24 @@ 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. + + 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. + """ + 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()