diff --git a/AGENTS.md b/AGENTS.md index cf2f2fd0..8b5bb824 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -96,11 +96,43 @@ Key files: - **`Monomial`** (`cpp/monoprop/core/Monomial.h`) = `Bitset<2*N>`: ONE basis operator, two bits per mode/qubit. Basis-agnostic — read as a Majorana product, or as a Pauli string (JW image). Collections: `MonomialList` (no coeffs) and `MonomialMap` (monomial → real coeff). -- **Row access** (`cpp/monoprop/detail/operator/RowAccess.h`): the one backend-agnostic vocabulary - (`materialize_row`, `assign_row`, `row_popcount`, `for_each_row_position`) over the dense - `MonomialList` and the packed `detail::OperatorIndex`. Any template parameterized on the row - store must include that header — the `OperatorIndex` overloads live in `monoprop::`, so ADL cannot - find them from a `monoprop::detail` argument. +- **The row-store seam** (`cpp/monoprop/detail/operator/RowAccess.h`): a dense monomial is a transient, + not the storage. Four accessors — `materialize_row`, `assign_row`, `row_popcount`, + `for_each_row_position` — and three backends answer them: `MonomialList`, `detail::OperatorIndex` + (packed position lists) and `detail::SparseRowStore` (fixed-width mode lanes plus one 2-bit-per-slot + `codes` word per row). Reach rows through the accessors, never through a backend's own API, and add any + fourth backend to `cpp/tests/row_accessor_tests.cpp`, which asserts that all of them agree through every + accessor. Any template parameterized on the row store must include that header — the overloads live in + `monoprop::`, so ADL cannot find them from a `monoprop::detail` argument. + `algebra/CodesAlgebra.h` is the structural algebra on a sparse row, one function per dense counterpart, + reading the `codes` word instead of looping over storage words, plus `sparse_toggle` — the product + `M ⊕ G` as one merge over two ascending lane arrays. It is exact, not an approximation: + `cpp/tests/codes_algebra_tests.cpp` and `codes_product_tests.cpp` assert agreement with each dense + version over the fixtures and randomized rows. Change one side and you must change the other. A product + can occupy more modes than either input; past its scratch capacity `sparse_toggle` reports `overflowed` + and the caller must fall back to the dense product — never truncate, because a truncated mode list still + carries a plausible-looking `codes` word. +- **Which backend, and where it is bound**: a propagator uses one of the two row stores, chosen once from + its mode count by `SparseRowStore::preferred_for_modes()` — a build-time constant + (`monoprop_SPARSE_ROW_MIN_MODES`, derived in `CMakeLists.txt` from whether `ARCH_FLAG` is actually + emitted rather than from the option that asks for it, and deliberately not a cache entry) because what + moves the crossover is the target ISA. `monoprop_ROW_STORE=dense|sparse` forces it process-wide; an + unrecognized value throws rather than falling back, since the point of setting it is to know which + backend ran. `MPOperator` holds one pointer per backend with exactly one non-null and binds the live one + via `with_store` — **once per layer, inside `build_layer`**, never per term: the scan asks the store for + a row per anticommuting term, so everything downstream is templated on the store + (`LayerBuildEngine`, `fused_find_and_collect`, `probe_incoming_queries`). Off that path, + use `MPOperator`'s forwarding accessors; there is no accessor handing out a store, because there is no + one type to hand out — which is why `MonomialPropagator` exposes `for_each_term()`/`num_local_terms()` + rather than the `indexing()` it used to. Every C++ case runs a second time under + `monoprop_ROW_STORE=sparse` (the `sparse-rows` ctest label) — every fixture is below the crossover, so + without that the sparse backend would ship untested; `cpp/tests/row_store_selection_tests.cpp` is what + fails if the variable stops reaching the propagator. The two backends agree on term sets and values but + not on term *order*, so compare them with `just diff-baseline-sparse` (tolerance), never + `just diff-baseline` (byte-wise). A benchmark run records the backend it resolved to in its artifact's + `meta` (`monoprop_row_store` as asked, `row_store_effective` as run) and `REPORT.md` shows both: under + the default `auto` the setting alone does not identify the backend, and the two differ in footprint and + in accumulation order. - **`Basis` / the `Algebra` policy** (`cpp/monoprop/algebra/`): the two algebras are sibling models (`MajoranaAlgebra`, `PauliAlgebra` in `algebra/Algebra.h`) over shared structural primitives (`algebra/AlgebraCommon.h`). The propagation backbone (the scan/fold in `detail/evolution/...`) is diff --git a/CMakeLists.txt b/CMakeLists.txt index c98b529e..980deaa6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -77,6 +77,20 @@ endif() include(${PROJECT_SOURCE_DIR}/cmake/compiler_flags/Sanitizers.cmake) include(${PROJECT_SOURCE_DIR}/cmake/compiler_flags/CXXFlags.cmake) +# The storage width at or above which a propagator picks the support-form row store over the dense one. +# Derived from whether ARCH_FLAG is actually emitted rather than from the option that asks for it, and +# deliberately not a cache entry: what moves the crossover is the target ISA, so a stale cached value +# would silently pick the wrong backend after a flag change. +# +# Thresholds are the first full 32-mode block where sparse is clearly faster than dense beyond +# run-to-run noise. Expect about +/-1 block variation across machines due to cache and popcount +# throughput. +if(ARCH_FLAG) + set(monoprop_SPARSE_ROW_MIN_MODES 768) +else() + set(monoprop_SPARSE_ROW_MIN_MODES 256) +endif() + # report on compiler flags in use message(STATUS "Configuring a ${CMAKE_BUILD_TYPE} build") string(TOUPPER ${CMAKE_BUILD_TYPE} _cmake_build_type_upper) @@ -92,6 +106,10 @@ message( " Build-type-specific : ${_cmake_build_type_specific_flags}" ) message(STATUS " Vectorization flag : ${ARCH_FLAG}") +message( + STATUS + " Sparse rows from : ${monoprop_SPARSE_ROW_MIN_MODES} modes" +) message( STATUS " Project defaults : ${CMAKE_CXX${CMAKE_CXX_STANDARD}_STANDARD_COMPILE_OPTION} ${monoprop_CXX_FLAGS}" diff --git a/README.md b/README.md index 3f03bbbe..dfff6b5c 100644 --- a/README.md +++ b/README.md @@ -117,8 +117,12 @@ uv sync --all-groups --all-extras -v # installs the workspace, incl. the benc uv run python -m pytest -m "not mpi" # Python tests (serial) just test-mpi # Python + C++ tests under MPI just test-wide # Python + C++ unit tests with a 64-bit TermIndex +just test-sparse-rows # Python tests with the support-form row backend forced ``` +The C++ suite runs against both row backends: `ctest` registers every case a second +time with `monoprop_ROW_STORE=sparse`, labelled `sparse-rows`. + See the [testing guide](https://docs.monoprop.algorithmiq.tech/testing) for the with/without-MPI details and the rank matrix. diff --git a/benches/conftest.py b/benches/conftest.py index e6284a69..58e09fd5 100644 --- a/benches/conftest.py +++ b/benches/conftest.py @@ -214,6 +214,7 @@ def _meta(nodes: int, ranks_per_node: int) -> dict[str, Any]: "nodes": nodes, "ranks_per_node": ranks_per_node, "monoprop_threads": os.environ.get("monoprop_NUM_THREADS", "default"), # noqa: SIM112 + "monoprop_row_store": os.environ.get("monoprop_ROW_STORE") or "auto", # noqa: SIM112 "cpu_count_logical": psutil.cpu_count(logical=True), "cpu_count_physical": psutil.cpu_count(logical=False), "hostname": socket.gethostname(), @@ -238,6 +239,27 @@ def _meta(nodes: int, ranks_per_node: int) -> dict[str, Any]: return meta +def _record_row_store(propagator: Any) -> None: + """Fold one propagator's resolved row backend into this run's metadata. + + ``monoprop_ROW_STORE`` says what was asked for, not what ran: unset lets the mode width pick, and + the crossover it picks against is a build-time constant. The two backends accumulate a term sum in + different orders and have different footprints, so a report has to name the one that ran. Widths + differ within a run, hence so can the backend: a disagreement records as ``"mixed"`` rather than + letting the last propagator speak for the others. + """ + if _rank() != 0: + return + # Read straight off the binding, with no getattr fallback: a benchmark whose whole job is to name + # the backend that ran must fail loudly against an extension that cannot say, not quietly record + # nothing. + resolved = "sparse" if propagator._simulator.rows_are_sparse else "dense" + seen = _RESULTS["meta"].get("row_store_effective") + _RESULTS["meta"]["row_store_effective"] = ( + resolved if seen in (None, resolved) else "mixed" + ) + + def _params(config: pytest.Config) -> dict[str, Any]: """Return the resolved random-problem hyperparameters (defaults included).""" return { @@ -353,6 +375,7 @@ def _record_model_stats( ) -> None: """Record term count, operator memory breakdown and footprint under ``key``.""" _record("opsize", key, {"terms": _reduce_sum(comm, propagator.size())}) + _record_row_store(propagator) # Placement is only observable while the propagator's threads are alive. _record_placement(comm) @@ -553,6 +576,7 @@ def built_graph( # Under MPI the operator is partitioned, so sum the partitions. _record("opsize", picture, {"terms": _reduce_sum(bench_comm, mp.size())}) + _record_row_store(mp) # Settled RSS once the build's transients are released -- the persistent # footprint the per-operation peak cannot see. diff --git a/cpp/include/monoprop/MonomialPropagator.h b/cpp/include/monoprop/MonomialPropagator.h index b684ff58..deaec980 100644 --- a/cpp/include/monoprop/MonomialPropagator.h +++ b/cpp/include/monoprop/MonomialPropagator.h @@ -156,14 +156,30 @@ class MonomialPropagator { /// graph_layers(), optimizer order) or a per-gate one (length n_gates()); on a tie, per-layer wins. auto set_parameter_mapping(const VecZ ¶meter_mapping) -> void; - /// This rank's monomial → coefficient index. Single-partition only — see require_single_partition_. - auto indexing() -> detail::OperatorIndex & { - require_single_partition_("indexing()"); - return *mp_op_.store; + /// This rank's terms as fn(monomial, coefficient index), in the index's own slot order. + /// Single-partition only — see require_single_partition_. No accessor for the store itself: which + /// backend holds the rows is a runtime choice (see MPOperator::with_store), so there is no one type + /// to hand out. + template + auto for_each_term(Fn &&fn) const -> void { + require_single_partition_("for_each_term()"); + mp_op_.for_each_term(std::forward(fn)); } - auto indexing() const -> const detail::OperatorIndex & { - require_single_partition_("indexing()"); - return *mp_op_.store; + /// This rank's term count. Single-partition only. + auto num_local_terms() const -> size_t { + require_single_partition_("num_local_terms()"); + return mp_op_.size(); + } + + /// Whether this propagator's rows live in the support-form backend. Which one it is is decided once + /// at construction (see use_sparse_rows_); this reports the answer rather than re-deriving it. + /// Partition-transparent: every partition of a facade makes the same choice from the same width and + /// the same environment, so partition 0 speaks for all of them. + auto rows_are_sparse() const -> bool { + if (is_partition_facade()) { + return first_partition_().rows_are_sparse(); + } + return mp_op_.rows_are_sparse(); } /// Per-layer (cos_inds, local_cycles, cross_rank_sin_send, cross_rank_sin_recv) for this @@ -272,7 +288,7 @@ class MonomialPropagator { /// Contract the graph into the operator (Heisenberg) or state (Schrodinger). `inplace` consumes the /// graph and updates internal state; otherwise nothing is mutated. Core term excluded either way. - /// Coefficients are positioned by the owning partition's indexing(), so on a facade the result is + /// Coefficients are positioned by the owning partition's row index, so on a facade the result is /// the per-partition blocks concatenated in partition order: the same multiset as an unpartitioned /// run, but not positionally stable across partition counts — and the count is auto-picked from the /// host's core count unless pinned. Use evolved_operator_terms() when positions must mean something. @@ -312,8 +328,13 @@ class MonomialPropagator { detail::MatchedEpochSet matched_scratch_; // A perf hint, never a correctness constraint: overflow spills losslessly. Sized to the cutoff's - // structural position bound when it has one. - auto packed_inline_width_() const -> size_t; + // structural position bound when it has one. Shared by both backends -- the bound is in physical + // slots, which is what an OperatorIndex inline width and a SparseRowStore slot count both count. + auto row_width_bound_() const -> size_t; + + // Which row backend to build on, decided once per propagator. See config::Settings::row_store for + // why an unrecognized monoprop_ROW_STORE is a throw rather than a silent fall back to auto. + auto use_sparse_rows_() const -> bool; // `requested` 0 ⇒ env/auto. Returns 1 for the ordinary single-partition path. static auto resolve_partition_count_(size_t requested, mpi::Comm comm) -> size_t; diff --git a/cpp/monoprop/Bitset.h b/cpp/monoprop/Bitset.h index 6ad135a3..4891ab7f 100644 --- a/cpp/monoprop/Bitset.h +++ b/cpp/monoprop/Bitset.h @@ -216,6 +216,17 @@ class Bitset { return os; } }; + +// The splitmix64 finalizer. Every hash in the engine ends here, and the value routes MPI ownership, so +// this must stay bit-identical wherever it is reached from. +[[nodiscard]] constexpr auto splitmix_finalize(uint64_t x) noexcept -> uint64_t { + x ^= x >> 30; + x *= 0xbf58476d1ce4e5b9ULL; + x ^= x >> 27; + x *= 0x94d049bb133111ebULL; + x ^= x >> 31; + return x; +} } // namespace monoprop template @@ -223,14 +234,7 @@ struct SplitmixHash; template struct SplitmixHash> { - static constexpr auto mix(uint64_t x) noexcept -> uint64_t { - x ^= x >> 30; - x *= 0xbf58476d1ce4e5b9ULL; - x ^= x >> 27; - x *= 0x94d049bb133111ebULL; - x ^= x >> 31; - return x; - } + static constexpr auto mix(uint64_t x) noexcept -> uint64_t { return monoprop::splitmix_finalize(x); } auto operator()(const monoprop::Bitset &bs) const noexcept -> size_t { constexpr size_t W = monoprop::Bitset::num_words(); diff --git a/cpp/monoprop/CMakeLists.txt b/cpp/monoprop/CMakeLists.txt index 1b233211..c6ca6d19 100644 --- a/cpp/monoprop/CMakeLists.txt +++ b/cpp/monoprop/CMakeLists.txt @@ -22,6 +22,7 @@ target_compile_definitions( PUBLIC $<$:monoprop_ENABLE_MPI> $<$:monoprop_WIDE_TERM_INDEX> + monoprop_SPARSE_ROW_MIN_MODES=${monoprop_SPARSE_ROW_MIN_MODES} ) # flags to prepend @@ -100,6 +101,7 @@ target_compile_definitions( INTERFACE $<$:monoprop_ENABLE_MPI> $<$:monoprop_WIDE_TERM_INDEX> + monoprop_SPARSE_ROW_MIN_MODES=${monoprop_SPARSE_ROW_MIN_MODES} ) target_compile_features(monoprop INTERFACE cxx_std_23) diff --git a/cpp/monoprop/algebra/CMakeLists.txt b/cpp/monoprop/algebra/CMakeLists.txt index ddc05bd5..be0c7083 100644 --- a/cpp/monoprop/algebra/CMakeLists.txt +++ b/cpp/monoprop/algebra/CMakeLists.txt @@ -6,6 +6,7 @@ target_sources( FILES "Algebra.h" "AlgebraCommon.h" + "CodesAlgebra.h" "MajoranaAlgebra.h" "PauliAlgebra.h" ) diff --git a/cpp/monoprop/algebra/CodesAlgebra.h b/cpp/monoprop/algebra/CodesAlgebra.h new file mode 100644 index 00000000..a2a06a67 --- /dev/null +++ b/cpp/monoprop/algebra/CodesAlgebra.h @@ -0,0 +1,317 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +// The structural algebra on a sparse row's `codes` word, one function per dense counterpart in +// AlgebraCommon.h / PauliAlgebra.h / MajoranaAlgebra.h. Each is exact, not approximate: agreement with +// the dense version over the tests/data fixtures and randomized rows is asserted in +// cpp/tests/codes_algebra_tests.cpp, and that test is the gate on ever making these the default. +// +// Why any of this is possible in one word: a mode's two physical positions 2m, 2m+1 become the 2-bit +// field of slot j, so quantities the dense form derives from a per-word masked shift chain over the +// whole register become popcounts of two masks of a single word, independent of the storage width. +// With n = popcount(row_occupied_bits(codes)) and d = popcount(row_paired_bits(codes)): +// +// or_sum (support/Pauli weight) = n popcount_sum (length) = n + d xor_sum = n - d +// is_paired <=> d == n pair_swap = swap the two bits of every field +// Y letters = fields equal to 0b01 +// +// Nothing here reads a mode lane except codes_interleave_phase, which is inherently a two-row +// operation, and codes_cutoff_sums when a logical width narrower than the storage width makes some +// modes inactive. +// +// Names are prefixed rather than overloading the dense functions: while both representations are live +// a call site should say which one it means, and overload resolution between `MonomialLike auto` and +// SparseRow would decide that silently. + +#include +#include +#include +#include + +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/algebra/PauliAlgebra.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +namespace monoprop::detail { + +// Set bits of `codes` in slots strictly below `slot`. The dense counterpart is the prefix popcount an +// interleave scan maintains word by word. +[[nodiscard]] inline auto codes_popcount_below(RowCodes codes, size_t slot) noexcept -> size_t { + // A shift by 2*kRowMaxSlots would be undefined, and "below every slot" is the whole word anyway. + if (slot >= kRowMaxSlots) { + return static_cast(std::popcount(codes)); + } + return static_cast(std::popcount(codes & ((RowCodes{1} << (2 * slot)) - 1))); +} + +// or_sum / popcount_sum / xor_sum for a row whose every slot is inside the active window. Two popcounts +// and no reference to the storage width, where the dense cutoff_sums runs a masked shift chain per word +// and needs CutoffMasks to avoid rederiving the masks per term. +[[nodiscard]] inline auto codes_cutoff_sums(RowCodes codes) noexcept -> CutoffSums { + const auto n = static_cast(std::popcount(row_occupied_bits(codes))); + const auto d = static_cast(std::popcount(row_paired_bits(codes))); + return {n - d, n + d, n}; +} + +// The same, restricted to the active window. inactive_mode_prefix is the count of leading *physical* +// modes the logical width excludes -- storage_num_modes - logical_num_modes, half of +// CutoffMasks::active_bit_offset -- and the dense form applies it as `mono >> active_bit_offset`. +// +// The inactive modes are exactly the low ones, so they are a prefix of the ascending slots and drop out +// with one shift. A propagator's rows never carry them (a term is built from logical indices, which map +// into the window), so the common case is the zero-prefix early exit; the general path exists because +// the dense function it must agree with accepts such a monomial. +[[nodiscard]] inline auto codes_cutoff_sums(const SparseRow &row, size_t inactive_mode_prefix) noexcept -> CutoffSums { + if (inactive_mode_prefix == 0) { + return codes_cutoff_sums(row.codes); + } + const size_t n = row.num_slots(); + size_t inactive_slots = 0; + while (inactive_slots < n && row.mode(inactive_slots) < inactive_mode_prefix) { + ++inactive_slots; + } + if (inactive_slots >= kRowMaxSlots) { + return {0, 0, 0}; + } + return codes_cutoff_sums(row.codes >> (2 * inactive_slots)); +} + +// Both cutoffs keep a fully paired row unconditionally, exactly as the dense ones do: those are the +// only terms contributing to an expectation value against a product reference state. +[[nodiscard]] inline auto codes_length_cutoff(const SparseRow &row, + unsigned int cutoff, + size_t inactive_mode_prefix) noexcept -> bool { + const auto sums = codes_cutoff_sums(row, inactive_mode_prefix); + return sums.xor_sum == 0 || sums.popcount_sum <= cutoff; +} + +[[nodiscard]] inline auto codes_support_cutoff(const SparseRow &row, + unsigned int cutoff, + size_t inactive_mode_prefix) noexcept -> bool { + const auto sums = codes_cutoff_sums(row, inactive_mode_prefix); + return sums.xor_sum == 0 || sums.or_sum <= cutoff; +} + +// The counterparts of CutoffEvaluator::passes_with_popcount, one per concrete cutoff functor. Same +// shortcut and same reasoning: the predicate is `xor_sum == 0 || measure <= cutoff`, so a popcount +// already at or below the bound proves keep without reading the row at all (or_sum <= popcount_sum makes +// that sound for the support cutoff too). `popcount_sum` is the whole-register count, which can only +// exceed the active-window one, so the shortcut stays conservative when a logical width is narrower than +// the storage width. +// +// There is no evaluator argument: which cutoff a propagator has is fixed for its lifetime, so the caller +// resolves it once per gate rather than re-branching per term. +[[nodiscard]] inline auto codes_length_passes_with_popcount(const SparseRow &row, + unsigned int cutoff, + size_t popcount_sum, + size_t inactive_mode_prefix) noexcept -> bool { + return popcount_sum <= cutoff || codes_length_cutoff(row, cutoff, inactive_mode_prefix); +} + +[[nodiscard]] inline auto codes_support_passes_with_popcount(const SparseRow &row, + unsigned int cutoff, + size_t popcount_sum, + size_t inactive_mode_prefix) noexcept -> bool { + return popcount_sum <= cutoff || codes_support_cutoff(row, cutoff, inactive_mode_prefix); +} + +// Every occupied mode holds both of its positions, i.e. every field is 0b11. Unoccupied modes are not +// slots at all and are trivially paired, which is why this needs no window argument -- and matches the +// dense is_paired, which likewise checks the whole register. +[[nodiscard]] inline auto codes_is_paired(RowCodes codes) noexcept -> bool { + return row_paired_bits(codes) == row_occupied_bits(codes); +} + +// The pair-swap involution J: swap the two physical bits of every mode (u <-> v). Occupancy is +// preserved -- 0b01 <-> 0b10 and 0b11 is fixed -- so the mode lanes are untouched and the row's whole +// transform is this one word operation, against a per-word masked shift-and-or on the dense side. +[[nodiscard]] constexpr auto codes_pair_swap(RowCodes codes) noexcept -> RowCodes { + return ((codes & kRowLoBits) << 1) | ((codes >> 1) & kRowLoBits); +} + +// Y letters: v=1, u=0 under the JW image, so the field is exactly 0b01. +[[nodiscard]] inline auto codes_pauli_y_count(RowCodes codes) noexcept -> size_t { + const RowCodes v = codes & kRowLoBits; + const RowCodes u = (codes >> 1) & kRowLoBits; + return static_cast(std::popcount(v & ~u)); +} + +// Whether two Pauli strings anticommute: the symplectic inner product x_P.z_G + z_P.x_G mod 2, which +// dense-side is p.parity_and(pair_swap(g)). Sparse-side the AND is over the modes both rows occupy, so +// it is a merge of the two ascending lane arrays. +[[nodiscard]] inline auto codes_pauli_anticommutes(const SparseRow &p, const SparseRow &g) noexcept -> bool { + const size_t np = p.num_slots(); + const size_t ng = g.num_slots(); + unsigned int parity = 0; + size_t i = 0; + for (size_t k = 0; k < ng; ++k) { + const size_t g_mode = g.mode(k); + while (i < np && p.mode(i) < g_mode) { + ++i; + } + if (i < np && p.mode(i) == g_mode) { + // popcount of the pair-swapped generator field ANDed with p's, both 2 bits wide. + const unsigned int swapped = ((g.code(k) & 1U) << 1) | ((g.code(k) >> 1) & 1U); + parity ^= static_cast(std::popcount(p.code(i) & swapped)) & 1U; + } + } + return parity != 0; +} + +// Ordering sign (-1)^S of maj.gen, S = #{set bits of maj strictly below each set bit of gen} mod 2 over +// physical bit positions. Slots ascend in the mode and a mode's low position is 2*mode, so ascending +// slots are ascending positions and one merge walk over the two rows suffices -- O(slots), where the +// dense form is either a prefix-XOR scan over every word or a per-layer full-width mask W plus a +// parity_and per term. The mask has no sparse counterpart worth building: W is dense by construction +// (roughly half the register), so this replaces it with the direct walk instead. +[[nodiscard]] inline auto codes_interleave_phase(const SparseRow &maj, const SparseRow &gen) noexcept -> int { + const size_t nm = maj.num_slots(); + const size_t ng = gen.num_slots(); + unsigned int parity = 0; + size_t below_slots = 0; // maj slots at modes strictly below the current generator mode + for (size_t k = 0; k < ng; ++k) { + const size_t g_mode = gen.mode(k); + // Monotone across k, since generator modes ascend: the whole walk is one pass over each row. + while (below_slots < nm && maj.mode(below_slots) < g_mode) { + ++below_slots; + } + const size_t below = codes_popcount_below(maj.codes, below_slots); + const unsigned int g_code = gen.code(k); + // maj's bits at g_mode itself, if it occupies it: position 2*g_mode is below 2*g_mode+1 and so + // counts for the generator's high bit only. + const unsigned int m_code = (below_slots < nm && maj.mode(below_slots) == g_mode) ? maj.code(below_slots) : 0U; + if ((g_code & 1U) != 0U) { + parity ^= static_cast(below) & 1U; + } + if ((g_code & 2U) != 0U) { + parity ^= static_cast(below + (m_code & 1U)) & 1U; + } + } + return parity == 0 ? 1 : -1; +} + +// Pauli's per-term rotation sign, the counterpart of pauli_rotation_sign. Same exponent, one merge walk +// instead of a masked pass over the generator's nonzero words: +// e = g_y + sum(y_mono - y_new) + 2 * sum(v_mono & x_gen), sign = (e mod 4 == 1 ? -1 : +1) +// where per mode the code's low bit is v (physical position 2*mode) and its high bit is u, so a Y letter +// is the field 0b01 and x = u ^ v is one bit. +// +// The dense version restricts its sums to the words the generator occupies, on the grounds that +// elsewhere the two Y counts cancel and x_gen is zero. Per *mode* that argument is exact and tighter: a +// mode the generator misses has new_mono's field equal to mono's, so the Y terms cancel, and x_gen = 0 +// kills the cross term. So this walks the generator's slots and reads mono's field at each, which also +// means new_mono never has to exist -- the sign comes out of the same merge the toggle does. +[[nodiscard]] inline auto codes_pauli_rotation_sign(const SparseRow &mono, const SparseRow &gen) noexcept -> int { + auto delta = static_cast(codes_pauli_y_count(gen.codes)); + long cross = 0; + const size_t nm = mono.num_slots(); + const size_t ng = gen.num_slots(); + size_t i = 0; + for (size_t k = 0; k < ng; ++k) { + const size_t g_mode = gen.mode(k); + while (i < nm && mono.mode(i) < g_mode) { + ++i; + } + const unsigned int a = (i < nm && mono.mode(i) == g_mode) ? mono.code(i) : 0U; + const unsigned int b = gen.code(k); + delta += (a == 0b01U) ? 1 : 0; + delta -= ((a ^ b) == 0b01U) ? 1 : 0; + cross += ((a & 1U) != 0U && ((b ^ (b >> 1)) & 1U) != 0U) ? 1 : 0; + } + return mod4(delta + (2 * cross)) == 1 ? -1 : 1; +} + +// The product row of a term and a generator, and the overlap the emit phase needs. `codes` and +// `num_slots` describe the row written into `out_lanes`. +struct SparseProduct { + RowCodes codes = 0; + size_t num_slots = 0; + size_t overlap = 0; // popcount(mono & gen); 0 and meaningless when overflowed + bool overflowed = false; +}; + +// mono (+) gen: per mode the fields XOR, a mode whose field cancels to zero disappears, and the overlap +// is the popcount of the fields' AND. This is the dense fused_xor_into in support form, and it is the +// operation the whole representation exists for -- one merge over two ascending lane arrays, O(slots), +// against a pass over every storage word. +// +// out_lanes.size() lanes are available and must not exceed kRowMaxSlots (a codes word's worth). The +// product can occupy more modes than either input: up to mono's slots plus the generator's, so a +// scratch row needs CutoffEvaluator::max_mode_bound() + the generator's locality, not just the bound. +// When even that is not enough the result is reported as overflowed rather than truncated -- a truncated +// mode list keeps a plausible-looking codes word, which is exactly how the Stage 3 bench measured a +// capacity bug as if it were a speedup. On overflow the caller must fall back to the dense product; +// `overlap` is partial and is deliberately not returned. +[[nodiscard]] inline auto sparse_toggle(const SparseRow &mono, + const SparseRow &gen, + std::span out_lanes) noexcept -> SparseProduct { + assert(out_lanes.size() <= kRowMaxSlots && "sparse_toggle capacity exceeds one codes word"); + const size_t nm = mono.num_slots(); + const size_t ng = gen.num_slots(); + SparseProduct result; + size_t used = 0; + bool over = false; + + const auto emit = [&](size_t mode, unsigned int code) { + if (used == out_lanes.size()) { + over = true; + return; + } + out_lanes[used] = static_cast(mode); + result.codes |= static_cast(code) << (2 * used); + ++used; + }; + + size_t i = 0; + size_t j = 0; + while (!over && i < nm && j < ng) { + const size_t m_mode = mono.mode(i); + const size_t g_mode = gen.mode(j); + if (m_mode < g_mode) { + emit(m_mode, mono.code(i)); + ++i; + } + else if (g_mode < m_mode) { + emit(g_mode, gen.code(j)); + ++j; + } + else { + const unsigned int a = mono.code(i); + const unsigned int b = gen.code(j); + result.overlap += static_cast(std::popcount(a & b)); + if (const unsigned int c = a ^ b; c != 0U) { + emit(m_mode, c); + } + ++i; + ++j; + } + } + while (!over && i < nm) { + emit(mono.mode(i), mono.code(i)); + ++i; + } + while (!over && j < ng) { + emit(gen.mode(j), gen.code(j)); + ++j; + } + if (over) { + return SparseProduct{0, 0, 0, true}; + } + result.num_slots = used; + return result; +} + +} // namespace monoprop::detail diff --git a/cpp/monoprop/detail/EnvConfig.h b/cpp/monoprop/detail/EnvConfig.h index 30120417..ea96fa79 100644 --- a/cpp/monoprop/detail/EnvConfig.h +++ b/cpp/monoprop/detail/EnvConfig.h @@ -15,7 +15,9 @@ #pragma once #include +#include #include +#include #include // Single home for runtime environment configuration. Kept dependency-free by design, because it is @@ -23,9 +25,15 @@ // // monoprop_NUM_THREADS positive int (1..1e6), else ignored → num_threads // monoprop_PARTITIONS int N | "auto" | "off"; parsed where it is used (resolve_partition_count_) +// monoprop_ROW_STORE "auto" (default) | "dense" | "sparse"; unset == auto → row_store namespace monoprop::config { +// Which row backend a propagator builds on. Auto is the measured crossover +// (SparseRowStore::preferred_for_modes); the two explicit values force one backend for every +// propagator in the process, which is how the suite is run either way -- see row_store below. +enum class RowStore : std::uint8_t { Auto, Dense, Sparse }; + namespace detail { inline auto parse_positive_int(const char *text) -> std::optional { @@ -43,17 +51,40 @@ inline auto parse_positive_int(const char *text) -> std::optional { return static_cast(value); } +inline auto parse_row_store(const char *text) -> std::optional { + using enum RowStore; + if (text == nullptr || text[0] == '\0' || std::strcmp(text, "auto") == 0) { + return Auto; + } + if (std::strcmp(text, "dense") == 0) { + return Dense; + } + if (std::strcmp(text, "sparse") == 0) { + return Sparse; + } + return std::nullopt; +} + } // namespace detail struct Settings { std::optional num_threads; + // nullopt means monoprop_ROW_STORE held something unrecognized -- reported rather than ignored, + // unlike every other setting here: this one exists to prove the sparse backend was exercised, so a + // typo that silently fell back to auto would mean believing a configuration ran that never did. The + // throw is raised by the propagator, which has the exception types; this header stays + // dependency-free. + std::optional row_store = RowStore::Auto; }; -// Parse the environment once; the Settings are cached and shared across TUs. +// Parse the environment once; the Settings are cached and shared across TUs. Cached deliberately: a +// setting must not change between two propagators in one process, since the row backend is part of a +// monomial's hash and so of every cross-propagator comparison. inline auto get() -> const Settings & { static const Settings settings = [] { Settings s; s.num_threads = detail::parse_positive_int(std::getenv("monoprop_NUM_THREADS")); + s.row_store = detail::parse_row_store(std::getenv("monoprop_ROW_STORE")); return s; }(); return settings; diff --git a/cpp/monoprop/detail/evolution/layer_build/Engine.h b/cpp/monoprop/detail/evolution/layer_build/Engine.h index d073f00c..2816cc6f 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Engine.h +++ b/cpp/monoprop/detail/evolution/layer_build/Engine.h @@ -45,7 +45,7 @@ namespace monoprop::detail { template inline auto append_inserted_endpoints(CosMask &cos_all, size_t combined_size, const MPOperator &op) -> void { const size_t cos_lo = combined_size; - const size_t cos_hi = op.store->size(); + const size_t cos_hi = op.size(); CosineWordBuilder end_b; for (size_t idx = cos_lo; idx < cos_hi; ++idx) { end_b.push_index(idx); @@ -295,7 +295,7 @@ struct ContractSink { }; // Owns build_layer's machinery over a compile-time Sink policy. combined_size = the pre-layer operator size. -template +template struct LayerBuildEngine { struct DeferredSelfMiss { Monomial mono; @@ -304,6 +304,9 @@ struct LayerBuildEngine { double v_src = 0.0; // ContractSink only: op_pre[src] captured at scan emit; 0 for GraphSink }; MPOperator &local_op; // scanned, looked up, and grown by the inserts + // The row backend, bound once per layer by build_layer (MPOperator::with_store) and held by + // reference here so nothing on the per-term path re-derives it. + Store &store; mpi::Comm comm; size_t R; size_t my_rank; @@ -321,6 +324,7 @@ struct LayerBuildEngine { Sink sink; LayerBuildEngine(MPOperator &local_op_, + Store &store_, mpi::Comm comm_, size_t R_, size_t my_rank_, @@ -328,6 +332,7 @@ struct LayerBuildEngine { size_t combined_size_, Sink &&sink_) : local_op(local_op_), + store(store_), comm(comm_), R(R_), my_rank(my_rank_), @@ -380,7 +385,8 @@ struct LayerBuildEngine { std::vector &send = sink.send_buffer(queries_r, src_val_r, combined_qv_); std::vector> inc_q; mpi::begin_alltoallv(send, comm).wait_into(inc_q); - auto resp = resolve_incoming(inc_q, local_op, R, is_leader_pass, matched, combined_size, sink); + auto resp = + resolve_incoming(inc_q, local_op, store, R, is_leader_pass, matched, combined_size, sink); std::vector resp_recv = response_recv_counts(); std::vector> inc_r; mpi::begin_alltoallv(resp, comm, /*skip_self=*/false, &resp_recv).wait_into(inc_r); @@ -437,9 +443,9 @@ struct LayerBuildEngine { } auto key_at = [&](size_t k) -> const Monomial & { return deferred_self_misses[k].mono; }; sink.prepare_deferred(n_miss); - insert_absent_terms(local_op, n_miss, key_at, [&](size_t k, size_t base) { + insert_absent_terms(local_op, store, n_miss, key_at, [&](size_t k, size_t base) { const auto &m = deferred_self_misses[k]; - assign_row(*local_op.store, base + k, m.mono); + assign_row(store, base + k, m.mono); sink.emit_deferred(k, base + k, m.src, m.phase, m.v_src); }); } @@ -469,7 +475,7 @@ struct LayerBuildEngine { size_t lo, size_t hi, bool is_leader_pass) -> void { - const size_t op_size = local_op.store->size(); + const size_t op_size = store.size(); std::array, kResolveBatch> keys; std::array phases; std::array srcs; @@ -493,7 +499,7 @@ struct LayerBuildEngine { if (m == 0) { break; } - local_op.store->find_batch(keys.data(), m, found.data()); + store.find_batch(keys.data(), m, found.data()); for (size_t j = 0; j < m; ++j) { double v_src = 0.0; if constexpr (Sink::wants_values) { @@ -560,79 +566,88 @@ auto build_layer(MPOperator &local_op, } assert(fused_scale_coeffs == nullptr || (local_coeffs && &local_coeffs->get() == fused_scale_coeffs)); - FusedScanResult fused = [&] { - double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; - return with_algebra(basis, [&]() { - return fused_find_and_collect(local_op, - gen, - cut_eval, - cut_st, - coeffs, - only_rotate_len_k, - R, - my_rank, - /*capture_values=*/use_fused, - sweep_ptr, - cos_build); - }); - }(); - - CosMask cos_all; - if (fused.cos_blocks.size() == 1) { - // The serial scan produces a single cosine block set — take it wholesale. - cos_all = std::move(fused.cos_blocks[0]); - } - else { - // Cosine block sets are disjoint and ascending; concatenate in order. - for (const auto &block : fused.cos_blocks) { - cos_all.total_count += block.total_count; - cos_all.blocks.insert(cos_all.blocks.end(), block.blocks.begin(), block.blocks.end()); + // The single place the row backend is bound: everything from the scan to the inserts is templated on + // it (the key batch, the row writes, the resolve), and binding it here means the choice costs one + // branch per layer instead of one per term. Both arms are instantiated, so this doubles the + // scan/engine template instantiations -- the same trade with_algebra already makes for Basis. + std::shared_ptr storage = local_op.with_store([&](S &store) -> std::shared_ptr { + FusedScanResult fused = [&] { + double *const sweep_ptr = fused_scale ? fused_scale_coeffs->data() : nullptr; + return with_algebra(basis, [&]() { + return fused_find_and_collect(local_op, + store, + gen, + cut_eval, + cut_st, + coeffs, + only_rotate_len_k, + R, + my_rank, + /*capture_values=*/use_fused, + sweep_ptr, + cos_build); + }); + }(); + + CosMask cos_all; + if (fused.cos_blocks.size() == 1) { + // The serial scan produces a single cosine block set — take it wholesale. + cos_all = std::move(fused.cos_blocks[0]); } - } - fused.cos_blocks = std::vector{}; - - auto run = [&](Sink sink) -> std::shared_ptr { - LayerBuildEngine eng(local_op, - comm, - R, - my_rank, - matched_scratch, - /*combined_size=*/local_op.store->size(), - std::move(sink)); - eng.run_exchange(/*is_leader_pass=*/true, - std::move(fused.leader_queries), - std::move(fused.leader_src), - std::move(fused.leader_val)); - eng.run_exchange(/*is_leader_pass=*/false, - std::move(fused.follower_queries), - std::move(fused.follower_src), - std::move(fused.follower_val)); - - return eng.finish(std::move(cos_all), out_cos); - }; - - std::shared_ptr storage; - if (use_fused) { - const double inv_cos = fused_scale ? 1.0 / cos_build : 1.0; // pre-cos recovery factor for hit v_tgt - storage = run(ContractSink{.R = R, - .my_rank = my_rank, - .fc = *fused_contract, - .op_coeffs = coeffs, - .fused_scale = fused_scale, - .inv_cos = inv_cos, - .schrodinger = schrodinger, - .basis = basis}); - } - else { - storage = run(GraphSink{R, my_rank}); - } + else { + // Cosine block sets are disjoint and ascending; concatenate in order. + for (const auto &block : fused.cos_blocks) { + cos_all.total_count += block.total_count; + cos_all.blocks.insert(cos_all.blocks.end(), block.blocks.begin(), block.blocks.end()); + } + } + fused.cos_blocks = std::vector{}; + + auto run = [&](Sink sink) -> std::shared_ptr { + LayerBuildEngine eng(local_op, + store, + comm, + R, + my_rank, + matched_scratch, + /*combined_size=*/store.size(), + std::move(sink)); + eng.run_exchange(/*is_leader_pass=*/true, + std::move(fused.leader_queries), + std::move(fused.leader_src), + std::move(fused.leader_val)); + eng.run_exchange(/*is_leader_pass=*/false, + std::move(fused.follower_queries), + std::move(fused.follower_src), + std::move(fused.follower_val)); + + return eng.finish(std::move(cos_all), out_cos); + }; + + std::shared_ptr layer; + if (use_fused) { + const double inv_cos = fused_scale ? 1.0 / cos_build : 1.0; // pre-cos recovery factor for hit v_tgt + layer = run(ContractSink{.R = R, + .my_rank = my_rank, + .fc = *fused_contract, + .op_coeffs = coeffs, + .fused_scale = fused_scale, + .inv_cos = inv_cos, + .schrodinger = schrodinger, + .basis = basis}); + } + else { + layer = run(GraphSink{R, my_rank}); + } + return layer; + }); // Recompute metadata rides with the layer so it survives every graph transform. scaled_count is the // post-insert operator size: the fold truncated to it reproduces the "all anticommuting" cos // bit-for-bit with no stored bitmap. Fused mode has no LayerCore to stamp. if (storage != nullptr) { storage->generator_words.assign(gen.data(), gen.data() + mpi_detail::kWords); - storage->scaled_count = static_cast(local_op.store->size()); + storage->scaled_count = static_cast(local_op.size()); } return storage; diff --git a/cpp/monoprop/detail/evolution/layer_build/Resolve.h b/cpp/monoprop/detail/evolution/layer_build/Resolve.h index ef90d91f..9bede256 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Resolve.h +++ b/cpp/monoprop/detail/evolution/layer_build/Resolve.h @@ -45,9 +45,12 @@ struct IncomingProbe { // Phases 1-2, read-only w.r.t. operator contents. QW = per-record stride: the plain query width, or // kQueryWordsFused for the fused resolver. The caller runs Phase 3, then insert_incoming_misses. -template > +// `store` is passed alongside `op` rather than taken off it: the caller is inside build_layer, which has +// already bound the concrete backend, and re-entering with_store() here would bind it a second time. +template > auto probe_incoming_queries(const std::vector &incoming, // serialized, one VecZ per sender MPOperator &op, + Store &store, size_t rank_count) -> IncomingProbe { constexpr size_t W = QW; IncomingProbe pr; @@ -83,8 +86,8 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on pr.phase_of[g] = ph; } { - const size_t op_size = op.store->size(); - op.store->find_batch(pr.mono.data(), pr.nq_total, pr.idx_of.data()); + const size_t op_size = store.size(); + store.find_batch(pr.mono.data(), pr.nq_total, pr.idx_of.data()); for (size_t g = 0; g < pr.nq_total; ++g) { if (pr.idx_of[g] >= op_size) { // kNotFound is size_t max → also lands here pr.idx_of[g] = kMissingIndex; @@ -93,7 +96,7 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on } // Phase 2 ((sender,query) prefix order): each miss takes the next index base+j. - pr.base = op.store->size(); + pr.base = store.size(); for (size_t g = 0; g < pr.nq_total; ++g) { if (pr.idx_of[g] == kMissingIndex) { pr.idx_of[g] = pr.base + pr.miss_g.size(); @@ -105,17 +108,18 @@ auto probe_incoming_queries(const std::vector &incoming, // serialized, on // Phase 4 (bulk insert of the distinct absent terms) into op slots [base, base+n_miss). Call after the // caller's Phase-3 scatter, which reads pre-insert op_coeffs for hits and needs base == op.size(). -template -auto insert_incoming_misses(MPOperator &op, const IncomingProbe &pr) -> void { +template +auto insert_incoming_misses(MPOperator &op, Store &store, const IncomingProbe &pr) -> void { const size_t n_miss = pr.miss_g.size(); if (n_miss == 0) { return; } insert_absent_terms( op, + store, n_miss, [&](size_t j) -> const Monomial & { return pr.mono[pr.miss_g[j]]; }, - [&](size_t j, size_t base) { assign_row(*op.store, base + j, pr.mono[pr.miss_g[j]]); }); + [&](size_t j, size_t base) { assign_row(store, base + j, pr.mono[pr.miss_g[j]]); }); } // resolve_incoming / process_responses are the picture-independent cross-rank exchange skeletons; what @@ -125,16 +129,18 @@ auto insert_incoming_misses(MPOperator &op, const IncomingProbe +template auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ per sender MPOperator &op, + Store &store, size_t rank_count, bool is_leader_pass, MatchedEpochSet &matched, size_t combined_size, // pre-layer op size: bounds the matched set Sink &sink) -> std::vector> { using Resp = typename Sink::Response; - const IncomingProbe pr = probe_incoming_queries(incoming, op, rank_count); + const IncomingProbe pr = + probe_incoming_queries(incoming, op, store, rank_count); std::vector> responses(rank_count); for (size_t s = 0; s < rank_count; ++s) { responses[s].assign(pr.goff[s + 1] - pr.goff[s], Sink::init_response()); @@ -156,7 +162,7 @@ auto resolve_incoming(const std::vector &incoming, // serialized, one VecZ } } - insert_incoming_misses(op, pr); + insert_incoming_misses(op, store, pr); return responses; } diff --git a/cpp/monoprop/detail/evolution/layer_build/Scan.h b/cpp/monoprop/detail/evolution/layer_build/Scan.h index d5b8a77e..2070ecf8 100644 --- a/cpp/monoprop/detail/evolution/layer_build/Scan.h +++ b/cpp/monoprop/detail/evolution/layer_build/Scan.h @@ -33,6 +33,7 @@ #include "monoprop/detail/mpi/MPIUtils.h" #include "monoprop/detail/operator/InvertedIndex.h" #include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/detail/operator/RowAccess.h" namespace monoprop::detail { @@ -159,15 +160,15 @@ inline auto rotation_dynamic_gate(std::optional only_rotate_len_k, // phase_factor is the basis-specific sign only: Majorana interleave_phase, still to be folded with // hermitian_phase at emit; Pauli pauli_rotation_sign, already rotation-ready. -template -[[gnu::always_inline]] inline auto emit_term_products(const OperatorIndex &ham, +template +[[gnu::always_inline]] inline auto emit_term_products(const Store &ham, size_t i, const typename A::GenContext &ctx, Monomial &new_mono, size_t &overlap, int &phase_factor) -> void { Monomial mono; - ham.for_each_position(i, [&](size_t pos) { mono.set(pos); }); + for_each_row_position(ham, i, [&](size_t pos) { mono.set(pos); }); const Monomial &gen = A::generator(ctx); new_mono = mono ^ gen; overlap = mono.count_and(gen); @@ -191,8 +192,12 @@ struct FusedScanResult { // deterministic. `fused_scale_coeffs` (no length cap only; must alias coeffs.data()) scales every anticommuting // coeff in place by `fused_scale_cos`=cos(2·build_angle), so no cosine set is built and a hit's stored // value is post-cos (resolve recovers it via 1/cos). -template +// Templated on the row backend rather than reading it off `op`: build_layer binds the concrete store +// once per layer (MPOperator::with_store) and hands it down, so a scan that touches one row per +// anticommuting term pays no branch for the choice. +template auto fused_find_and_collect(const MPOperator &op, + const Store &store, const Monomial &gen, const CutoffEvaluator &cutoff_eval, const CutoffContext &cut_st, @@ -236,7 +241,7 @@ auto fused_find_and_collect(const MPOperator &op, return res; } const uint64_t *const row_parity_ptr = g_odd ? inverted_index.row_parity_words() : nullptr; - const size_t n = op.store->size(); + const size_t n = store.size(); // The fused sweep writes fused_scale_coeffs[i] for every anticommuting i < n, so it must be the // very array the reads come from and cover the full operator — a violation corrupts 1/cos recovery. assert(fused_scale_coeffs == nullptr || (fused_scale_coeffs == coeffs.data() && coeffs.size() >= n)); @@ -277,7 +282,7 @@ auto fused_find_and_collect(const MPOperator &op, Monomial new_mono; size_t overlap = 0; int phase_factor = 0; - emit_term_products(*op.store, i, ectx, new_mono, overlap, phase_factor); + emit_term_products(store, i, ectx, new_mono, overlap, phase_factor); // Structural cutoff on the partner M⊕G, unless upper_atol rescues it (CutoffContext::is_above_upper). const size_t new_pop = mono_pop + gen_pop - 2 * overlap; const bool struct_pass = cutoff_eval.passes_with_popcount(new_mono, new_pop); @@ -353,7 +358,7 @@ auto fused_find_and_collect(const MPOperator &op, if (cut_st.is_below_sin(abs_c)) { continue; } - const size_t mono_pop = op.store->popcount(i); + const size_t mono_pop = row_popcount(store, i); const bool is_follower = (w.foll >> tz) & 1u; emit(mono_pop, i, abs_c, v_src, is_follower); } @@ -369,7 +374,7 @@ auto fused_find_and_collect(const MPOperator &op, if (cut_st.is_below_sin(abs_c)) { continue; } - const size_t mono_pop = op.store->popcount(i); + const size_t mono_pop = row_popcount(store, i); const bool is_follower = (w.foll >> tz) & 1u; emit(mono_pop, i, abs_c, v_src, is_follower); } @@ -380,7 +385,7 @@ auto fused_find_and_collect(const MPOperator &op, for (uint64_t m = w.overlap; m; m &= m - 1) { const size_t tz = static_cast(std::countr_zero(m)); const size_t i = w.base + tz; - const size_t mono_pop = op.store->popcount(i); + const size_t mono_pop = row_popcount(store, i); if (mono_pop > static_cast(*only_rotate_len_k)) { continue; } diff --git a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl index 34814906..2d795e1f 100644 --- a/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl +++ b/cpp/monoprop/detail/monomial_propagator/MonomialPropagator.inl @@ -166,12 +166,20 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope auto op = schrodinger_ ? generate_paired_op(sc / 2 + sc % 2, logical_num_modes_) : local_heisenberg_terms; const size_t expected_local_terms = std::max(1, op.size() / std::max(1, num_ranks)); - // Must run before the store: packed_inline_width_() derives the packed-row width from cutoff_fn_. + // Must run before the store: row_width_bound_() derives the row width from cutoff_fn_. regenerate_cutoff_fn_(); - mp_op_.store = std::make_unique>(packed_inline_width_()); - mp_op_.store->reserve(expected_local_terms); - // Store replaced: drop the stale lazy inverted index so it rebuilds against the new store. - mp_op_.inverted_index_.reset(); + // Replaces the store MPOperator's default: same width, but now with the cutoff-derived row width, + // which is only knowable after regenerate_cutoff_fn_() above. set_store() drops the stale lazy + // inverted index with it. Nothing has been inserted yet, so there are no rows to migrate. + const size_t bound = row_width_bound_(); + if (use_sparse_rows_()) { + mp_op_.set_store(std::make_unique>( + detail::SparseRowStore::slots_for_bound(bound))); + } + else { + mp_op_.set_store(std::make_unique>(bound)); + } + mp_op_.reserve_terms(expected_local_terms); size_t i = 0; // The initial monomials are distinct, so emplace (insert-if-absent) is an assigning insert here. @@ -179,7 +187,7 @@ MonomialPropagator::MonomialPropagator(const OperatorDict &initial_ope const auto &mono = materialize_row(op, r); if (my_rank == find_rank(mono, num_ranks)) { mp_op_.append_term(mono); - mp_op_.store->emplace(mono, i++); + mp_op_.index_term(mono, i++); } } @@ -343,18 +351,34 @@ auto MonomialPropagator::partitioned_graph_memory_usage_() const -> Gr } template -auto MonomialPropagator::packed_inline_width_() const -> size_t { - constexpr size_t kMax = detail::OperatorIndex::kMaxInlinePositions; +auto MonomialPropagator::row_width_bound_() const -> size_t { + // Schrodinger's initial fill is the whole paired basis, whose rows are far longer than the bound + // derived from cutoff_fn_ below -- both backends would spill most of the initial rows. constexpr size_t kDefault = detail::OperatorIndex::kDefaultInlinePositions; if (schrodinger_) { return kDefault; } // The bound is already in physical slots (CutoffEvaluator::max_slot_bound), so nothing to scale. - const auto bound = detail::CutoffEvaluator(cutoff_fn_).max_slot_bound(); - if (!bound) { - return kDefault; + return detail::CutoffEvaluator(cutoff_fn_).max_slot_bound().value_or(kDefault); +} + +template +auto MonomialPropagator::use_sparse_rows_() const -> bool { + const auto &settings = config::get(); + if (!settings.row_store) { + throw PropagatorConfigError( + R"(monoprop_ROW_STORE must be "auto", "dense" or "sparse". Unset it to pick by system size.)"); + } + using enum monoprop::config::RowStore; + switch (*settings.row_store) { + case Dense: + return false; + case Sparse: + return true; + case Auto: + break; } - return std::min(*bound, kMax); + return detail::SparseRowStore::preferred_for_modes(NumModes); } template @@ -1107,11 +1131,11 @@ template auto MonomialPropagator::evolved_operator_terms(const VecD ¶meters, double atol) -> std::vector>> { using Term = std::pair>; - // `p` is always unpartitioned here (a partition, or *this), so indexing() is available. + // `p` is always unpartitioned here (a partition, or *this), so for_each_term() is available. const auto collect = [&](MonomialPropagator &p) -> std::vector { std::vector terms; const VecD evolved = p.contract_partially(parameters, false); - p.indexing().for_each([&](const auto &mono, size_t idx) { + p.for_each_term([&](const auto &mono, size_t idx) { if (idx >= evolved.size()) { return; } diff --git a/cpp/monoprop/detail/operator/CMakeLists.txt b/cpp/monoprop/detail/operator/CMakeLists.txt index 939009ff..66003ed3 100644 --- a/cpp/monoprop/detail/operator/CMakeLists.txt +++ b/cpp/monoprop/detail/operator/CMakeLists.txt @@ -9,4 +9,5 @@ target_sources( "OperatorIndex.h" "RowAccess.h" "RowHashTable.h" + "SparseRowStore.h" ) diff --git a/cpp/monoprop/detail/operator/MPOperator.h b/cpp/monoprop/detail/operator/MPOperator.h index 774baf6c..b7a45159 100644 --- a/cpp/monoprop/detail/operator/MPOperator.h +++ b/cpp/monoprop/detail/operator/MPOperator.h @@ -31,6 +31,7 @@ #include "monoprop/core/Monomial.h" #include "monoprop/detail/operator/InvertedIndex.h" #include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" // Forward-declared to break an include cycle with algebra/Algebra.h. namespace monoprop { @@ -61,9 +62,20 @@ class OperatorTermNotFound : public std::runtime_error { template struct MPOperator { - // The store is non-copyable/non-movable, so it is heap-owned by unique_ptr (keeping MPOperator - // itself cheaply movable). Always non-null. - std::unique_ptr> store{std::make_unique>()}; + // The row store is one of two backends, chosen per propagator from its mode count + // (SparseRowStore::preferred_for_modes) and then fixed for the propagator's lifetime. Exactly one + // of these is non-null. + // + // Two pointers rather than one, and not a virtual interface, because the scan asks the store for a + // row per anticommuting term: a branch or an indirect call on that path is not affordable. + // with_store() binds the concrete type once per layer instead -- the same shape as with_algebra() + // for a runtime Basis, and the reason build_layer is a template. Everything off that path goes + // through the forwarding accessors below, which pay one well-predicted branch. + // + // Heap-owned because neither store is copyable or movable (single-writer, and their views borrow + // their arrays), which keeps MPOperator itself cheaply movable. + std::unique_ptr> dense_rows{std::make_unique>()}; + std::unique_ptr> sparse_rows = nullptr; VecD op_coeffs; // Only fully-paired terms score nonzero (see score_new_state_rows_), which on production models is // ~0.07% of the rows -- a dense vector here is 99.9% zeros. state_rows_ is strictly ascending: rows are @@ -85,7 +97,8 @@ struct MPOperator { MPOperator &operator=(MPOperator &&) noexcept = default; MPOperator(const MPOperator &other) - : store(other.store->clone()), + : dense_rows(other.dense_rows ? other.dense_rows->clone() : nullptr), + sparse_rows(other.sparse_rows ? other.sparse_rows->clone() : nullptr), op_coeffs(other.op_coeffs), state_rows_(other.state_rows_), state_vals_(other.state_vals_), @@ -96,23 +109,76 @@ struct MPOperator { basis(other.basis), inverted_index_(other.inverted_index_) {} - auto size() const -> size_t { return store->size(); } + // Binds the live store to a concrete type for the duration of the call. Both arms are instantiated, + // so `f` must be a generic lambda and must return the same type from each. + template + [[gnu::always_inline]] auto with_store(F &&f) -> decltype(auto) { + if (sparse_rows) { + return f(*sparse_rows); + } + return f(*dense_rows); + } + template + [[gnu::always_inline]] auto with_store(F &&f) const -> decltype(auto) { + if (sparse_rows) { + return f(*sparse_rows); + } + return f(*dense_rows); + } + + // Installs a backend, dropping the lazy inverted index with it: the index addresses the old rows, + // and leaving it would let a stale one answer for the new store until its row count happened to + // disagree. One overload per backend rather than a tag, so a call site names the choice. + auto set_store(std::unique_ptr> rows) -> void { + dense_rows = std::move(rows); + sparse_rows.reset(); + inverted_index_.reset(); + } + auto set_store(std::unique_ptr> rows) -> void { + sparse_rows = std::move(rows); + dense_rows.reset(); + inverted_index_.reset(); + } + [[nodiscard]] auto rows_are_sparse() const -> bool { return sparse_rows != nullptr; } + + auto size() const -> size_t { + return with_store([](const auto &rows) { return rows.size(); }); + } // Does not keep the lazy inverted index in sync: appends happen during setup, before the index is // first materialized, so a later append just makes inverted_index() rebuild via its staleness guard. - auto append_term(const Monomial &mono) -> void { store->push_back(mono); } + auto append_term(const Monomial &mono) -> void { + with_store([&mono](auto &rows) { rows.push_back(mono); }); + } + + // Setup-path forwards, kept here rather than exposing a store, so nothing outside has to know which + // backend is live. + auto reserve_terms(size_t n) -> void { + with_store([n](auto &rows) { rows.reserve(n); }); + } + auto index_term(const Monomial &mono, size_t row) -> void { + with_store([&mono, row](auto &rows) { rows.emplace(mono, row); }); + } + [[nodiscard]] auto find(const Monomial &mono) const -> std::optional { + return with_store([&mono](const auto &rows) { return rows.find(mono); }); + } + // This rank's terms as fn(monomial, row), in the index's slot order. Materializes each row. + template + auto for_each_term(Fn &&fn) const -> void { + with_store([&fn](const auto &rows) { rows.for_each(std::forward(fn)); }); + } - // Resync the inverted index after a bulk growth of `store`, preserving has_value() ⟹ rows()==store.size(). + // Resync the inverted index after a bulk growth of the store, preserving has_value() ⟹ rows()==size(). auto reindex_after_growth(size_t base, size_t n) -> void { if (inverted_index_.has_value()) { - inverted_index_->append_rows(*store, base, n); + with_store([this, base, n](const auto &rows) { inverted_index_->append_rows(rows, base, n); }); } } auto inverted_index() const -> const InvertedIndex & { - if (!inverted_index_.has_value() || inverted_index_->rows() != store->size()) { + if (!inverted_index_.has_value() || inverted_index_->rows() != size()) { inverted_index_.emplace(); - inverted_index_->rebuild(*store); + with_store([this](const auto &rows) { inverted_index_->rebuild(rows); }); } return *inverted_index_; } @@ -131,7 +197,7 @@ struct MPOperator { const auto before = init_op_map.size(); erase_if(init_op_map, [this](const auto &kv) { - const auto found = store->find(kv.first); + const auto found = find(kv.first); if (found) { op_coeffs[*found] = kv.second; } @@ -197,7 +263,7 @@ struct MPOperator { for (const auto &[k, v] : op_dict) { // Unchecked by design: the only caller bounds-checks against its logical_num_modes_. const auto mono = indices_to_bitset(k); - const auto rank_evolved_op = store->find(mono); + const auto rank_evolved_op = find(mono); const auto rank_init_op = init_op_map.find(mono); const auto coeff = algebra_encode_coeff(basis, v, mono); @@ -237,14 +303,16 @@ struct MPOperator { VecZ new_inds(size() - state_scored_rows_); std::iota(new_inds.begin(), new_inds.end(), state_scored_rows_); // NOLINT(modernize-use-ranges) - const auto paired_inds = is_fully_paired(new_inds, *store); - state_rows_.reserve(state_rows_.size() + paired_inds.size()); - state_vals_.reserve(state_vals_.size() + paired_inds.size()); - - // The algebra picks the diagonal ⟨b|·|b⟩ phase of each fully-paired term. - algebra_score_state(basis, paired_inds, initial_state, *store, [this](size_t row, double phase) { - state_rows_.push_back(static_cast(row)); - state_vals_.push_back(phase); + with_store([this, &new_inds](const auto &rows) { + const auto paired_inds = is_fully_paired(new_inds, rows); + state_rows_.reserve(state_rows_.size() + paired_inds.size()); + state_vals_.reserve(state_vals_.size() + paired_inds.size()); + + // The algebra picks the diagonal ⟨b|·|b⟩ phase of each fully-paired term. + algebra_score_state(basis, paired_inds, initial_state, rows, [this](size_t row, double phase) { + state_rows_.push_back(static_cast(row)); + state_vals_.push_back(phase); + }); }); state_scored_rows_ = size(); @@ -263,13 +331,14 @@ struct MPOperator { // Callers must pass pairwise-distinct, currently-absent keys: bulk_insert then skips duplicate probes and // slot k deterministically lands at base+k. Call after any pass that reads pre-insert op state // (op.size() must equal the returned base). -template -inline auto insert_absent_terms(MPOperator &op, size_t n, KeyAt &&key_at, PerSlot &&per_slot) -> size_t { - const size_t base = op.store->grow_rows_geometric(n); +template +inline auto insert_absent_terms(MPOperator &op, Store &store, size_t n, KeyAt &&key_at, PerSlot &&per_slot) + -> size_t { + const size_t base = store.grow_rows_geometric(n); for (size_t k = 0; k < n; ++k) { per_slot(k, base); } - op.store->bulk_insert(n, base, std::forward(key_at)); + store.bulk_insert(n, base, std::forward(key_at)); op.reindex_after_growth(base, n); return base; } @@ -329,13 +398,16 @@ struct MPOperatorMemoryBreakdown final { template inline auto estimate_memory_usage(const MPOperator &op) -> MPOperatorMemoryBreakdown { MPOperatorMemoryBreakdown breakdown; - breakdown.operator_terms_bytes = op.store->memory_bytes(); + op.with_store([&breakdown](const auto &rows) { + breakdown.operator_terms_bytes = rows.memory_bytes(); + breakdown.indexing_bytes = rows.index_estimated_memory_bytes(); + breakdown.operator_terms_slack_bytes = rows.slack_bytes(); + }); breakdown.op_coeffs_bytes = op.op_coeffs.capacity() * sizeof(double); // Every representation of the state at once: the sparse scored set plus the dense vector. breakdown.state_coeffs_bytes = op.state_coeffs.capacity() * sizeof(double) + op.state_rows_.capacity() * sizeof(TermIndex) + op.state_vals_.capacity() * sizeof(double); - breakdown.indexing_bytes = op.store->index_estimated_memory_bytes(); breakdown.init_operator_bytes = unordered_flat_map_storage_bytes(op.init_op_map); breakdown.init_operator_entries = op.init_op_map.size(); breakdown.initial_state_bytes = op.initial_state.capacity() * sizeof(size_t); @@ -346,7 +418,6 @@ inline auto estimate_memory_usage(const MPOperator &op) -> MPOperatorM breakdown.inverted_index_sparse_bytes = tiers[1]; breakdown.inverted_index_dense_columns = tiers[2]; } - breakdown.operator_terms_slack_bytes = op.store->slack_bytes(); // State phases are unit-magnitude, so at rest the scored count IS the nonzero count; a live vector needs a scan. breakdown.state_coeffs_nonzero = op.state_coeffs.empty() diff --git a/cpp/monoprop/detail/operator/OperatorIndex.h b/cpp/monoprop/detail/operator/OperatorIndex.h index 397b22d4..f8cd9096 100644 --- a/cpp/monoprop/detail/operator/OperatorIndex.h +++ b/cpp/monoprop/detail/operator/OperatorIndex.h @@ -89,8 +89,7 @@ class OperatorIndex { auto grow_rows_geometric(size_t n) -> size_t { const size_t base = size_; if (capacity() < base + n) { - const size_t cap = capacity(); - reserve_rows(std::max(base + n, cap + (cap / 2) + 1)); + reserve_rows(geometric_row_capacity(base, n, capacity())); } // Default-init grow, not a zeroing resize: every freshly grown row is overwritten by set() // before any read, so a tail zero-fill would be wasted bandwidth. @@ -155,9 +154,7 @@ class OperatorIndex { return overflow_.at(i).count(); } [[nodiscard]] auto memory_bytes() const -> size_t { - size_t total = rows_.capacity() * sizeof(PosT); - total += overflow_.size() * (sizeof(value_type) + sizeof(size_t) + 24); - return total; + return (rows_.capacity() * sizeof(PosT)) + spilled_rows_bytes(overflow_); } auto find(const key_type &key) const -> std::optional { diff --git a/cpp/monoprop/detail/operator/RowAccess.h b/cpp/monoprop/detail/operator/RowAccess.h index 3e7ed0cf..5d249975 100644 --- a/cpp/monoprop/detail/operator/RowAccess.h +++ b/cpp/monoprop/detail/operator/RowAccess.h @@ -14,10 +14,11 @@ #pragma once -// One row-reader/writer vocabulary over both operator backends: the dense MonomialList and the packed -// detail::OperatorIndex. Templates parameterized on the row store (`Rows`) call these unqualified, so -// every such template must include this header — ADL cannot reach monoprop:: from an argument in -// monoprop::detail, and a later declaration is not found for an already-parsed template definition. +// One row-reader/writer vocabulary over all three operator backends: the dense MonomialList, the packed +// detail::OperatorIndex and the fixed-width-lane detail::SparseRowStore. Templates parameterized on the +// row store (`Rows`) call these unqualified, so every such template must include this header — ADL +// cannot reach monoprop:: from an argument in monoprop::detail, and a later declaration is not found for +// an already-parsed template definition. #include #include @@ -25,11 +26,12 @@ #include "monoprop/core/Monomial.h" #include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" namespace monoprop { -// materialize_row() returns a const ref (dense backend, zero-copy) or a fresh value (packed backend), -// so callers must bind with `const auto&` to extend the temporary's lifetime. +// materialize_row() returns a const ref (dense backend, zero-copy) or a fresh value (packed/sparse +// backend), so callers must bind with `const auto&` to extend the temporary's lifetime. template [[nodiscard]] inline auto materialize_row(const std::vector> &op, size_t i) -> const Monomial & { @@ -54,6 +56,9 @@ inline auto for_each_row_position(const std::vector> &op, siz } } +// The packed and support-form backends answer the same four accessors, one overload set each. They are +// written out per backend rather than behind a concept because every call site names its width +// explicitly (`materialize_row(op, i)`), which a concept-constrained `Op` cannot deduce. template [[nodiscard]] inline auto materialize_row(const detail::OperatorIndex &op, size_t i) -> Monomial { return op.row(i); @@ -71,4 +76,34 @@ inline auto for_each_row_position(const detail::OperatorIndex &op, siz op.for_each_position(i, std::forward(fn)); } +template +[[nodiscard]] inline auto materialize_row(const detail::SparseRowStore &op, size_t i) -> Monomial { + return op.row(i); +} +template +inline auto assign_row(detail::SparseRowStore &op, size_t i, const Monomial &mono) -> void { + op.set(i, mono); +} +// A row written from a key that is already in the store's own form -- what the insert of an absent term +// does once the query record it came from was read in that form. Only the support-form store has a form +// of its own, so these two have no counterpart on the other backends: there is no such thing as an +// OperatorIndex-shaped key that is not simply a monomial. +template +inline auto assign_row(detail::SparseRowStore &op, size_t i, const detail::SparseRow &row) -> void { + op.set(i, row); +} +template +inline auto assign_row(detail::SparseRowStore &op, size_t i, const detail::SparseRowKey<2 * NumModes> &key) + -> void { + op.set(i, key); +} +template +[[nodiscard]] inline auto row_popcount(const detail::SparseRowStore &op, size_t i) -> size_t { + return op.popcount(i); +} +template +inline auto for_each_row_position(const detail::SparseRowStore &op, size_t i, Fn &&fn) -> void { + op.for_each_position(i, std::forward(fn)); +} + } // namespace monoprop diff --git a/cpp/monoprop/detail/operator/RowHashTable.h b/cpp/monoprop/detail/operator/RowHashTable.h index acb86332..80273968 100644 --- a/cpp/monoprop/detail/operator/RowHashTable.h +++ b/cpp/monoprop/detail/operator/RowHashTable.h @@ -14,6 +14,7 @@ #pragma once +#include #include #include #include @@ -27,6 +28,23 @@ namespace monoprop::detail { +// The next row-array capacity for a geometric (1.5x) grow by `n` rows from `base` (the pre-growth size), +// given the current capacity. Never exact-fit: an exact fit would realloc the whole operator every layer. +// Shared by OperatorIndex and SparseRowStore, whose grow_rows_geometric() differ only in which arrays +// that capacity gets applied to. +[[nodiscard]] inline auto geometric_row_capacity(size_t base, size_t n, size_t capacity) noexcept -> size_t { + return std::max(base + n, capacity + (capacity / 2) + 1); +} + +// What a row store's spilled rows cost outside its own arrays: the map node per entry (key, mapped +// value and ~24 bytes of std::unordered_map node and bucket overhead). Shared for the same reason as +// the capacity rule above -- the node-overhead estimate is a single number that must not be corrected +// in one store and not the other, which would skew operator_memory_breakdown() for one backend only. +template +[[nodiscard]] inline auto spilled_rows_bytes(const OverflowMap &overflow) -> size_t { + return overflow.size() * (sizeof(typename OverflowMap::mapped_type) + sizeof(size_t) + 24); +} + class TermIndexCeilingReached : public std::runtime_error { public: using std::runtime_error::runtime_error; diff --git a/cpp/monoprop/detail/operator/SparseRowStore.h b/cpp/monoprop/detail/operator/SparseRowStore.h new file mode 100644 index 00000000..bf4c42e2 --- /dev/null +++ b/cpp/monoprop/detail/operator/SparseRowStore.h @@ -0,0 +1,765 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/RowHashTable.h" + +// Logical mode count at or above which the sparse rows are the cheaper backend. Build-time and not a +// runtime knob because what moves the crossover is the target ISA, which is fixed when the translation +// unit is compiled: dense costs one pass per storage word and sparse is flat in the width, so without a +// vector popcount the dense pass degrades an order of magnitude sooner. Set from CMake off the arch +// flags actually emitted (see the top-level CMakeLists for the measured values). +// +// Hard error rather than a fallback default: the value is a usage requirement of monoprop-objs, so a +// translation unit reaching here without it did not inherit that target's requirements, and any +// fallback would differ from the value the rest of the library was compiled with. kMinModes reaches +// inline definitions (MPOperator, with_store), so under LTO that disagreement is an ODR violation +// resolving to one arbitrary answer -- a silently wrong backend choice rather than a build failure. +#ifndef monoprop_SPARSE_ROW_MIN_MODES +#error \ + "monoprop_SPARSE_ROW_MIN_MODES is undefined: link against the monoprop-objs target rather than adding its include paths by hand." +#endif + +namespace monoprop::detail { + +// A sparse row's two storage types, at namespace scope rather than inside the store: the algebra that +// reads rows (CodesAlgebra.h) must not depend on the container that owns them. +using RowMode = uint16_t; +using RowCodes = uint64_t; + +// Two bits per slot in one RowCodes. A wider row is representable in the mode lanes but would put the +// algebra back on a multi-word loop, which is the whole cost the support form removes. +inline constexpr size_t kRowMaxSlots = 32; +inline constexpr RowCodes kRowLoBits = 0x5555555555555555ULL; // bit 2j of every slot + +// Bit 2j of each slot: set iff slot j is occupied at all. popcount is n, the support measure. +[[nodiscard]] constexpr auto row_occupied_bits(RowCodes codes) noexcept -> RowCodes { + return (codes | (codes >> 1)) & kRowLoBits; +} +// Bit 2j of each slot: set iff slot j holds both of its positions. popcount is d. +[[nodiscard]] constexpr auto row_paired_bits(RowCodes codes) noexcept -> RowCodes { + return codes & (codes >> 1) & kRowLoBits; +} +[[nodiscard]] constexpr auto row_slot_count(RowCodes codes) noexcept -> size_t { + return static_cast(std::popcount(row_occupied_bits(codes))); +} + +// Non-owning view of one row: ascending mode lanes plus the codes word. The lane array is only read +// below num_slots(), which the codes word determines -- so a view stays valid over a padded row and +// carries no length of its own. It borrows the store's arrays, so it must not outlive them, and a row +// mutation invalidates it the way an iterator would. +struct SparseRow { + const RowMode *modes = nullptr; + RowCodes codes = 0; + + [[nodiscard]] auto num_slots() const noexcept -> size_t { return row_slot_count(codes); } + [[nodiscard]] auto mode(size_t j) const noexcept -> size_t { return static_cast(modes[j]); } + // The 2-bit field of slot j: 0b01 is physical position 2*mode alone, 0b10 is 2*mode+1 alone, 0b11 + // the paired mode. + [[nodiscard]] auto code(size_t j) const noexcept -> unsigned int { + return static_cast((codes >> (2 * j)) & 0b11U); + } +}; + +// A row *key* that may be too wide for a codes word: `spilled` non-null means the key is that dense +// monomial and `row` is unread. These are the two shapes a stored row already has, and a query needs +// both for the same reason a stored row does -- a query is M ⊕ G, and a fully paired product escapes +// the cutoff, so nothing bounds its support. +// +// Deliberately not folded into SparseRow, which is what the per-term algebra reads: that one is a view +// of something that fits, and a branch on every read of it is the cost the support form exists to +// avoid. The branch belongs here, on the probe path, where it runs once per query. +template +struct SparseRowKey { + SparseRow row; + const Bitset *spilled = nullptr; + + [[nodiscard]] auto is_spilled() const noexcept -> bool { return spilled != nullptr; } +}; + +// Visits a *dense* monomial as (mode, code) slots, ascending: the same sequence a SparseRow over the +// same monomial yields. Positions arrive ascending, so a mode's two positions are adjacent and one pass +// closes each slot before opening the next. +template +inline auto for_each_mode_slot(const Bitset &mono, Fn &&fn) -> void { + size_t pos = mono.find_first(); + while (pos < mono.size()) { + const size_t mode = pos >> 1; + unsigned int code = 1U << (pos & 1U); + pos = mono.find_next(pos); + if (pos < mono.size() && (pos >> 1) == mode) { + code |= 1U << (pos & 1U); + pos = mono.find_next(pos); + } + fn(mode, code); + } +} + +// Occupied modes in a dense monomial, via the same slot walk as sparse_row_hash/dense_row_equals below -- +// what a spilled row's occupied_modes() reports, and what a per-gate generator's mode count also needs +// (see sparse_record_capacity in layer_build/TermProduct.h). +template +[[nodiscard]] inline auto occupied_mode_count(const Bitset &mono) -> size_t { + size_t n = 0; + for_each_mode_slot(mono, [&n](size_t, unsigned int) { ++n; }); + return n; +} + +// The row hash, as an accumulator over (mode, code) slots. One definition with two walkers -- a sparse +// row and a dense monomial -- because a keyed store must hold both and hash them identically: a fully +// paired term escapes the cutoff, so a row can occupy more modes than any codes word holds and has to +// spill to the dense side map, where it still needs to be findable. +// +// That requirement is why the mix is *sequential* rather than an XOR-fold of slot-indexed terms, which +// is what the plan's "fixed-width mix over the padded row" would have been. Sequential mixing is +// positional without packing a slot index into the mixed word, so it does not care how many slots there +// are -- and it depends on neither the row capacity nor the padding, so two stores tuned to different +// capacities agree, which matters because the hash decides probe order. +class SparseRowHasher { +public: + auto add(size_t mode, unsigned int code) noexcept -> void { + h_ = monoprop::splitmix_finalize(h_ ^ ((static_cast(mode) << 2) | code)); + } + [[nodiscard]] auto value() const noexcept -> size_t { return static_cast(h_); } + +private: + // Nonzero, so an empty row does not hash to zero and every slot count starts from a mixed state. + uint64_t h_ = 0x9E3779B97F4A7C15ULL; +}; + +[[nodiscard]] inline auto sparse_row_hash(const SparseRow &row) noexcept -> size_t { + SparseRowHasher hasher; + const size_t n = row.num_slots(); + for (size_t j = 0; j < n; ++j) { + hasher.add(row.mode(j), row.code(j)); + } + return hasher.value(); +} + +template +[[nodiscard]] inline auto sparse_row_hash(const Bitset &mono) noexcept -> size_t { + SparseRowHasher hasher; + for_each_mode_slot(mono, [&hasher](size_t mode, unsigned int code) { hasher.add(mode, code); }); + return hasher.value(); +} + +// Dispatches to whichever shape the key holds, so a batch of keys hashes identically whether or not any +// of them spilled. The two arms must agree with the store's own row hash, which is what makes a spilled +// row findable by either form. +template +[[nodiscard]] inline auto sparse_row_hash(const SparseRowKey &key) noexcept -> size_t { + return key.is_spilled() ? sparse_row_hash(*key.spilled) : sparse_row_hash(key.row); +} + +// Whether a dense monomial and a sparse row hold the same slots, without materializing either. Used +// where one side is a spilled row (no codes word) and the other is a query. +template +[[nodiscard]] inline auto dense_row_equals(const Bitset &mono, const SparseRow &row) -> bool { + const size_t n = row.num_slots(); + size_t j = 0; + bool equal = true; + for_each_mode_slot(mono, [&equal, &j, &n, &row](size_t mode, unsigned int code) { + if (!equal) { + return; + } + if (j >= n || row.mode(j) != mode || row.code(j) != code) { + equal = false; + return; + } + ++j; + }); + return equal && j == n; +} + +// Writes a sparse row's occupied slots into `mono`, which must already be cleared -- a fresh +// Monomial, or one a caller reset itself before refilling it. The shared body behind every +// dense materialization of a SparseRow. +template +inline auto fill_from_sparse_row(const SparseRow &row, Bitset &mono) -> void { + const size_t n = row.num_slots(); + for (size_t j = 0; j < n; ++j) { + const unsigned int code = row.code(j); + if ((code & 1U) != 0U) { + mono.set(2 * row.mode(j)); + } + if ((code & 2U) != 0U) { + mono.set((2 * row.mode(j)) + 1); + } + } +} + +// Materializes a sparse row as a fresh dense monomial. +template +[[nodiscard]] inline auto sparse_row_to_monomial(const SparseRow &row) -> Bitset { + Bitset mono; + fill_from_sparse_row(row, mono); + return mono; +} + +// Operator-term store in support form: each row is a fixed-width list of the *modes* it occupies plus +// one word holding two bits per occupied mode. It is the third backend behind the four TypeAliases.h +// row accessors, alongside std::vector and OperatorIndex, and agrees with both through them +// (cpp/tests/row_accessor_tests.cpp). +// +// Layout, structure-of-arrays: modes_ is `slots_per_row_` ModeT lanes per row, ascending, padded with +// kPadLane; the codes array is one word per row, at the narrowest of three widths that holds +// 2 * slots_per_row_ bits (see CodesWidth). The two live in separate arrays on purpose -- the cutoff and +// pairing algebra reads only the codes, so evaluating it over a run of rows is a sequential walk that +// never touches a mode list, and narrowing that array puts proportionally more rows on each line of it. +// +// A codes word packs slot j (the j-th occupied mode, ascending) into bits 2j and 2j+1: bit 2j marks physical +// position 2*mode, bit 2j+1 marks 2*mode+1. The whole cutoff algebra follows from that one word -- +// with occupied = (codes | codes>>1) & 0x5555..., paired = codes & (codes>>1) & 0x5555..., +// n = popcount(occupied) and d = popcount(paired) give or_sum = n, popcount_sum = n + d and +// xor_sum = n - d, independent of the storage width. popcount() below is the first consumer; the rest +// arrives with the algebra port. +// +// Rows wider than slots_per_row_ spill losslessly to a side map. They are not a corner case to be +// ruled out by sizing: a fully-paired term escapes the cutoff (xor_sum == 0 is kept unconditionally), +// so support is genuinely unbounded no matter what the cutoff is. They are rare -- ~0.07% of rows on +// production models, per MPOperator.h -- and an all-0b11 row needs only its mode list, so a second +// cheap row kind is available if that ever stops being true. +// +// The keyless index over the rows is the shared RowHashTable, the same one OperatorIndex uses, so both +// stores produce the same slot layout for the same insertion sequence. What differs is only what a key +// is: rows hash through sparse_row_hash and confirm through a codes compare plus a lane memcmp, where +// OperatorIndex hashes a whole Bitset. That hash is *not* the dense one, so a store swap changes probe +// order, MPI owner routing and therefore floating-point accumulation order -- the deliberate +// re-baseline, not a regression. +// +// static_assert cannot express it, so: SparseRowStore is interchangeable with OperatorIndex through the +// TypeAliases.h accessors and through find/emplace/bulk_insert/find_batch, and cpp/tests are what hold +// that. It is not a subclass of anything and nothing dispatches on it. +// +// Single-writer, like OperatorIndex: one partition, one thread; parallelism is cross-partition. +template +class SparseRowStore { +public: + using value_type = Monomial; + using key_type = Monomial; + using mapped_type = size_t; + using ModeT = RowMode; + using CodesT = RowCodes; + + static constexpr size_t kMaxSlots = kRowMaxSlots; + static constexpr size_t kDefaultSlots = 8; + + // The top two ModeT values are markers, so a valid mode index is at most kPadLane - 2. kPadLane + // fills the unused lanes of a short row (fixed, so two equal rows have equal lanes); kOverflowLane + // sits in lane 0 of a spilled row, where it cannot be confused with the empty row's kPadLane. + static constexpr ModeT kPadLane = std::numeric_limits::max(); + static constexpr ModeT kOverflowLane = static_cast(kPadLane - 1); + static constexpr size_t kMaxModes = static_cast(kOverflowLane); // exclusive bound + + static constexpr size_t kIndexCeiling = RowHashTable::kIndexCeiling; + static constexpr size_t kNotFound = RowHashTable::kNotFound; + + // The Stage 3 crossover, as a predicate rather than a bare number so the rule has one home. + static constexpr size_t kMinModes = monoprop_SPARSE_ROW_MIN_MODES; + [[nodiscard]] static constexpr auto preferred_for_modes(size_t num_modes) noexcept -> bool { + return num_modes >= kMinModes; + } + + // kPadLane and kOverflowLane take the top two ModeT values, so the widest system this backend can + // represent stops two modes short of ModeT's range -- a static_assert rather than a throw, since the + // width is a template argument. + static_assert(NumModes <= kMaxModes, "SparseRowStore: NumModes exceeds what a mode lane can address"); + + // slots_per_row is the per-row mode capacity -- any value is correct, since over-long rows spill; + // size it from CutoffEvaluator::max_mode_bound(). + explicit SparseRowStore(size_t slots_per_row = kDefaultSlots) + : slots_per_row_(std::clamp(slots_per_row, 1, kMaxSlots)), + codes_width_(codes_width_for(slots_per_row_)) {} + + SparseRowStore(const SparseRowStore &) = delete; + SparseRowStore &operator=(const SparseRowStore &) = delete; + SparseRowStore(SparseRowStore &&) = delete; + SparseRowStore &operator=(SparseRowStore &&) = delete; + +private: + // Storage width of the codes array. A codes word carries two bits per slot, so a store sized from a + // cutoff bound only ever sets 2 * slots_per_row_ of the 64 bits a CodesT has -- 12 at cutoff 6 and 16 + // at cutoff 8, the two shipping models. Narrowing the storage recovers the rest, 6 bytes per row at + // both, which was this backend's whole per-row gap to OperatorIndex's (1 + inline_width) payloads. + // + // Only the array narrows. Every reader still sees a CodesT, zero-extended on load, so CodesAlgebra.h, + // sparse_row_hash, SparseRow and the cutoff algebra are untouched -- and with them the term set, the + // values and the probe order. Rows are payload, never a hash input and never serialized, so nothing + // here is visible to a baseline diff; the footprint gate is the memory_bytes() case in + // cpp/tests/sparse_row_store_tests.cpp. + enum class CodesWidth : uint8_t { Narrow, Medium, Wide }; + + [[nodiscard]] static constexpr auto codes_width_for(size_t slots) noexcept -> CodesWidth { + if ((2 * slots) <= 16) { + return CodesWidth::Narrow; + } + if ((2 * slots) <= 32) { + return CodesWidth::Medium; + } + return CodesWidth::Wide; + } + + // Sits mid-class for the reason OperatorIndex::with_rows does: a deduced return type is not available + // to a caller that appears earlier in the class body. + // + // Binds the codes storage type for one call. codes_width_ is fixed at construction, so the branch is + // a load-and-test on a member that never changes -- predicted, and one per row read rather than per + // slot. Not hoisted into the store type, for the same reason the row payload is not: the codes width + // is not part of the seam the scan is templated on (see with_store in MPOperator), and making it so + // would triple every downstream instantiation to save a predicted branch. + template + [[gnu::always_inline]] auto with_codes(this Self &&self, F &&f) -> decltype(auto) { + switch (self.codes_width_) { + case CodesWidth::Narrow: + return f(self.codes16_); + case CodesWidth::Medium: + return f(self.codes32_); + default: + return f(self.codes64_); + } + } + + [[nodiscard]] auto load_codes(size_t i) const -> CodesT { + return with_codes([i](const auto &codes) -> CodesT { return static_cast(codes[i]); }); + } + + // The narrowing cast is exact rather than checked: a row reaches here only after it is known to fit + // slots_per_row_ slots, and slot j occupies bits 2j and 2j+1, so no bit at or above 2 * slots_per_row_ + // is ever set. The assert is what holds that when a caller hands over a SparseRow it built itself. + auto store_codes(size_t i, CodesT codes) -> void { + with_codes([i, codes](auto &store) { + using ElemT = typename std::remove_cvref_t::value_type; + assert(codes == static_cast(static_cast(codes)) && "codes word wider than its storage"); + store[i] = static_cast(codes); + }); + } + + auto resize_codes(size_t n) -> void { + with_codes([n](auto &codes) { codes.resize(n); }); + } + auto reserve_codes(size_t n) -> void { + with_codes([n](auto &codes) { codes.reserve(n); }); + } + [[nodiscard]] auto codes_capacity() const -> size_t { + return with_codes([](const auto &codes) { return codes.capacity(); }); + } + // The three arrays differ only in element size, so everything that counts bytes rather than reading a + // word is plain arithmetic off this and needs no type bound. + [[nodiscard]] auto codes_bytes() const noexcept -> size_t { + switch (codes_width_) { + case CodesWidth::Narrow: + return sizeof(uint16_t); + case CodesWidth::Medium: + return sizeof(uint32_t); + default: + return sizeof(CodesT); + } + } + +public: + [[nodiscard]] static constexpr auto num_bits() noexcept -> size_t { return 2 * NumModes; } + [[nodiscard]] auto slots_per_row() const noexcept -> size_t { return slots_per_row_; } + // The backend-neutral spelling of the line above, so a caller holding either store asks the same + // question of both (see MPOperator::row_width). + [[nodiscard]] auto row_width() const noexcept -> size_t { return slots_per_row_; } + [[nodiscard]] auto size() const noexcept -> size_t { return size_; } + + // Called only on an idle store, so it needs no synchronization. + [[nodiscard]] auto clone() const -> std::unique_ptr { + auto out = std::make_unique(slots_per_row_); + out->modes_ = modes_; + // Exactly one of the three is non-empty, and out shares slots_per_row_ so it shares codes_width_. + out->codes16_ = codes16_; + out->codes32_ = codes32_; + out->codes64_ = codes64_; + out->size_ = size_; + out->overflow_ = overflow_; + out->table_ = table_; // RowHashTable is rule-of-zero copyable; a plain copy preserves slot order exactly. + return out; + } + + // Same term set at a different slots_per_row_, e.g. after a cutoff change moves the bound rows are + // sized from. Every row's monomial is re-flowed through set() at the new stride, which decides + // inline-vs-overflow the same way a fresh insert would; the hash index is copied as-is, since + // fold_hash (via sparse_row_hash) depends only on the monomial, never on slots_per_row_, so no rehash + // is needed. Row index i is preserved for every row -- load-bearing, since callers key op_coeffs, + // state_rows_/state_vals_ and the evolution graph by this same index. + [[nodiscard]] auto resized(size_t new_slots_per_row) const -> std::unique_ptr { + auto out = std::make_unique(new_slots_per_row); + out->modes_.resize(size_ * out->slots_per_row_); + out->resize_codes(size_); + out->size_ = size_; + // Reflow via view()/overflow_ directly rather than row(i): row() would materialize a fresh + // Bitset from the slots and set() would immediately re-walk it to rebuild them, a double pass + // this store's own non-allocating set(SparseRow) / set(value_type) overloads make unnecessary. + for (size_t i = 0; i < size_; ++i) { + if (spilled(i)) { + out->set(i, overflow_.at(i)); + } + else { + out->set(i, view(i)); + } + } + out->table_ = table_; // RowHashTable is rule-of-zero copyable; a plain copy preserves slot order exactly. + return out; + } + + auto reserve(size_t n) -> void { + reserve_rows_(n); + table_.reserve(n); + } + + // Returns the pre-growth size (the caller's insert base). Growth is geometric (1.5x), never + // exact-fit: an exact fit would realloc the whole operator every layer. Rows only -- the table grows + // on its own load factor, and pre-sizing it per layer would rehash for nothing. + auto grow_rows_geometric(size_t n) -> size_t { + const size_t base = size_; + if (capacity() < base + n) { + reserve_rows_(geometric_row_capacity(base, n, capacity())); + } + // Default-init grow, not a zeroing resize: every freshly grown row is overwritten by set() + // before any read, so a tail zero-fill would be wasted bandwidth. + modes_.resize((base + n) * slots_per_row_); + resize_codes(base + n); + size_ = base + n; + return base; + } + + auto push_back(const value_type &mono) -> void { set(grow_rows_geometric(1), mono); } + + // Row i may be grown-but-uninitialized or hold a prior value, so nothing in the row is pre-read; a + // stale overflow entry at i, if any, is dropped. + auto set(size_t i, const value_type &mono) -> void { + ModeT *lanes = &modes_[i * slots_per_row_]; + CodesT codes = 0; + size_t used = 0; + bool overflows = false; + // The slot walk is shared with the hash, so the two cannot disagree about what a row's slots are. + // Lanes come out ascending because the walk is. + for_each_mode_slot(mono, [&overflows, &used, &lanes, &codes, this](size_t mode, unsigned int code) { + if (overflows) { + return; + } + if (used == slots_per_row_) { + overflows = true; + return; + } + lanes[used] = static_cast(mode); + codes |= static_cast(code) << (2 * used); + ++used; + }); + if (overflows) { + lanes[0] = kOverflowLane; + store_codes(i, 0); + overflow_[i] = mono; + return; + } + if (!overflow_.empty()) { + overflow_.erase(i); + } + for (size_t j = used; j < slots_per_row_; ++j) { + lanes[j] = kPadLane; + } + store_codes(i, codes); + } + + // The row form of set(), and the write the support form exists for: the lanes are already ascending + // and the codes word already says what each holds, so this copies `n` lanes and one word where the + // dense overload walks the monomial's storage words. + // + // A row wider than this store's capacity still has to spill, and a spilled row is held densely, so + // that arm materializes. It cannot be asserted away: the capacity is sized from the cutoff and a + // fully paired term escapes the cutoff. + auto set(size_t i, const SparseRow &row) -> void { + const size_t n = row.num_slots(); + // Contiguity from slot 0 is the representation's invariant -- num_slots() counts occupied slots + // and the lanes are read from 0 -- so a row with a hole would silently lose its high slots here. + assert((n >= kRowMaxSlots || (row.codes >> (2 * n)) == 0) && "SparseRow slots must be contiguous from slot 0"); + ModeT *lanes = &modes_[i * slots_per_row_]; + if (n > slots_per_row_) { + lanes[0] = kOverflowLane; + store_codes(i, 0); + overflow_[i] = to_monomial_(row); + return; + } + // Hygiene, not correctness: spilled() reads lane 0, so a stale entry here is already unreachable + // -- it would just keep a monomial alive for the store's lifetime. + if (!overflow_.empty()) { + overflow_.erase(i); + } + if (n != 0) { + std::memcpy(lanes, row.modes, n * sizeof(ModeT)); + } + // Padding is load-bearing for the empty row and only for it: with n == 0 nothing above writes a + // lane, so lane 0 would keep a previous occupant's kOverflowLane and the row would read as spilled. + for (size_t j = n; j < slots_per_row_; ++j) { + lanes[j] = kPadLane; + } + store_codes(i, row.codes); + } + + // Whichever shape the key holds. The spilled arm is the dense set(), so a key that arrived too wide + // for a codes word lands in the side map exactly as the dense path would have put it. + auto set(size_t i, const SparseRowKey<2 * NumModes> &key) -> void { + if (key.is_spilled()) { + set(i, *key.spilled); + return; + } + set(i, key.row); + } + + [[nodiscard]] auto row(size_t i) const -> value_type { + if (spilled(i)) { + return overflow_.at(i); + } + return sparse_row_to_monomial<2 * NumModes>(view(i)); + } + + // Ascending, matching the dense backends: slots are stored ascending in the mode, and within a mode + // position 2*mode precedes 2*mode+1. + template + auto for_each_position(size_t i, Fn &&fn) const -> void { + if (spilled(i)) { + const auto &m = overflow_.at(i); + for (size_t b = m.find_first(); b < m.size(); b = m.find_next(b)) { + fn(b); + } + return; + } + for_each_slot(i, [&fn](size_t mode, unsigned int code) { + if ((code & 1U) != 0U) { + fn(2 * mode); + } + if ((code & 2U) != 0U) { + fn((2 * mode) + 1); + } + }); + } + + // Visits (mode, code) per occupied slot, ascending in the mode; code is the 2-bit field, so 0b01 is + // position 2*mode alone, 0b10 is 2*mode+1 alone and 0b11 is the paired mode. Row i must not be + // spilled -- a spilled row has no slots, and its lane 0 marker would read as a mode. + template + auto for_each_slot(size_t i, Fn &&fn) const -> void { + assert(!spilled(i) && "SparseRowStore::for_each_slot on a spilled row"); + const ModeT *lanes = &modes_[i * slots_per_row_]; + const CodesT codes = load_codes(i); + for (size_t j = 0; j < slots_per_row_ && lanes[j] != kPadLane; ++j) { + fn(static_cast(lanes[j]), static_cast((codes >> (2 * j)) & 0b11U)); + } + } + + // The row's codes word. Meaningless for a spilled row -- ask spilled(i) first; the algebra port + // will need the same guard, which is why the spill is kept rare rather than made general. + [[nodiscard]] auto codes(size_t i) const -> CodesT { return load_codes(i); } + + // What the codes algebra reads. Borrows this store's arrays, so it is invalidated by anything that + // reallocates them (grow_rows_geometric, reserve) or rewrites row i; row i must not be spilled. + [[nodiscard]] auto view(size_t i) const -> SparseRow { + assert(!spilled(i) && "SparseRowStore::view on a spilled row"); + return SparseRow{&modes_[i * slots_per_row_], load_codes(i)}; + } + + [[nodiscard]] auto spilled(size_t i) const -> bool { return modes_[i * slots_per_row_] == kOverflowLane; } + + // Occupied modes -- the support measure, or_sum. + [[nodiscard]] auto slot_count(size_t i) const -> size_t { + if (spilled(i)) { + return occupied_modes(overflow_.at(i)); + } + return row_slot_count(load_codes(i)); + } + + // Set bits -- the length measure, popcount_sum = n + d straight off the codes word. + [[nodiscard]] auto popcount(size_t i) const -> size_t { + if (spilled(i)) { + return overflow_.at(i).count(); + } + const CodesT codes = load_codes(i); + return row_slot_count(codes) + static_cast(std::popcount(row_paired_bits(codes))); + } + + // --- the keyless index over those rows --------------------------------------------------------- + // + // A key is a SparseRow or a Bitset, and both hash through sparse_row_hash, so the two are + // interchangeable at a call site. Prefer the row: it is what the scan holds, and it is the only form + // that needs no slot walk to hash. The Bitset form is what a caller still holding a monomial uses. + + auto find(const SparseRow &key) const -> std::optional { return find_hashed_(key); } + auto find(const key_type &key) const -> std::optional { return find_hashed_(key); } + auto find(const SparseRowKey<2 * NumModes> &key) const -> std::optional { return find_hashed_(key); } + + // Insert-or-no-op. The row at `value` must already be written -- the confirm reads it. + template + auto emplace(const Key &key, mapped_type value) -> void { + table_.emplace(fold_hash(key), value, [&key, this](size_t i) { return row_eq_key(i, key); }); + } + + // Insert n distinct rows with consecutive indices [base, base+n). Rows must already be written. + template + auto bulk_insert(size_t n, mapped_type base, KeyFn &&key_at) -> void { + table_.insert_distinct_range(base, n, [&key_at](size_t k) { return fold_hash(key_at(k)); }); + } + + // out[i] = row index of keys[i], or kNotFound. Same result as n find() calls; see + // RowHashTable::find_batch for why the row prefetch sits between probe and confirm. Both arrays are + // prefetched: the confirm reads the codes word first and the lanes only if it matches, but they are + // separate allocations and so separate cache misses. + template + auto find_batch(const Key *keys, size_t n, size_t *out) const -> void { + table_.find_batch( + keys, + n, + out, + [](const Key &key) { return fold_hash(key); }, + [this](size_t i) { + with_codes([i](const auto &codes) { __builtin_prefetch(&codes[i], 0, 0); }); + __builtin_prefetch(&modes_[i * slots_per_row_], 0, 0); + }, + [this](size_t i, const Key &key) { return row_eq_key(i, key); }); + } + + // Rows in table order (for_each_slot walks the slot array, i.e. hash/probe order, not ascending row + // index -- see the class comment above), as fn(row_index). Not the row itself: a spilled row has no + // view, so what a caller wants off the index is the index. + template + auto for_each_index(Fn &&fn) const -> void { + table_.for_each_slot([&fn](TermIndex idx, uint32_t) { fn(static_cast(idx)); }); + } + + // OperatorIndex's signature, fn(monomial, row_index), so the two stores are interchangeable at the + // one call site that wants both. Materializes each row, which for_each_index does not -- prefer that + // where the index alone will do. + template + auto for_each(Fn &&fn) const -> void { + for_each_index([&fn, this](size_t i) { fn(row(i), i); }); + } + + [[nodiscard]] auto indexed_count() const noexcept -> size_t { return table_.count(); } + + [[nodiscard]] auto index_estimated_memory_bytes() const -> size_t { + return sizeof(SparseRowStore) + table_.slot_bytes(); + } + + [[nodiscard]] auto memory_bytes() const -> size_t { + return (modes_.capacity() * sizeof(ModeT)) + (codes_capacity() * codes_bytes()) + spilled_rows_bytes(overflow_); + } + + // Diagnostic: the part of memory_bytes() that is unused geometric-growth capacity. + [[nodiscard]] auto slack_bytes() const -> size_t { + const size_t lanes = modes_.capacity() - std::min(modes_.capacity(), size_ * slots_per_row_); + const size_t words = codes_capacity() - std::min(codes_capacity(), size_); + return (lanes * sizeof(ModeT)) + (words * codes_bytes()); + } + + // Slot count for a cutoff bound in modes (CutoffEvaluator::max_mode_bound()), clamped to what one + // codes word holds. A bound above kMaxSlots is not an error: the rows that exceed it spill. + [[nodiscard]] static auto slots_for_bound(size_t mode_bound) noexcept -> size_t { + return std::clamp(mode_bound, 1, kMaxSlots); + } + + // Slot count for a row in flight rather than a row at rest: a product occupies up to the term's modes + // plus the generator's, so a scan scratch row and a wire record both need the cutoff's mode bound plus + // the widest generator's locality. Every rank derives this from the same circuit and cutoff, so they + // agree on it without communication -- which is what lets it fix a wire stride. + [[nodiscard]] static auto scratch_slots_for(size_t mode_bound, size_t max_generator_modes) noexcept -> size_t { + return std::clamp(mode_bound + max_generator_modes, 1, kMaxSlots); + } + +private: + // Spilled rows have no codes word, so their support is counted the dense way. + [[nodiscard]] static auto occupied_modes(const value_type &mono) -> size_t { return occupied_mode_count(mono); } + + // The 32-bit fold the table stores as its equality pre-filter, over the full-width row hash. + template + static auto fold_hash(const Key &key) noexcept -> uint32_t { + return RowHashTable::fold(sparse_row_hash(key)); + } + + template + auto find_hashed_(const Key &key) const -> std::optional { + return table_.find(fold_hash(key), [&key, this](size_t i) { return row_eq_key(i, key); }); + } + + // The find confirm, against a row query. Codes first: one word compare rejects nearly every + // pre-filter false positive, and it is what fixes the lane compare's length -- equal codes means + // equal slot counts, so only the occupied lanes can differ. Comparing the padded width instead (as + // the plan sketched) would be the same cost but would silently mismatch any query a caller left + // unpadded, and the padding is capacity-dependent where a key must not be. + [[nodiscard]] auto row_eq_key(size_t i, const SparseRow &key) const -> bool { + if (spilled(i)) { + return dense_row_equals(overflow_.at(i), key); + } + if (load_codes(i) != key.codes) { + return false; + } + const size_t k = row_slot_count(key.codes); + return k == 0 || std::memcmp(&modes_[i * slots_per_row_], key.modes, k * sizeof(ModeT)) == 0; + } + + // The same against a monomial query, which is the only form that can match a spilled row exactly. + [[nodiscard]] auto row_eq_key(size_t i, const key_type &key) const -> bool { + if (spilled(i)) { + return overflow_.at(i) == key; + } + return dense_row_equals(key, view(i)); + } + + [[nodiscard]] auto row_eq_key(size_t i, const SparseRowKey<2 * NumModes> &key) const -> bool { + return key.is_spilled() ? row_eq_key(i, *key.spilled) : row_eq_key(i, key.row); + } + + // A row at this store's width. Only the spill arms need it: everything else reads slots in place. + [[nodiscard]] auto to_monomial_(const SparseRow &row) const -> value_type { + return sparse_row_to_monomial<2 * NumModes>(row); + } + + auto reserve_rows_(size_t n) -> void { + modes_.reserve(n * slots_per_row_); + reserve_codes(n); + } + + [[nodiscard]] auto capacity() const -> size_t { return codes_capacity(); } + + DefaultInitVector modes_ = {}; + // Exactly one is ever non-empty, selected by codes_width_ -- the same one-live-arm shape + // OperatorIndex uses for its row payload. + DefaultInitVector codes16_ = {}; + DefaultInitVector codes32_ = {}; + DefaultInitVector codes64_ = {}; + size_t size_ = 0; + size_t slots_per_row_ = kDefaultSlots; + // Declared after slots_per_row_: the constructor derives it from the clamped value. + CodesWidth codes_width_ = codes_width_for(kDefaultSlots); + // Lossless side-map for rows occupying more than slots_per_row_ modes. + std::unordered_map overflow_ = {}; + RowHashTable table_ = {}; +}; + +} // namespace monoprop::detail diff --git a/cpp/tests/RandomMonomial.h b/cpp/tests/RandomMonomial.h new file mode 100644 index 00000000..222e59fc --- /dev/null +++ b/cpp/tests/RandomMonomial.h @@ -0,0 +1,49 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include +#include + +#include "monoprop/core/Monomial.h" + +namespace test_utils { + +// A random monomial over `NumModes` modes occupying at most `max_slots` of them, each with a uniformly +// random non-empty code (one Majorana of the mode, the other, or the pair). +// +// One definition, deliberately: this is the input distribution of the whole randomized sparse/codes test +// surface, and what it biases toward -- paired slots in particular, which is what drives spills and +// product overflow -- decides what those tests actually cover. A per-file copy would let one of them be +// tuned and the rest silently left behind. Kept out of TestUtilities.h, which pulls in Boost.Test, +// MonomialPropagator and MPI; the files that want this want nothing else. +template +inline auto random_monomial(std::mt19937_64 &rng, size_t max_slots) -> monoprop::Monomial { + monoprop::Monomial mono; + const size_t occupied = rng() % (max_slots + 1); + for (size_t k = 0; k < occupied; ++k) { + const size_t mode = rng() % NumModes; + const auto code = 1U + static_cast(rng() % 3U); + if ((code & 1U) != 0U) { + mono.set(2 * mode); + } + if ((code & 2U) != 0U) { + mono.set((2 * mode) + 1); + } + } + return mono; +} + +} // namespace test_utils diff --git a/cpp/tests/TestOperator.h b/cpp/tests/TestOperator.h new file mode 100644 index 00000000..e3d4359e --- /dev/null +++ b/cpp/tests/TestOperator.h @@ -0,0 +1,47 @@ +// Copyright 2026 Algorithmiq +// +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#pragma once + +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/MPOperator.h" +#include "monoprop/detail/operator/RowAccess.h" + +namespace test_utils { + +// An MPOperator holding `terms` whose rows are also *findable*: append_term writes a row and nothing +// else, so find()/find_batch see nothing until the hash index is populated, which only the +// insert_absent_terms path does. That sequence is the one correct incantation for "an operator a resolve +// can look terms up in", so it lives here rather than being copied into each test that needs one. +template +inline auto indexed_operator(const monoprop::MonomialList &terms, + monoprop::Basis basis = monoprop::Basis::Majorana) + -> monoprop::detail::MPOperator { + monoprop::detail::MPOperator op; + op.basis = basis; + op.with_store([&](auto &rows) { + monoprop::detail::insert_absent_terms( + op, + rows, + terms.size(), + [&](size_t k) -> const monoprop::Monomial & { return terms[k]; }, + [&](size_t k, size_t base) { monoprop::assign_row(rows, base + k, terms[k]); }); + }); + return op; +} + +} // namespace test_utils diff --git a/cpp/tests/boost-test.cmake b/cpp/tests/boost-test.cmake index 7d0d52c9..b8f960cd 100644 --- a/cpp/tests/boost-test.cmake +++ b/cpp/tests/boost-test.cmake @@ -4,6 +4,12 @@ set( CACHE STRING "Semicolon-separated list of ranks for MPI test variants" ) +set( + monoprop_MPI_SPARSE_ROWS_TEST_PROCS + "2" + CACHE STRING + "Semicolon-separated list of ranks for the sparse-row-backend MPI test variants (kept separate from monoprop_MPI_TEST_PROCS so growing dense-backend MPI coverage does not silently multiply how many sparse-row mpiexec launches CI pays for)" +) set(_monoprop_mpiexec "${MPIEXEC_EXECUTABLE}") if(NOT _monoprop_mpiexec) @@ -71,6 +77,7 @@ function(discover_tests TARGET) "TEST_LIST=${_TEST_LIST}" -D "CTEST_FILE=${ctest_tests_file}" -D "TEST_ENABLE_MPI_VARIANTS=${_enable_mpi_variants}" -D "TEST_MPI_NUMPROCS=${monoprop_MPI_TEST_PROCS}" -D + "TEST_MPI_SPARSE_ROWS_NUMPROCS=${monoprop_MPI_SPARSE_ROWS_TEST_PROCS}" -D "MPIEXEC_EXECUTABLE=${_monoprop_mpiexec}" -D "MPIEXEC_NUMPROC_FLAG=${_monoprop_mpiexec_numproc_flag}" -D "MPIEXEC_PREFLAGS=${MPIEXEC_PREFLAGS}" -D diff --git a/cpp/tests/boostAddTests.cmake b/cpp/tests/boostAddTests.cmake index 74b5ff00..04d9096e 100644 --- a/cpp/tests/boostAddTests.cmake +++ b/cpp/tests/boostAddTests.cmake @@ -4,6 +4,9 @@ endif() if(NOT DEFINED TEST_MPI_NUMPROCS) set(TEST_MPI_NUMPROCS "2") endif() +if(NOT DEFINED TEST_MPI_SPARSE_ROWS_NUMPROCS) + set(TEST_MPI_SPARSE_ROWS_NUMPROCS "2") +endif() if(TEST_ENABLE_MPI_VARIANTS AND NOT MPIEXEC_EXECUTABLE) message( WARNING @@ -37,6 +40,28 @@ if(TEST_ENABLE_MPI_VARIANTS) list(REMOVE_DUPLICATES _mpi_ranks) endif() +# Validated the same way as _mpi_ranks above, but kept in its own list (TEST_MPI_SPARSE_ROWS_NUMPROCS) +# rather than reusing _mpi_ranks: growing dense-backend rank coverage must not silently multiply how +# many sparse-row mpiexec launches CI pays for. +set(_mpi_sparse_ranks) +if(TEST_ENABLE_MPI_VARIANTS) + if("${TEST_MPI_SPARSE_ROWS_NUMPROCS}" STREQUAL "") + set(TEST_MPI_SPARSE_ROWS_NUMPROCS 2) + endif() + + foreach(_rank IN LISTS TEST_MPI_SPARSE_ROWS_NUMPROCS) + if(NOT _rank MATCHES "^[1-9][0-9]*$") + message( + FATAL_ERROR + "Invalid MPI rank '${_rank}' in TEST_MPI_SPARSE_ROWS_NUMPROCS='${TEST_MPI_SPARSE_ROWS_NUMPROCS}'. Use positive integers." + ) + endif() + endforeach() + + set(_mpi_sparse_ranks ${TEST_MPI_SPARSE_ROWS_NUMPROCS}) + list(REMOVE_DUPLICATES _mpi_sparse_ranks) +endif() + set(extra_args ${TEST_EXTRA_ARGS}) set(properties ${TEST_PROPERTIES}) set(serial_env ${TEST_SERIAL_ENVIRONMENT}) @@ -209,6 +234,28 @@ foreach(LINE ${LINES}) ENVIRONMENT ${serial_env} ) + # Run the same case again with the support-form row backend forced. The suite is below + # SparseRowStore::preferred_for_modes()'s crossover, so the automatic choice would compile + # that backend but never run it, even though it is the one used for wide systems. Running each + # case separately makes any divergence easy to identify. + # + # Keep serial_env for the same reason as the variant above: this is another world-size-1 run + # of the same case, so without it half of `-L serial` would pay the MPI_Init setup cost that + # the other half avoids. + register_variant("${test}_sparse_rows" + COMMAND + "${TEST_EXECUTABLE}" + "--run_test=${test}" + "--report_level=detailed" + "--catch_system_errors=yes" + ${extra_args} + LABELS + serial + sparse-rows + ENVIRONMENT + ${serial_env} + "monoprop_ROW_STORE=sparse" + ) endif() endforeach() @@ -247,6 +294,46 @@ if(TEST_ENABLE_MPI_VARIANTS AND MPIEXEC_EXECUTABLE) "OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1" ) endforeach() + + # MPI counterpart of the "_sparse_rows" serial variant above. The cross-rank resolve inserts absent + # terms into whichever backend is live, and the sparse one is what wide (MPI-scale) systems actually + # resolve to -- so it needs its own multi-rank coverage, not just the single-rank one above. Runs over + # _mpi_sparse_ranks, not _mpi_ranks, so it stays cheap by default regardless of how wide the dense rank + # list grows. + foreach(_mpi_rank IN LISTS _mpi_sparse_ranks) + set(mpi_cmd "${MPIEXEC_EXECUTABLE}") + list( + APPEND mpi_cmd + "${MPIEXEC_NUMPROC_FLAG}" + "${_mpi_rank}" + ) + if(MPIEXEC_PREFLAGS) + list(APPEND mpi_cmd ${MPIEXEC_PREFLAGS}) + endif() + list( + APPEND mpi_cmd + "${TEST_EXECUTABLE}" + "--report_level=detailed" + "--catch_system_errors=yes" + ${extra_args} + ) + if(MPIEXEC_POSTFLAGS) + list(APPEND mpi_cmd ${MPIEXEC_POSTFLAGS}) + endif() + + register_variant("${TEST_TARGET}_mpi_${_mpi_rank}_sparse_rows" + COMMAND + ${mpi_cmd} + LABELS + mpi + "mpi-${_mpi_rank}" + sparse-rows + ENVIRONMENT + "OMPI_ALLOW_RUN_AS_ROOT=1" + "OMPI_ALLOW_RUN_AS_ROOT_CONFIRM=1" + "monoprop_ROW_STORE=sparse" + ) + endforeach() endif() # Create a list of all discovered tests, which users may use to e.g. set diff --git a/cpp/tests/codes_algebra_tests.cpp b/cpp/tests/codes_algebra_tests.cpp new file mode 100644 index 00000000..9dad9e3e --- /dev/null +++ b/cpp/tests/codes_algebra_tests.cpp @@ -0,0 +1,288 @@ +// 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. + +// The differential test between CodesAlgebra.h and the dense implementations it must replace. Every +// function is checked to agree *exactly* -- these are integer and sign quantities, so there is no +// tolerance to spend -- over real fixture monomials and over randomized rows, at storage widths both +// equal to and wider than the logical width. Making the codes form the default is gated on this. + +#include + +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/algebra/CodesAlgebra.h" +#include "monoprop/algebra/MajoranaAlgebra.h" +#include "monoprop/algebra/PauliAlgebra.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "TestData.h" +#include "TestUtilities.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +// Rebuild a dense monomial from a row's mode lanes and an arbitrary codes word, so a codes-side +// transform (pair_swap) can be compared against its dense counterpart. sparse_row_to_monomial is the +// store's own materialization, which is the point: re-implementing it here would leave this oracle +// agreeing with a slot convention the store no longer uses. Only valid where the substituted codes word +// has the same occupancy as the row's own, which is the case for every transform here. +template +auto to_monomial(const SparseRow &row, RowCodes codes) -> Monomial { + return sparse_row_to_monomial<2 * NumModes>(SparseRow{.modes = row.modes, .codes = codes}); +} + +// Which outcomes the comparison actually reached. Every branch of every ported function must be +// exercised by the inputs, or agreement is vacuous -- an interleave phase stubbed to `return 1` agrees +// with the dense version on any input set that happens to contain only even permutations. +struct Seen { + bool paired = false; + bool unpaired = false; + bool y_letters = false; + bool phase_plus = false; + bool phase_minus = false; + bool anticommutes = false; + bool commutes = false; + bool cutoff_kept = false; + bool cutoff_dropped = false; + + auto operator|=(const Seen &o) -> Seen & { + paired |= o.paired; + unpaired |= o.unpaired; + y_letters |= o.y_letters; + phase_plus |= o.phase_plus; + phase_minus |= o.phase_minus; + anticommutes |= o.anticommutes; + commutes |= o.commutes; + cutoff_kept |= o.cutoff_kept; + cutoff_dropped |= o.cutoff_dropped; + return *this; + } +}; + +auto require_discriminating(const Seen &seen) -> void { + BOOST_TEST(seen.paired); + BOOST_TEST(seen.unpaired); + BOOST_TEST(seen.y_letters); + BOOST_TEST(seen.phase_plus); + BOOST_TEST(seen.phase_minus); + BOOST_TEST(seen.anticommutes); + BOOST_TEST(seen.commutes); + BOOST_TEST(seen.cutoff_kept); + BOOST_TEST(seen.cutoff_dropped); +} + +// Every single-row function at once, against the dense version of each. +template +auto check_row(const Monomial &mono, const SparseRow &row, size_t logical_num_modes, Seen &seen) -> void { + const size_t inactive_prefix = NumModes - logical_num_modes; + + const auto dense_sums = cutoff_sums(mono, logical_num_modes); + const auto codes_sums = codes_cutoff_sums(row, inactive_prefix); + BOOST_TEST(codes_sums.or_sum == dense_sums.or_sum); + BOOST_TEST(codes_sums.popcount_sum == dense_sums.popcount_sum); + BOOST_TEST(codes_sums.xor_sum == dense_sums.xor_sum); + + // Both sides of each cutoff's two branches: a cutoff below the term's measure exercises the + // fully-paired escape, one above it the plain comparison. + for (const unsigned int cutoff : {0U, 1U, 2U, 4U, 8U, 64U}) { + const bool kept = codes_length_cutoff(row, cutoff, inactive_prefix); + BOOST_TEST(kept == length_cutoff(mono, cutoff, logical_num_modes)); + BOOST_TEST(codes_support_cutoff(row, cutoff, inactive_prefix) + == support_cutoff(mono, cutoff, logical_num_modes)); + seen.cutoff_kept |= kept; + seen.cutoff_dropped |= !kept; + } + + const bool paired = codes_is_paired(row.codes); + BOOST_TEST(paired == is_paired(mono)); + seen.paired |= paired; + seen.unpaired |= !paired; + + const size_t y = codes_pauli_y_count(row.codes); + BOOST_TEST(y == pauli_y_count(mono)); + seen.y_letters |= y > 0; + + BOOST_TEST((to_monomial(row, codes_pair_swap(row.codes)) == pair_swap(mono))); +} + +// Encode monomials into one store and hand back both the store and the dense originals. All rows are +// pushed before any view is taken: a view borrows the store's arrays, so growth would dangle it. +template +struct Encoded { + std::vector> dense; + SparseRowStore store{SparseRowStore::kMaxSlots}; + + auto add(const Monomial &mono) -> void { + dense.push_back(mono); + store.push_back(mono); + } +}; + +template +auto check_all(Encoded &enc, size_t logical_num_modes) -> Seen { + Seen seen; + BOOST_REQUIRE(enc.dense.size() == enc.store.size()); + for (size_t i = 0; i < enc.dense.size(); ++i) { + BOOST_REQUIRE_MESSAGE(!enc.store.spilled(i), "row " << i << " spilled; the algebra needs a codes word"); + check_row(enc.dense[i], enc.store.view(i), logical_num_modes, seen); + } + // The two-row functions, over every ordered pair for small sets and a stride otherwise: they are + // O(n^2) in the row count and the fixtures carry hundreds of terms. + const size_t n = enc.dense.size(); + const size_t stride = n > 24 ? (n / 24) + 1 : 1; + for (size_t i = 0; i < n; i += stride) { + for (size_t k = 0; k < n; k += stride) { + const auto maj = enc.store.view(i); + const auto gen = enc.store.view(k); + const int phase = codes_interleave_phase(maj, gen); + BOOST_TEST(phase == interleave_phase(enc.dense[i], enc.dense[k])); + const bool anti = codes_pauli_anticommutes(maj, gen); + BOOST_TEST(anti == pauli_anticommutes(enc.dense[i], enc.dense[k])); + seen.phase_plus |= phase > 0; + seen.phase_minus |= phase < 0; + seen.anticommutes |= anti; + seen.commutes |= !anti; + } + } + return seen; +} + +// The fixtures' Hamiltonian keys and generator index lists are the real-world monomials: Hermitian +// Majorana products, so the set includes fully paired rows, which are the inputs both cutoffs treat +// specially. NumModes is the storage width; the fixture's own mode count is the logical one. +template +auto check_fixture(const std::string &name) -> Seen { + const auto data = test_utils::load_case_data(name); + BOOST_REQUIRE(data.num_modes > 0); + BOOST_REQUIRE(NumModes >= data.num_modes); + const size_t max_index = 2 * data.num_modes; + constexpr size_t kMaxSlots = SparseRowStore::kMaxSlots; + + Encoded enc; + for (const auto &[inds, coeff] : data.hamiltonian) { + if (inds.size() > kMaxSlots) { + continue; // would spill; the store's own tests cover that path + } + enc.add(indices_to_bitset_checked(inds, max_index)); + } + for (const auto &inds : data.majoranas) { + if (inds.size() > kMaxSlots) { + continue; + } + enc.add(indices_to_bitset_checked(inds, max_index)); + } + BOOST_REQUIRE_MESSAGE(enc.dense.size() > 1, "fixture " << name << " yielded no monomials to compare"); + return check_all(enc, data.num_modes); +} + +// Randomized rows over one (storage, logical) pair. +template +auto check_random(std::mt19937_64 &rng, size_t logical_num_modes) -> Seen { + Encoded enc; + for (size_t trial = 0; trial < 120; ++trial) { + Monomial mono; + const size_t occupied = rng() % (SparseRowStore::kMaxSlots + 1); + // Every fourth row is forced fully paired: that is the branch both cutoffs short-circuit on and + // the only input is_paired accepts. + const bool force_paired = (trial % 4) == 0; + for (size_t k = 0; k < occupied; ++k) { + const size_t mode = rng() % NumModes; + const unsigned int code = force_paired ? 0b11U : 1U + static_cast(rng() % 3U); + if ((code & 1U) != 0U) { + mono.set(2 * mode); + } + if ((code & 2U) != 0U) { + mono.set((2 * mode) + 1); + } + } + enc.add(mono); + } + return check_all(enc, logical_num_modes); +} + +} // namespace + +// Whole register: storage width equals the logical width, so every mode is active and the codes form +// takes its zero-prefix path. The storage width is a template parameter here, so each fixture is +// instantiated at its own mode count. +BOOST_AUTO_TEST_CASE(codes_algebra_matches_dense_on_fixtures_whole_register) { + Seen seen; + seen |= check_fixture<8>("random_exact.msgpack"); + seen |= check_fixture<12>("lih_fermionic_spin_exact.msgpack"); + seen |= check_fixture<16>("S0_8e8o_majoranic_c6.msgpack"); + seen |= check_fixture<16>("majorana_lattice_layer_30.msgpack"); + require_discriminating(seen); +} + +// The production layout: the logical modes occupy the top of a wider register and the low physical modes +// are inactive. This is the case cutoff_sums applies active_bit_offset for, and the one the codes form +// has to reproduce by dropping a slot prefix. +BOOST_AUTO_TEST_CASE(codes_algebra_matches_dense_on_fixtures_padded_storage) { + Seen seen; + seen |= check_fixture<32>("random_exact.msgpack"); + seen |= check_fixture<32>("lih_fermionic_spin_exact.msgpack"); + seen |= check_fixture<32>("S0_8e8o_majoranic_c6.msgpack"); + seen |= check_fixture<64>("majorana_lattice_layer_30.msgpack"); + require_discriminating(seen); +} + +// Randomized rows reach occupancies and code patterns the fixtures do not: single-position modes in +// every combination, rows at the slot capacity, empty rows, and inactive modes actually populated -- +// which a propagator never produces but the dense functions accept, so the two must still agree. +BOOST_AUTO_TEST_CASE(codes_algebra_matches_dense_on_randomized_rows) { + std::mt19937_64 rng(20260812U); + Seen seen; + seen |= check_random<32>(rng, 32); + seen |= check_random<32>(rng, 16); + seen |= check_random<32>(rng, 29); + seen |= check_random<64>(rng, 64); + seen |= check_random<64>(rng, 32); + seen |= check_random<64>(rng, 61); + seen |= check_random<128>(rng, 128); + seen |= check_random<128>(rng, 64); + seen |= check_random<128>(rng, 125); + require_discriminating(seen); +} + +// The identities in the header, spelled out on hand-built words so a regression names the broken one. +BOOST_AUTO_TEST_CASE(codes_algebra_word_identities) { + // Slots 0..2 = 0b11, 0b10, 0b01: one paired mode, one upper-only, one lower-only. + constexpr RowCodes codes = 0b01'10'11ULL; + BOOST_TEST(row_slot_count(codes) == 3U); + const auto sums = codes_cutoff_sums(codes); + BOOST_TEST(sums.or_sum == 3U); // n + BOOST_TEST(sums.popcount_sum == 4U); // n + d, d = 1 + BOOST_TEST(sums.xor_sum == 2U); // n - d + BOOST_TEST(!codes_is_paired(codes)); + BOOST_TEST(codes_is_paired(0b11'11ULL)); + BOOST_TEST(codes_is_paired(0U)); // an empty row is vacuously paired, as dense is_paired agrees + BOOST_TEST(codes_pauli_y_count(codes) == 1U); + BOOST_TEST(codes_pair_swap(codes) == 0b10'01'11ULL); + BOOST_TEST(codes_pair_swap(codes_pair_swap(codes)) == codes); // an involution + + BOOST_TEST(codes_popcount_below(codes, 0U) == 0U); + BOOST_TEST(codes_popcount_below(codes, 1U) == 2U); + BOOST_TEST(codes_popcount_below(codes, 2U) == 3U); + BOOST_TEST(codes_popcount_below(codes, 3U) == 4U); + // Past the last slot the answer is the whole word, and the shift that would express it is undefined. + BOOST_TEST(codes_popcount_below(codes, SparseRowStore<32>::kMaxSlots) == 4U); +} diff --git a/cpp/tests/codes_product_tests.cpp b/cpp/tests/codes_product_tests.cpp new file mode 100644 index 00000000..ebbfbb52 --- /dev/null +++ b/cpp/tests/codes_product_tests.cpp @@ -0,0 +1,267 @@ +// 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. + +// The scan's per-term kernel in support form, against the dense one it must replace. What +// emit_term_products computes per term is the product M(+)G, the overlap popcount(M&G), and the basis +// rotation sign; this asserts all three agree exactly, for both algebras, on real generators from the +// fixtures and on randomized rows -- including the capacity overflow, which must be reported rather than +// silently truncating a mode list. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/Algebra.h" +#include "monoprop/algebra/CodesAlgebra.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "RandomMonomial.h" +#include "TestData.h" +#include "TestUtilities.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +// A row plus the lane storage behind it, so a test can hold several at once. The store's own rows borrow +// its arrays; these do not, which is what lets a product row be built and then compared. +struct OwnedRow { + std::vector lanes; + RowCodes codes = 0; + + explicit OwnedRow(size_t capacity) : lanes(capacity, 0) {} + + [[nodiscard]] auto view() const -> SparseRow { return SparseRow{lanes.data(), codes}; } + + template + static auto encode(const Bitset &mono, size_t capacity) -> OwnedRow { + OwnedRow row(capacity); + size_t used = 0; + for_each_mode_slot(mono, [&](size_t mode, unsigned int code) { + BOOST_REQUIRE_MESSAGE(used < capacity, "test row exceeded its capacity"); + row.lanes[used] = static_cast(mode); + row.codes |= static_cast(code) << (2 * used); + ++used; + }); + return row; + } +}; + +struct Seen { + bool cancelled_a_mode = false; // a mode present in both, cancelling to nothing + bool nonzero_overlap = false; + bool zero_overlap = false; + bool majorana_minus = false; + bool majorana_plus = false; + bool pauli_minus = false; + bool pauli_plus = false; + + auto operator|=(const Seen &o) -> Seen & { + cancelled_a_mode |= o.cancelled_a_mode; + nonzero_overlap |= o.nonzero_overlap; + zero_overlap |= o.zero_overlap; + majorana_minus |= o.majorana_minus; + majorana_plus |= o.majorana_plus; + pauli_minus |= o.pauli_minus; + pauli_plus |= o.pauli_plus; + return *this; + } +}; + +// One term against one generator, every quantity emit_term_products would produce. +template +auto check_product(const Monomial &mono, const Monomial &gen, size_t capacity, Seen &seen) + -> void { + const auto mono_row = OwnedRow::encode(mono, capacity); + const auto gen_row = OwnedRow::encode(gen, capacity); + + // The dense reference, exactly as the scan computes it. + const Monomial dense_product = mono ^ gen; + const size_t dense_overlap = mono.count_and(gen); + + std::vector out_lanes(capacity, 0); + const auto product = sparse_toggle(mono_row.view(), gen_row.view(), std::span(out_lanes)); + BOOST_REQUIRE(!product.overflowed); + + const SparseRow product_row{out_lanes.data(), product.codes}; + BOOST_TEST(product.overlap == dense_overlap); + BOOST_TEST(product.num_slots == row_slot_count(product.codes)); + BOOST_TEST((sparse_row_to_monomial<2 * NumModes>(product_row) == dense_product)); + + // The two rotation signs. Majorana's dense form goes through the per-layer interleave mask, which is + // the hot path the sparse walk replaces, so compare against that and not only against + // interleave_phase. + const auto majorana_ctx = MajoranaAlgebra::make_gen_context(gen); + const int dense_majorana = MajoranaAlgebra::rotation_sign(majorana_ctx, mono, dense_product); + const int sparse_majorana = codes_interleave_phase(mono_row.view(), gen_row.view()); + BOOST_TEST(sparse_majorana == dense_majorana); + + const auto pauli_ctx = PauliAlgebra::make_gen_context(gen); + const int dense_pauli = PauliAlgebra::rotation_sign(pauli_ctx, mono, dense_product); + const int sparse_pauli = codes_pauli_rotation_sign(mono_row.view(), gen_row.view()); + BOOST_TEST(sparse_pauli == dense_pauli); + + seen.nonzero_overlap |= dense_overlap > 0; + seen.zero_overlap |= dense_overlap == 0; + seen.majorana_minus |= dense_majorana < 0; + seen.majorana_plus |= dense_majorana > 0; + seen.pauli_minus |= dense_pauli < 0; + seen.pauli_plus |= dense_pauli > 0; + // A cancelling mode is the case a naive union would get wrong: fewer product slots than the union of + // the two inputs' modes. + size_t shared_cancelling = 0; + for_each_mode_slot(mono, [&](size_t mode, unsigned int code) { + const auto row = gen_row.view(); + for (size_t j = 0; j < row.num_slots(); ++j) { + if (row.mode(j) == mode && row.code(j) == code) { + ++shared_cancelling; + } + } + }); + seen.cancelled_a_mode |= shared_cancelling > 0; +} + +template +auto check_random_products(std::mt19937_64 &rng, size_t trials, size_t mono_slots, size_t gen_slots, Seen &seen) + -> void { + for (size_t trial = 0; trial < trials; ++trial) { + const auto mono = test_utils::random_monomial(rng, mono_slots); + const auto gen = test_utils::random_monomial(rng, gen_slots); + check_product(mono, gen, SparseRowStore::kMaxSlots, seen); + } +} + +// The fixtures' Majorana generator list against their Hamiltonian keys. +template +auto check_fixture_products(const std::string &name, Seen &seen) -> size_t { + const auto data = test_utils::load_case_data(name); + const size_t max_index = 2 * data.num_modes; + + std::vector> terms; + for (const auto &[inds, coeff] : data.hamiltonian) { + if (inds.size() <= 12) { + terms.push_back(indices_to_bitset_checked(inds, max_index)); + } + } + std::vector> gens; + for (const auto &inds : data.majoranas) { + if (inds.size() <= 12) { + gens.push_back(indices_to_bitset_checked(inds, max_index)); + } + } + BOOST_REQUIRE(!terms.empty()); + BOOST_REQUIRE(!gens.empty()); + + size_t pairs = 0; + const size_t stride = terms.size() > 40 ? (terms.size() / 40) + 1 : 1; + for (size_t i = 0; i < terms.size(); i += stride) { + for (const auto &gen : gens) { + check_product(terms[i], gen, SparseRowStore::kMaxSlots, seen); + ++pairs; + } + } + return pairs; +} + +} // namespace + +// Randomized terms against randomized generators. Generators are drawn from the same distribution and +// deliberately overlap the terms, since a disjoint generator exercises neither the overlap count nor the +// cancelling-mode branch. +BOOST_AUTO_TEST_CASE(codes_product_matches_dense_on_randomized_rows) { + std::mt19937_64 rng(20260812U); + Seen seen; + check_random_products<32>(rng, 400, 6, 4, seen); + check_random_products<64>(rng, 400, 6, 4, seen); + check_random_products<300>(rng, 400, 6, 4, seen); + BOOST_TEST(seen.cancelled_a_mode); + BOOST_TEST(seen.nonzero_overlap); + BOOST_TEST(seen.zero_overlap); + BOOST_TEST(seen.majorana_minus); + BOOST_TEST(seen.majorana_plus); + BOOST_TEST(seen.pauli_minus); + BOOST_TEST(seen.pauli_plus); +} + +// Small mode counts, so terms and generators collide constantly: nearly every product goes through the +// equal-mode branch, and many modes cancel outright. +BOOST_AUTO_TEST_CASE(codes_product_matches_dense_under_heavy_overlap) { + std::mt19937_64 rng(4242U); + Seen seen; + check_random_products<6>(rng, 2000, 6, 6, seen); + BOOST_TEST(seen.cancelled_a_mode); + BOOST_TEST(seen.nonzero_overlap); + BOOST_TEST(seen.majorana_minus); + BOOST_TEST(seen.pauli_minus); +} + +// Real generators and real terms. +BOOST_AUTO_TEST_CASE(codes_product_matches_dense_on_fixture_generators) { + Seen seen; + size_t pairs = 0; + pairs += check_fixture_products<8>("random_exact.msgpack", seen); + pairs += check_fixture_products<12>("lih_fermionic_spin_exact.msgpack", seen); + BOOST_TEST(pairs > 100U); + BOOST_TEST(seen.nonzero_overlap); + BOOST_TEST(seen.majorana_plus); + BOOST_TEST(seen.pauli_plus); +} + +// The product occupies up to the term's modes plus the generator's, which is why a scratch row is sized +// max_mode_bound() + generator locality. Past that the answer must be "overflowed", never a truncated +// mode list beside a plausible codes word -- that combination is what makes a capacity bug read as a +// speedup. +BOOST_AUTO_TEST_CASE(codes_product_reports_capacity_overflow) { + constexpr size_t kNumModes = 32; + Monomial mono; + Monomial gen; + for (const size_t mode : {0U, 1U, 2U}) { // three disjoint modes each + mono.set(2 * mode); + } + for (const size_t mode : {10U, 11U, 12U}) { + gen.set(2 * mode); + } + const auto mono_row = OwnedRow::encode(mono, 8); + const auto gen_row = OwnedRow::encode(gen, 8); + + // Six modes in the product, so five lanes is one short and six is exactly enough. + for (const size_t capacity : {1U, 3U, 5U}) { + std::vector lanes(capacity, 0); + const auto product = sparse_toggle(mono_row.view(), gen_row.view(), std::span(lanes)); + BOOST_TEST(product.overflowed); + BOOST_TEST(product.codes == 0U); + BOOST_TEST(product.num_slots == 0U); + } + std::vector lanes(6, 0); + const auto product = sparse_toggle(mono_row.view(), gen_row.view(), std::span(lanes)); + BOOST_TEST(!product.overflowed); + BOOST_TEST(product.num_slots == 6U); + BOOST_TEST(product.overlap == 0U); + + // A cancelling term needs *fewer* lanes than the union, so capacity is about the product and not + // about the inputs: gen against itself is empty. + std::vector same(1, 0); + const auto cancelled = sparse_toggle(gen_row.view(), gen_row.view(), std::span(same)); + BOOST_TEST(!cancelled.overflowed); + BOOST_TEST(cancelled.num_slots == 0U); + BOOST_TEST(cancelled.codes == 0U); + BOOST_TEST(cancelled.overlap == 3U); +} diff --git a/cpp/tests/evolution_detail_tests.cpp b/cpp/tests/evolution_detail_tests.cpp index 91f2662b..47db218d 100644 --- a/cpp/tests/evolution_detail_tests.cpp +++ b/cpp/tests/evolution_detail_tests.cpp @@ -30,6 +30,8 @@ #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/RowAccess.h" +#include "TestOperator.h" + using namespace monoprop; using monoprop::detail::CutoffContext; using monoprop::detail::MatchedEpochSet; @@ -43,15 +45,8 @@ struct RecordingSink { auto self_hit(size_t src, size_t found, int /*phase*/, double /*v_src*/) -> void { hits.emplace_back(src, found); } }; -// append_term writes a row only; find_batch needs the hash index, which insert_absent_terms populates. auto indexed_op(const std::vector> &terms) -> detail::MPOperator<8> { - detail::MPOperator<8> op; - detail::insert_absent_terms<8>( - op, - terms.size(), - [&](size_t k) -> const Monomial<8> & { return terms[k]; }, - [&](size_t k, size_t base) { assign_row<8>(*op.store, base + k, terms[k]); }); - return op; + return test_utils::indexed_operator<8>(terms); } } // namespace @@ -147,13 +142,14 @@ BOOST_AUTO_TEST_CASE(self_resolve_mark_bounded_by_combined_size) { // than past the end of epoch_, where it would be silent undefined behaviour. matched.begin_gate(op.size()); - detail::LayerBuildEngine<8, RecordingSink> eng(op, - mpi::Comm{}, - /*R_=*/1, - /*my_rank_=*/0, - matched, - combined_size, - RecordingSink{}); + detail::LayerBuildEngine<8, RecordingSink, detail::OperatorIndex<8>> eng(op, + *op.dense_rows, + mpi::Comm{}, + /*R_=*/1, + /*my_rank_=*/0, + matched, + combined_size, + RecordingSink{}); detail::query_push<8>(eng.queries_r[0], terms[1], 1); detail::query_push<8>(eng.queries_r[0], terms[5], -1); eng.src_idx_r[0] = {0, 2}; diff --git a/cpp/tests/mp_operator_tests.cpp b/cpp/tests/mp_operator_tests.cpp index 2c822c38..2dcc6efe 100644 --- a/cpp/tests/mp_operator_tests.cpp +++ b/cpp/tests/mp_operator_tests.cpp @@ -31,22 +31,15 @@ #include "monoprop/detail/operator/MPOperator.h" #include "monoprop/detail/operator/RowAccess.h" +#include "TestOperator.h" + using namespace monoprop; using cd = std::complex; namespace { -// Build an MPOperator whose store rows are also indexed (findable). append_term writes a row only; -// find() needs the hash index, which only the insert_absent_terms path populates. auto build_indexed_op(const std::vector> &terms, Basis basis = Basis::Majorana) -> detail::MPOperator<8> { - detail::MPOperator<8> op; - op.basis = basis; - detail::insert_absent_terms<8>( - op, - terms.size(), - [&](size_t k) -> const Monomial<8> & { return terms[k]; }, - [&](size_t k, size_t base) { assign_row<8>(*op.store, base + k, terms[k]); }); - return op; + return test_utils::indexed_operator<8>(terms, basis); } // Independent expected state vector: score paired rows with the basis' state phase, 0 otherwise. @@ -54,7 +47,7 @@ auto expected_state(detail::MPOperator<8> &op, Basis basis, const VecZ &initial_ const auto state_mask = initial_state_mask<8>(initial_state); VecD expected(op.size(), 0.0); for (size_t i = 0; i < op.size(); ++i) { - const auto row = materialize_row<8>(*op.store, i); + const auto row = op.with_store([i](const auto &rows) { return materialize_row<8>(rows, i); }); if (is_paired<8>(row)) { expected[i] = algebra_state_phase<8>(basis, row, state_mask); } @@ -295,19 +288,22 @@ BOOST_AUTO_TEST_CASE(mp_operator_insert_absent_terms_grows_and_indexes) { indices_to_bitset<8>({6, 7}), indices_to_bitset<8>({0, 3})}; - const size_t base = detail::insert_absent_terms<8>( - op, - fresh.size(), - [&](size_t k) -> const Monomial<8> & { return fresh[k]; }, - [&](size_t k, size_t b) { assign_row<8>(*op.store, b + k, fresh[k]); }); + const size_t base = op.with_store([&](auto &rows) { + return detail::insert_absent_terms<8>( + op, + rows, + fresh.size(), + [&](size_t k) -> const Monomial<8> & { return fresh[k]; }, + [&](size_t k, size_t b) { assign_row<8>(rows, b + k, fresh[k]); }); + }); BOOST_CHECK_EQUAL(base, 2U); BOOST_CHECK_EQUAL(op.size(), 5U); for (const auto &f : fresh) { - BOOST_CHECK(op.store->find(f).has_value()); + BOOST_CHECK(op.find(f).has_value()); } - BOOST_CHECK(op.store->find(e0).has_value()); // existing rows intact - BOOST_CHECK(op.store->find(e1).has_value()); + BOOST_CHECK(op.find(e0).has_value()); // existing rows intact + BOOST_CHECK(op.find(e1).has_value()); } BOOST_AUTO_TEST_CASE(mp_operator_append_term_after_materialization_rebuilds_inverted_index) { @@ -411,7 +407,7 @@ BOOST_AUTO_TEST_CASE(mp_operator_copy_constructor_clones_store_and_coeffs) { BOOST_CHECK(copy.state_rows_ == op.state_rows_); BOOST_CHECK(copy.state_vals_ == op.state_vals_); BOOST_CHECK(copy.materialize_state() == op.materialize_state()); - BOOST_CHECK(copy.store->find(indices_to_bitset<8>({0, 1})).has_value()); + BOOST_CHECK(copy.find(indices_to_bitset<8>({0, 1})).has_value()); // Mutating the copy must not touch the original (independent stores). copy.append_term(indices_to_bitset<8>({4, 5})); BOOST_CHECK_EQUAL(op.size(), 2U); diff --git a/cpp/tests/partition_equivalence_tests.cpp b/cpp/tests/partition_equivalence_tests.cpp index 9c8a984e..cebc2a9c 100644 --- a/cpp/tests/partition_equivalence_tests.cpp +++ b/cpp/tests/partition_equivalence_tests.cpp @@ -155,7 +155,7 @@ BOOST_AUTO_TEST_CASE(partition_raw_accessors_reject_a_facade) { sim.build_graph(data.majoranas, data.param_inds, data.gen_coeffs); BOOST_CHECK_THROW(static_cast(sim.graph()), std::runtime_error); BOOST_CHECK_THROW(static_cast(sim.mp_op()), std::runtime_error); - BOOST_CHECK_THROW(static_cast(sim.indexing()), std::runtime_error); + BOOST_CHECK_THROW(static_cast(sim.num_local_terms()), std::runtime_error); BOOST_CHECK_THROW(static_cast(sim.graph_data()), std::runtime_error); auto solo = majorana_sim(data, 1); diff --git a/cpp/tests/pauli_build_layer_tests.cpp b/cpp/tests/pauli_build_layer_tests.cpp index af069fcd..c666da57 100644 --- a/cpp/tests/pauli_build_layer_tests.cpp +++ b/cpp/tests/pauli_build_layer_tests.cpp @@ -87,7 +87,7 @@ auto dense_operator(MonomialPropagator &mp) -> std::vector { const size_t d = size_t{1} << N; std::vector m(d * d, cd(0, 0)); const auto &coeffs = mp.mp_op().get_operator(); - mp.indexing().for_each([&](const Monomial &mono, size_t idx) { + mp.for_each_term([&](const Monomial &mono, size_t idx) { if (idx >= coeffs.size()) { return; } @@ -406,7 +406,7 @@ BOOST_AUTO_TEST_CASE(pauli_build_layer_replay_fold_consumers) { const auto Gb = indices_to_bitset(slots_of_string("XII")); std::set expected; (void)mp.mp_op().get_operator(); // materialize the store size - mp.indexing().for_each([&](const Monomial &mono, size_t idx) { + mp.for_each_term([&](const Monomial &mono, size_t idx) { if (pauli_anticommutes(mono, Gb)) { expected.insert(idx); } diff --git a/cpp/tests/row_accessor_tests.cpp b/cpp/tests/row_accessor_tests.cpp index 123426a5..76c27203 100644 --- a/cpp/tests/row_accessor_tests.cpp +++ b/cpp/tests/row_accessor_tests.cpp @@ -12,7 +12,8 @@ // See the License for the specific language governing permissions and // limitations under the License. -// The dense-vector and packed OperatorIndex backends must agree through every RowAccess.h accessor. +// The dense-vector, packed OperatorIndex and sparse SparseRowStore backends must agree through every +// RowAccess.h accessor. #include @@ -22,6 +23,7 @@ #include "monoprop/algebra/MajoranaAlgebra.h" #include "monoprop/detail/operator/OperatorIndex.h" #include "monoprop/detail/operator/RowAccess.h" +#include "monoprop/detail/operator/SparseRowStore.h" using namespace monoprop; @@ -34,10 +36,13 @@ auto positions_of(const auto &backend, size_t i) -> std::vector { return out; } +// `slots` is the sparse backend's per-row mode capacity: pass one below a row's occupied-mode count to +// drive that row down the overflow path, which must stay invisible through the accessors. template -auto check_backends_agree(const std::vector> &raw_rows) -> void { +auto check_backends_agree(const std::vector> &raw_rows, size_t slots = 8) -> void { std::vector> dense; detail::OperatorIndex packed; + detail::SparseRowStore sparse(slots); for (const auto &bits : raw_rows) { Monomial m; for (size_t b : bits) { @@ -45,14 +50,19 @@ auto check_backends_agree(const std::vector> &raw_rows) -> v } dense.push_back(m); packed.push_back(m); + sparse.push_back(m); } BOOST_REQUIRE(packed.size() == dense.size()); + BOOST_REQUIRE(sparse.size() == dense.size()); for (size_t i = 0; i < dense.size(); ++i) { BOOST_TEST((materialize_row(dense, i) == materialize_row(packed, i))); + BOOST_TEST((materialize_row(dense, i) == materialize_row(sparse, i))); BOOST_TEST(row_popcount(dense, i) == row_popcount(packed, i)); + BOOST_TEST(row_popcount(dense, i) == row_popcount(sparse, i)); BOOST_TEST(row_popcount(dense, i) == materialize_row(dense, i).count()); BOOST_TEST(positions_of(dense, i) == positions_of(packed, i)); + BOOST_TEST(positions_of(dense, i) == positions_of(sparse, i)); } } @@ -66,25 +76,45 @@ BOOST_AUTO_TEST_CASE(row_accessor_backends_agree_multi_word) { check_backends_agree<96>({{0, 64, 191}, {5, 63, 64, 65}, {}, {128, 190}}); } +// Four occupied modes against a two-slot capacity: the first two rows spill, the empty row and the +// one-mode row do not, so the same store serves both kinds. +BOOST_AUTO_TEST_CASE(row_accessor_backends_agree_sparse_overflow) { + check_backends_agree<32>({{0, 3, 5, 8, 20, 21}, {1, 2, 40, 41, 62, 63}, {}, {10, 11}}, 2); +} + BOOST_AUTO_TEST_CASE(row_accessor_assign_row_overwrites) { constexpr size_t N = 32; std::vector> dense; detail::OperatorIndex packed; + // Two slots, and the original occupies three modes: the row starts spilled and the overwrite must + // pull it back inline rather than leaving the stale side-map entry to shadow it. + detail::SparseRowStore sparse(2); Monomial original; original.set(1); original.set(2); + original.set(40); + original.set(41); + original.set(60); dense.push_back(original); packed.push_back(original); + sparse.push_back(original); + BOOST_TEST(sparse.spilled(0)); + // Three set bits over two modes, so it fits the sparse store's two slots. Monomial replacement; replacement.set(10); - replacement.set(20); + replacement.set(11); replacement.set(30); assign_row(dense, 0, replacement); assign_row(packed, 0, replacement); + assign_row(sparse, 0, replacement); BOOST_TEST((materialize_row(dense, 0) == replacement)); BOOST_TEST((materialize_row(packed, 0) == replacement)); + BOOST_TEST((materialize_row(sparse, 0) == replacement)); BOOST_TEST(row_popcount(packed, 0) == 3U); + BOOST_TEST(row_popcount(sparse, 0) == 3U); + BOOST_TEST(!sparse.spilled(0)); BOOST_TEST(positions_of(dense, 0) == positions_of(packed, 0)); + BOOST_TEST(positions_of(dense, 0) == positions_of(sparse, 0)); } diff --git a/cpp/tests/row_store_selection_tests.cpp b/cpp/tests/row_store_selection_tests.cpp new file mode 100644 index 00000000..2767c181 --- /dev/null +++ b/cpp/tests/row_store_selection_tests.cpp @@ -0,0 +1,92 @@ +// 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. + +// Which row backend a propagator ends up on, and that monoprop_ROW_STORE is wired to it. +// +// This is the guard against a false negative in the rest of the suite: every case runs a second time +// under monoprop_ROW_STORE=sparse (see cpp/tests/boostAddTests.cmake), and every fixture is far below +// the automatic crossover -- so if the variable reached nothing, those extra passes would be the dense +// backend passing twice and nobody would notice. + +#include + +#include + +#include "monoprop/MonomialPropagator.h" +#include "monoprop/detail/EnvConfig.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "TestData.h" +#include "TestUtilities.h" + +using namespace monoprop; + +namespace { + +constexpr size_t kNumModes = 8; + +auto small_propagator() -> MonomialPropagator { + const auto data = test_utils::load_case_data("random_exact.msgpack"); + return test_utils::build_simulator(data); +} + +} // namespace + +// The propagator's backend must be what the environment asked for -- and, unset, what the crossover +// says. Both arms are live: the default ctest variant takes the first, the sparse-rows variant the +// second, so this case is the one that fails if the variable is ignored. +BOOST_AUTO_TEST_CASE(row_store_selection_follows_the_environment) { + const auto propagator = small_propagator(); + // A fixture-sized system: below every shipped crossover, so `auto` must pick dense here. If this + // ever fails, the suite's sparse coverage has stopped being a second configuration. + BOOST_REQUIRE(!monoprop::detail::SparseRowStore::preferred_for_modes(kNumModes)); + + // Dereferenced: an unrecognized value is nullopt, and the propagator above would have thrown on it + // before this line -- see the rejects-an-unrecognized-value case below. + BOOST_REQUIRE(config::get().row_store.has_value()); + switch (*config::get().row_store) { + case config::RowStore::Sparse: + BOOST_TEST(propagator.rows_are_sparse()); + break; + case config::RowStore::Dense: + case config::RowStore::Auto: + BOOST_TEST(!propagator.rows_are_sparse()); + break; + } +} + +// The automatic rule, independent of any propagator: the crossover is a whole 32-mode block, so no +// storage width can straddle it. +BOOST_AUTO_TEST_CASE(row_store_auto_crossover_is_a_whole_storage_block) { + constexpr size_t kMin = monoprop::detail::SparseRowStore::kMinModes; + BOOST_TEST(kMin % 32 == 0U); + BOOST_TEST(!monoprop::detail::SparseRowStore::preferred_for_modes(kMin - 1)); + BOOST_TEST(monoprop::detail::SparseRowStore::preferred_for_modes(kMin)); + BOOST_TEST(monoprop::detail::SparseRowStore::preferred_for_modes(kMin + 32)); +} + +// An unrecognized value must be rejected, not silently treated as `auto`: the whole reason to set the +// variable is to know which backend ran, and a typo that fell back would read as a passing run of a +// configuration that never happened. Parsed here rather than through config::get(), which caches the +// process environment once and so cannot be re-read per case. +BOOST_AUTO_TEST_CASE(row_store_env_parses_only_the_three_values) { + BOOST_TEST((config::detail::parse_row_store(nullptr) == config::RowStore::Auto)); + BOOST_TEST((config::detail::parse_row_store("") == config::RowStore::Auto)); + BOOST_TEST((config::detail::parse_row_store("auto") == config::RowStore::Auto)); + BOOST_TEST((config::detail::parse_row_store("dense") == config::RowStore::Dense)); + BOOST_TEST((config::detail::parse_row_store("sparse") == config::RowStore::Sparse)); + for (const char *bad : {"Sparse", "SPARSE", "spars", "sparse ", "1", "on", "packed"}) { + BOOST_TEST(!config::detail::parse_row_store(bad).has_value()); + } +} diff --git a/cpp/tests/simulator_copy_tests.cpp b/cpp/tests/simulator_copy_tests.cpp index be79c160..4ba939b2 100644 --- a/cpp/tests/simulator_copy_tests.cpp +++ b/cpp/tests/simulator_copy_tests.cpp @@ -22,7 +22,7 @@ // Copy-constructing a simulator must produce a fully independent deep copy -- the mechanism behind // Python __deepcopy__. The operator store is non-copyable, so the copy rebuilds it via clone() and -// find()/indexing() have to work on the copy's own rows. The MPI communicator handle is shared. +// find()/for_each_term() have to work on the copy's own rows. The MPI communicator handle is shared. using namespace test_utils; using namespace monoprop; @@ -97,11 +97,10 @@ BOOST_FIXTURE_TEST_CASE(copy_constructed_simulator_index_valid, ExampleDataFix) auto copy = sim; - const auto &idx = copy.indexing(); - BOOST_TEST(idx.size() == sim.indexing().size()); + BOOST_TEST(copy.num_local_terms() == sim.num_local_terms()); bool all_found = true; - idx.for_each([&](const auto &mono, size_t i) { - const auto f = idx.find(mono); + copy.for_each_term([&](const auto &mono, size_t i) { + const auto f = copy.mp_op().find(mono); if (!f || *f != i) { all_found = false; } diff --git a/cpp/tests/sparse_index_tests.cpp b/cpp/tests/sparse_index_tests.cpp new file mode 100644 index 00000000..beaeaaaf --- /dev/null +++ b/cpp/tests/sparse_index_tests.cpp @@ -0,0 +1,314 @@ +// 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. + +// SparseRowStore's keyless index must answer lookups exactly as OperatorIndex's does -- same hits, same +// misses, same insert-or-no-op -- because a propagator swaps one for the other. The hash *value* differs +// by design (that is the documented re-baseline); nothing else may. + +#include + +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "RandomMonomial.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +// Distinct rows only: bulk_insert and the emplace-is-idempotent check both require it, and a duplicate +// would make "found at index i" ambiguous. +template +struct RowSet { + std::vector> rows; + std::set> seen; + + auto add(const Monomial &mono) -> bool { + if (!seen.insert(positions(mono)).second) { + return false; + } + rows.push_back(mono); + return true; + } + [[nodiscard]] auto contains(const Monomial &mono) const -> bool { + return seen.contains(positions(mono)); + } + + static auto positions(const Monomial &mono) -> std::vector { + std::vector out; + for (size_t b = mono.find_first(); b < mono.size(); b = mono.find_next(b)) { + out.push_back(b); + } + return out; + } +}; + +// A spilled row has no codes word, so it can only be hashed by walking the dense monomial. That walk and +// the row walk are two producers of one hash, and if they ever disagreed a spilled row would become +// unfindable -- silently, and only for the small share of rows that spill. +template +auto check_hash_agrees(std::mt19937_64 &rng) -> size_t { + using Store = SparseRowStore; + Store store(Store::kMaxSlots); + std::vector> rows; + for (size_t t = 0; t < 200; ++t) { + rows.push_back(test_utils::random_monomial(rng, Store::kMaxSlots)); + store.push_back(rows.back()); + } + for (size_t i = 0; i < rows.size(); ++i) { + BOOST_REQUIRE(!store.spilled(i)); + BOOST_TEST(sparse_row_hash(store.view(i)) == sparse_row_hash(rows[i])); + } + return rows.size(); +} + +// The lookup contract, against OperatorIndex as the oracle. +template +auto check_lookups_match_packed(std::mt19937_64 &rng, size_t slots, size_t &spilled_rows, size_t &misses) -> void { + constexpr size_t kNumBits = 2 * NumModes; + SparseRowStore sparse(slots); + OperatorIndex packed; + RowSet set; + for (size_t t = 0; t < 300; ++t) { + const auto mono = test_utils::random_monomial(rng, 12); + if (!set.add(mono)) { + continue; + } + const size_t i = sparse.size(); + sparse.push_back(mono); + packed.push_back(mono); + sparse.emplace(mono, i); + packed.emplace(mono, i); + spilled_rows += sparse.spilled(i) ? 1 : 0; + } + BOOST_REQUIRE(sparse.size() == set.rows.size()); + + for (size_t i = 0; i < set.rows.size(); ++i) { + const auto by_mono = sparse.find(set.rows[i]); + BOOST_REQUIRE(by_mono.has_value()); + BOOST_TEST(*by_mono == i); + BOOST_TEST(*by_mono == packed.find(set.rows[i]).value()); + if (!sparse.spilled(i)) { + const auto by_row = sparse.find(sparse.view(i)); + BOOST_REQUIRE(by_row.has_value()); + BOOST_TEST(*by_row == i); + } + } + + // Absent keys must miss, not land on a hash neighbour. + for (size_t t = 0; t < 300; ++t) { + Monomial mono; + for (size_t k = 0; k < 1 + (rng() % 6); ++k) { + mono.set(rng() % kNumBits); + } + if (set.contains(mono)) { + continue; + } + BOOST_TEST(!sparse.find(mono).has_value()); + BOOST_TEST(!packed.find(mono).has_value()); + ++misses; + } + + // emplace is insert-or-no-op, so replaying every key must add no slot. + const size_t indexed = sparse.indexed_count(); + for (size_t i = 0; i < set.rows.size(); ++i) { + sparse.emplace(set.rows[i], i); + } + BOOST_TEST(sparse.indexed_count() == indexed); + + const auto copy = sparse.clone(); + BOOST_TEST(copy->indexed_count() == indexed); + for (size_t i = 0; i < set.rows.size(); ++i) { + const auto found = copy->find(set.rows[i]); + BOOST_REQUIRE(found.has_value()); + BOOST_TEST(*found == i); + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_index_hash_agrees_between_row_and_monomial) { + std::mt19937_64 rng(20260812U); + size_t compared = 0; + compared += check_hash_agrees<32>(rng); + compared += check_hash_agrees<64>(rng); + compared += check_hash_agrees<300>(rng); + BOOST_TEST(compared == 600U); +} + +// The row capacity is a tuning parameter. If it leaked into the hash, changing it would move probe order +// and MPI owner routing, so two stores tuned differently must agree on every row they can both hold. +BOOST_AUTO_TEST_CASE(sparse_index_hash_is_independent_of_row_capacity) { + std::mt19937_64 rng(4321U); + constexpr size_t kNumModes = 32; + SparseRowStore narrow(8); + SparseRowStore wide(SparseRowStore::kMaxSlots); + size_t compared = 0; + for (size_t t = 0; t < 300; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + narrow.push_back(mono); + wide.push_back(mono); + const size_t i = narrow.size() - 1; + if (narrow.spilled(i)) { + continue; + } + BOOST_TEST(sparse_row_hash(narrow.view(i)) == sparse_row_hash(wide.view(i))); + ++compared; + } + BOOST_TEST(compared > 200U); +} + +// Three widths and three row capacities, so that spilled and inline rows both occur. +BOOST_AUTO_TEST_CASE(sparse_index_lookups_match_the_packed_backend) { + std::mt19937_64 rng(99U); + size_t spilled_rows = 0; + size_t misses = 0; + for (const size_t slots : {4U, 8U, 32U}) { + check_lookups_match_packed<32>(rng, slots, spilled_rows, misses); + check_lookups_match_packed<64>(rng, slots, spilled_rows, misses); + check_lookups_match_packed<300>(rng, slots, spilled_rows, misses); + } + // Both row kinds have to have occurred, or the spill path above was never taken. + BOOST_TEST(spilled_rows > 0U); + BOOST_TEST(misses > 100U); +} + +// find_batch is a pipelined re-implementation of find, so it needs its own equality check -- against +// find, for both key forms, over a query list mixing hits and misses. +BOOST_AUTO_TEST_CASE(sparse_index_find_batch_matches_find) { + std::mt19937_64 rng(555U); + constexpr size_t kNumModes = 64; + constexpr size_t kNumBits = 2 * kNumModes; + SparseRowStore store(8); + RowSet set; + for (size_t t = 0; t < 400; ++t) { + const auto mono = test_utils::random_monomial(rng, 10); + if (!set.add(mono)) { + continue; + } + const size_t i = store.size(); + store.push_back(mono); + store.emplace(mono, i); + } + + // More than one group of 16, and interleaved absentees so a group holds both. + std::vector> queries; + for (size_t i = 0; i < set.rows.size(); ++i) { + queries.push_back(set.rows[i]); + if (i % 3 == 0) { + Monomial absent; + absent.set(rng() % kNumBits); + if (!set.contains(absent)) { + queries.push_back(absent); + } + } + } + BOOST_REQUIRE(queries.size() > 16U); + + std::vector batched(queries.size()); + store.find_batch(queries.data(), queries.size(), batched.data()); + for (size_t j = 0; j < queries.size(); ++j) { + const auto scalar = store.find(queries[j]); + BOOST_TEST(batched[j] == (scalar ? *scalar : SparseRowStore::kNotFound)); + } + + // The row-key form, which is what the scan will hand it. + std::vector row_queries; + std::vector expected; + for (size_t i = 0; i < store.size(); ++i) { + if (!store.spilled(i)) { + row_queries.push_back(store.view(i)); + expected.push_back(i); + } + } + BOOST_REQUIRE(row_queries.size() > 16U); + std::vector row_batched(row_queries.size()); + store.find_batch(row_queries.data(), row_queries.size(), row_batched.data()); + BOOST_TEST(row_batched == expected, boost::test_tools::per_element()); +} + +// bulk_insert skips the duplicate probe, so it is only correct on provably distinct keys; what it must +// still produce is an index that finds every one of them. +BOOST_AUTO_TEST_CASE(sparse_index_bulk_insert_indexes_every_row) { + std::mt19937_64 rng(777U); + constexpr size_t kNumModes = 32; + SparseRowStore store(8); + RowSet set; + for (size_t t = 0; t < 200; ++t) { + const auto mono = test_utils::random_monomial(rng, 6); + if (set.add(mono)) { + store.push_back(mono); + } + } + store.bulk_insert(store.size(), 0, [&](size_t k) { return set.rows[k]; }); + BOOST_TEST(store.indexed_count() == store.size()); + for (size_t i = 0; i < set.rows.size(); ++i) { + const auto found = store.find(set.rows[i]); + BOOST_REQUIRE(found.has_value()); + BOOST_TEST(*found == i); + } +} + +// The sparse hash is a different function from the dense one, which is exactly why the store swap needs +// a re-baseline. Pinned so that "the results moved" is never a surprise, and so a future change that +// accidentally reunified them would be noticed. +BOOST_AUTO_TEST_CASE(sparse_index_hash_differs_from_the_dense_hash) { + std::mt19937_64 rng(8888U); + constexpr size_t kNumModes = 32; + size_t differing = 0; + size_t total = 0; + for (size_t t = 0; t < 200; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + if (mono.count() == 0) { + continue; + } + differing += sparse_row_hash(mono) != monomial_hash(mono) ? 1 : 0; + ++total; + } + BOOST_REQUIRE(total > 100U); + BOOST_TEST(differing == total); +} + +// A pre-filter that collided often would still be correct but would degrade every probe into a lane +// compare. Over distinct rows the 32-bit folds should be near-injective. +BOOST_AUTO_TEST_CASE(sparse_index_hash_folds_are_near_injective) { + std::mt19937_64 rng(31337U); + constexpr size_t kNumModes = 64; + SparseRowStore store(SparseRowStore::kMaxSlots); + RowSet set; + for (size_t t = 0; t < 4000; ++t) { + const auto mono = test_utils::random_monomial(rng, 10); + if (set.add(mono)) { + store.push_back(mono); + } + } + std::set folds; + for (const auto &mono : set.rows) { + const size_t full = sparse_row_hash(mono); + folds.insert(static_cast(full ^ (static_cast(full) >> 32))); + } + BOOST_REQUIRE(set.rows.size() > 3000U); + // Birthday-bound expectation for n draws from 2^32 is ~n^2/2^33 collisions, i.e. under 1 for n=4000; + // allowing 4 keeps this from being flaky while still failing on a hash that structurally collides. + BOOST_TEST(set.rows.size() - folds.size() <= 4U); +} diff --git a/cpp/tests/sparse_row_store_tests.cpp b/cpp/tests/sparse_row_store_tests.cpp new file mode 100644 index 00000000..db400120 --- /dev/null +++ b/cpp/tests/sparse_row_store_tests.cpp @@ -0,0 +1,532 @@ +// 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. + +// SparseRowStore's own invariants. The three-way agreement with the dense and packed backends through +// the RowAccess.h accessors lives in row_accessor_tests.cpp; what is checked here is the part that has +// no counterpart in the other backends -- the codes word, and the row sizing that feeds it. + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/algebra/AlgebraCommon.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "RandomMonomial.h" + +using namespace monoprop; +using namespace monoprop::detail; + +// Interchangeable with OperatorIndex means the same ownership rules, so the store cannot be silently +// copied out of MPOperator's unique_ptr. +static_assert(!std::is_move_constructible_v>, "SparseRowStore must remain non-movable"); +static_assert(!std::is_copy_constructible_v>, "SparseRowStore must remain non-copyable"); +// A mode lane addresses at most kMaxModes modes, and the widest store the suite builds must fit. +static_assert(SparseRowStore::kMaxModes>::kMaxModes > 0, "kMaxModes must be instantiable"); + +namespace { + +// The support-form identities, spelled out here against the dense cutoff_sums. +auto sums_from_codes(RowCodes codes) -> CutoffSums { + const auto n = static_cast(std::popcount(row_occupied_bits(codes))); + const auto d = static_cast(std::popcount(row_paired_bits(codes))); + return {n - d, n + d, n}; +} + +// Owning lanes plus codes, so a test can hold a row the way the scan's scratch does. The lanes come out +// ascending and contiguous from slot 0, which is what the representation requires of every producer. +struct RowBuffer { + std::vector lanes; + RowCodes codes = 0; + + [[nodiscard]] auto view() const -> SparseRow { return SparseRow{lanes.data(), codes}; } +}; + +// Only for rows that fit one codes word (<= kRowMaxSlots slots) -- the shape a SparseRow can hold at all. +template +auto row_of(const Bitset &mono) -> RowBuffer { + RowBuffer out; + for_each_mode_slot(mono, [&](size_t mode, unsigned int code) { + out.codes |= static_cast(code) << (2 * out.lanes.size()); + out.lanes.push_back(static_cast(mode)); + }); + return out; +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_row_store_codes_encode_slot_pairs) { + constexpr size_t kNumModes = 32; + SparseRowStore store(8); + + // Modes 1 (both positions), 4 (upper only) and 9 (lower only), so the codes word must read + // 0b11, 0b10, 0b01 from slot 0 up. + Monomial mono; + mono.set(2); + mono.set(3); + mono.set(9); + mono.set(18); + store.push_back(mono); + + BOOST_TEST(!store.spilled(0)); + // Slot 0 = mode 1 (0b11) in bits 0-1, slot 1 = mode 4 (0b10) in bits 2-3, slot 2 = mode 9 (0b01) + // in bits 4-5. + BOOST_TEST(store.codes(0) == 0b01'10'11ULL); + BOOST_TEST(store.slot_count(0) == 3U); + BOOST_TEST(store.popcount(0) == 4U); + + std::vector modes; + std::vector codes; + store.for_each_slot(0, [&](size_t mode, unsigned int code) { + modes.push_back(mode); + codes.push_back(code); + }); + BOOST_REQUIRE(modes.size() == 3U); + BOOST_TEST(modes == (std::vector{1U, 4U, 9U}), boost::test_tools::per_element()); + BOOST_TEST(codes == (std::vector{0b11U, 0b10U, 0b01U}), boost::test_tools::per_element()); +} + +BOOST_AUTO_TEST_CASE(sparse_row_store_empty_row_has_empty_codes) { + SparseRowStore<32> store(4); + store.push_back(Monomial<32>{}); + BOOST_TEST(!store.spilled(0)); + BOOST_TEST(store.codes(0) == 0U); + BOOST_TEST(store.slot_count(0) == 0U); + BOOST_TEST(store.popcount(0) == 0U); + BOOST_TEST(store.row(0) == Monomial<32>{}); +} + +namespace { + +// The identities the whole support form rests on: or_sum = n, popcount_sum = n + d, xor_sum = n - d, +// against the dense cutoff_sums over the full storage window. Randomized rather than enumerated +// because what can break them is a particular occupancy pattern, not a particular width. +template +auto check_codes_reproduce_cutoff_sums(std::mt19937_64 &rng) -> void { + // Capacity above any row built below, so nothing spills and every row exercises the codes path. + SparseRowStore store(SparseRowStore::kMaxSlots); + for (size_t trial = 0; trial < 200; ++trial) { + const auto mono = test_utils::random_monomial(rng, SparseRowStore::kMaxSlots); + store.push_back(mono); + const size_t i = store.size() - 1; + BOOST_REQUIRE(!store.spilled(i)); + + const auto dense = cutoff_sums(mono, NumModes); + const auto sparse = sums_from_codes(store.codes(i)); + BOOST_TEST(sparse.or_sum == dense.or_sum); + BOOST_TEST(sparse.popcount_sum == dense.popcount_sum); + BOOST_TEST(sparse.xor_sum == dense.xor_sum); + BOOST_TEST(store.slot_count(i) == dense.or_sum); + BOOST_TEST(store.popcount(i) == dense.popcount_sum); + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_row_store_codes_reproduce_cutoff_sums) { + std::mt19937_64 rng(20260812U); + check_codes_reproduce_cutoff_sums<32>(rng); + check_codes_reproduce_cutoff_sums<64>(rng); + check_codes_reproduce_cutoff_sums<128>(rng); + check_codes_reproduce_cutoff_sums<512>(rng); +} + +// Spilled rows have no codes word, so the two measures must still come off the dense monomial. +BOOST_AUTO_TEST_CASE(sparse_row_store_spilled_rows_report_the_same_measures) { + constexpr size_t kNumModes = 64; + SparseRowStore store(2); + Monomial mono; + for (const size_t b : {0U, 1U, 4U, 20U, 21U, 99U}) { // modes 0 (paired), 2, 10 (paired), 49 + mono.set(b); + } + store.push_back(mono); + + BOOST_TEST(store.spilled(0)); + BOOST_TEST(store.row(0) == mono); + const auto dense = cutoff_sums(mono, kNumModes); + BOOST_TEST(store.slot_count(0) == dense.or_sum); + BOOST_TEST(store.popcount(0) == dense.popcount_sum); +} + +BOOST_AUTO_TEST_CASE(sparse_row_store_clone_preserves_rows_and_spills) { + constexpr size_t kNumModes = 32; + SparseRowStore store(2); + Monomial inline_row; + inline_row.set(4); + inline_row.set(5); + Monomial spilled_row; + for (const size_t b : {0U, 6U, 10U, 30U}) { + spilled_row.set(b); + } + store.push_back(inline_row); + store.push_back(spilled_row); + + const auto copy = store.clone(); + BOOST_REQUIRE(copy->size() == 2U); + BOOST_TEST(copy->num_bits() == 2 * kNumModes); + BOOST_TEST(copy->slots_per_row() == 2U); + BOOST_TEST(!copy->spilled(0)); + BOOST_TEST(copy->codes(0) == store.codes(0)); + BOOST_TEST(copy->row(0) == inline_row); + BOOST_TEST(copy->spilled(1)); + BOOST_TEST(copy->row(1) == spilled_row); +} + +// resized() is the row-width migration: every row keeps its index and its content, whichever way the +// width moved -- including a row that crosses the overflow boundary, since set() re-decides that per +// row rather than trusting the old classification. +BOOST_AUTO_TEST_CASE(sparse_row_store_resized_preserves_rows_when_widening) { + constexpr size_t kNumModes = 32; + SparseRowStore store(2); // width 2: the 3-slot row below starts spilled + Monomial inline_row; + inline_row.set(4); + inline_row.set(5); + Monomial spilled_row; + for (const size_t b : {0U, 6U, 10U, 30U}) { + spilled_row.set(b); + } + store.push_back(inline_row); + store.push_back(spilled_row); + + const auto wide = store.resized(4); // now fits inline + BOOST_REQUIRE(wide->size() == 2U); + BOOST_TEST(wide->slots_per_row() == 4U); + BOOST_TEST(!wide->spilled(0)); + BOOST_TEST(wide->row(0) == inline_row); + BOOST_TEST(!wide->spilled(1)); + BOOST_TEST(wide->row(1) == spilled_row); + + // Independent of the source: mutating store after the fact must not reach wide. + store.set(0, Monomial{}); + BOOST_TEST(wide->row(0) == inline_row); +} + +BOOST_AUTO_TEST_CASE(sparse_row_store_resized_preserves_rows_when_narrowing) { + constexpr size_t kNumModes = 32; + SparseRowStore store(4); // width 4: both rows below fit inline + Monomial a; + a.set(4); + a.set(5); + Monomial b; + for (const size_t bit : {0U, 6U, 10U, 30U}) { + b.set(bit); + } + store.push_back(a); + store.push_back(b); + BOOST_REQUIRE(!store.spilled(1)); + + const auto narrow = store.resized(2); // row 1 must now spill + BOOST_REQUIRE(narrow->size() == 2U); + BOOST_TEST(narrow->slots_per_row() == 2U); + BOOST_TEST(!narrow->spilled(0)); + BOOST_TEST(narrow->row(0) == a); + BOOST_TEST(narrow->spilled(1)); + BOOST_TEST(narrow->row(1) == b); +} + +// K comes from the cutoff in modes, and is the same number for both cutoff kinds -- halving the slot +// bound for a support cutoff would truncate rows a length cutoff of the same size admits. +BOOST_AUTO_TEST_CASE(sparse_row_store_slots_come_from_the_cutoff_in_modes) { + constexpr size_t kNumModes = 32; + using Store = SparseRowStore; + + CutoffFn length = LengthCutoff{.cutoff = 6U}; + CutoffFn support = SupportCutoff{.cutoff = 6U}; + BOOST_TEST(CutoffEvaluator(length).max_slot_bound().value() == 6U); + BOOST_TEST(CutoffEvaluator(support).max_slot_bound().value() == 12U); + + BOOST_TEST(Store::slots_for_bound(6U) == 6U); + // A bound past one codes word clamps rather than throwing: the excess rows spill. + BOOST_TEST(Store::slots_for_bound(100U) == Store::kMaxSlots); + BOOST_TEST(Store::slots_for_bound(0U) == 1U); +} + +// The switch rule itself: a build-time constant, because the crossover follows the target ISA rather +// than anything known at run time. The values are the measured crossovers, so what this pins is that +// the CMake default reached the compiler at all -- a missing definition would silently fall back. +BOOST_AUTO_TEST_CASE(sparse_row_store_preference_threshold_matches_the_build) { + using Store = SparseRowStore<32>; + static_assert(Store::kMinModes > 0, "the sparse crossover must be a positive mode count"); + BOOST_TEST(!Store::preferred_for_modes(Store::kMinModes - 1)); + BOOST_TEST(Store::preferred_for_modes(Store::kMinModes)); + BOOST_TEST(Store::preferred_for_modes(Store::kMinModes + 1)); + // 32 modes is where sparse measured 1.9x behind dense even on baseline x86-64; no build should be + // switching there. + BOOST_TEST(!Store::preferred_for_modes(32U)); +} + +namespace { + +// The row form of set() against the dense one, which is the only definition of what it must produce. +// Both fill the same row of two stores built alike; every observable of the two must agree, including +// the hash the table keys on -- a row written one way has to be findable by a key written the other. +template +auto check_row_form_set_matches_dense(std::mt19937_64 &rng) -> void { + using Store = SparseRowStore; + using Key = SparseRowKey<2 * NumModes>; + Store dense_written(Store::kMaxSlots); + Store row_written(Store::kMaxSlots); + for (size_t trial = 0; trial < 200; ++trial) { + const auto mono = test_utils::random_monomial(rng, Store::kMaxSlots); + const RowBuffer row = row_of(mono); + + dense_written.push_back(mono); + const size_t i = row_written.grow_rows_geometric(1); + row_written.set(i, row.view()); + + BOOST_REQUIRE(i == dense_written.size() - 1); + BOOST_TEST(row_written.spilled(i) == dense_written.spilled(i)); + BOOST_TEST(row_written.codes(i) == dense_written.codes(i)); + BOOST_TEST(row_written.slot_count(i) == dense_written.slot_count(i)); + BOOST_TEST(row_written.popcount(i) == dense_written.popcount(i)); + BOOST_TEST(row_written.row(i) == mono); + // The three key forms are interchangeable only if they hash alike; the table's probe order + // (and so MPI owner routing) is downstream of this. + BOOST_TEST(sparse_row_hash(row.view()) == sparse_row_hash(mono)); + BOOST_TEST(sparse_row_hash(Key{.row = row.view()}) == sparse_row_hash(mono)); + BOOST_TEST(sparse_row_hash(Key{.spilled = &mono}) == sparse_row_hash(mono)); + } + // Written rows are findable by either form, whichever way they went in. + for (size_t i = 0; i < row_written.size(); ++i) { + const auto mono = row_written.row(i); + row_written.emplace(row_of(mono).view(), i); + } + for (size_t i = 0; i < row_written.size(); ++i) { + const auto mono = row_written.row(i); + const RowBuffer row = row_of(mono); + const auto by_mono = row_written.find(mono); + const auto by_row = row_written.find(row.view()); + const auto by_key = row_written.find(Key{.row = row.view()}); + BOOST_REQUIRE(by_mono.has_value()); + BOOST_TEST((by_row == by_mono)); + BOOST_TEST((by_key == by_mono)); + // Distinct monomials may repeat across trials, so the found row must equal this one rather + // than be this index. + BOOST_TEST(row_written.row(*by_mono) == mono); + } +} + +} // namespace + +BOOST_AUTO_TEST_CASE(sparse_row_store_row_form_set_matches_the_dense_one) { + std::mt19937_64 rng(20260813U); + check_row_form_set_matches_dense<32>(rng); + check_row_form_set_matches_dense<96>(rng); + check_row_form_set_matches_dense<512>(rng); +} + +// A row wider than the store's capacity has to spill, exactly as the dense set() spills it. The capacity +// is sized from the cutoff and a fully paired term escapes the cutoff, so this arm is reachable by +// construction and cannot be sized away. +BOOST_AUTO_TEST_CASE(sparse_row_store_row_form_set_spills_a_row_past_the_capacity) { + constexpr size_t kNumModes = 64; + using Key = SparseRowKey<2 * kNumModes>; + SparseRowStore store(2); + Monomial mono; + for (const size_t b : {0U, 1U, 4U, 20U, 21U, 99U}) { // modes 0 (paired), 2, 10 (paired), 49 + mono.set(b); + } + const RowBuffer row = row_of(mono); + BOOST_REQUIRE(row.lanes.size() == 4U); + + const size_t i = store.grow_rows_geometric(1); + store.set(i, row.view()); + BOOST_TEST(store.spilled(i)); + BOOST_TEST(store.row(i) == mono); + BOOST_TEST(store.slot_count(i) == 4U); + BOOST_TEST(store.popcount(i) == 6U); + + // A spilled row is found by the dense key; the row key finds it through either shape it carries. + store.emplace(mono, i); + BOOST_TEST((store.find(mono) == std::optional(i))); + BOOST_TEST((store.find(Key{.spilled = &mono}) == std::optional(i))); + BOOST_TEST((store.find(Key{.row = row.view()}) == std::optional(i))); +} + +// Rows are overwritten in place by the miss inserts, so a row must not inherit the previous occupant's +// spill -- in either direction. +BOOST_AUTO_TEST_CASE(sparse_row_store_row_form_set_clears_a_stale_spill) { + constexpr size_t kNumModes = 32; + SparseRowStore store(2); + Monomial wide; + for (const size_t b : {0U, 6U, 10U, 30U}) { + wide.set(b); + } + Monomial narrow; + narrow.set(4); + narrow.set(5); + + const size_t i = store.grow_rows_geometric(1); + store.set(i, row_of(wide).view()); + BOOST_REQUIRE(store.spilled(i)); + + store.set(i, row_of(narrow).view()); + BOOST_TEST(!store.spilled(i)); + BOOST_TEST(store.row(i) == narrow); + BOOST_TEST(store.codes(i) == 0b11ULL); // one slot, mode 2, both positions + + store.set(i, row_of(wide).view()); + BOOST_TEST(store.spilled(i)); + BOOST_TEST(store.row(i) == wide); + + // The empty row is the case the lane padding exists for: it writes no lane of its own, so without the + // pad the previous occupant's overflow marker survives in lane 0 and the row reads as spilled while + // its side-map entry is gone. + store.set(i, row_of(Monomial{}).view()); + BOOST_TEST(!store.spilled(i)); + BOOST_TEST(store.codes(i) == 0U); + BOOST_TEST(store.slot_count(i) == 0U); + BOOST_TEST(store.row(i) == Monomial{}); +} + +// A batch of keys is homogeneous by type but not by shape: the spilled ones are what a query record +// escapes to. find_batch must resolve both, in one pass, against a store holding both kinds of row. +BOOST_AUTO_TEST_CASE(sparse_row_key_batch_resolves_both_shapes) { + constexpr size_t kNumModes = 64; + using Store = SparseRowStore; + using Key = SparseRowKey<2 * kNumModes>; + Store store(3); + std::mt19937_64 rng(20260814U); + + std::vector> monos; + for (size_t k = 0; k < 64; ++k) { + // Up to 6 slots against a capacity of 3, so roughly half the rows spill. + auto mono = test_utils::random_monomial(rng, 6); + if (std::find(monos.begin(), monos.end(), mono) != monos.end()) { + continue; + } + store.push_back(mono); + store.emplace(mono, store.size() - 1); + monos.push_back(mono); + } + BOOST_REQUIRE(monos.size() > 8U); + + std::vector rows; + std::vector keys; + rows.reserve(monos.size()); + keys.reserve(monos.size()); + for (const auto &mono : monos) { + rows.push_back(row_of(mono)); + } + size_t spilled_keys = 0; + for (size_t k = 0; k < monos.size(); ++k) { + // The shape a record would have carried: within the capacity it stays a row, past it the record + // escapes to the dense monomial. + if (rows[k].lanes.size() > store.slots_per_row()) { + keys.push_back(Key{.spilled = &monos[k]}); + ++spilled_keys; + } + else { + keys.push_back(Key{.row = rows[k].view()}); + } + } + BOOST_REQUIRE(spilled_keys > 0U); + BOOST_REQUIRE(spilled_keys < monos.size()); + + std::vector found(keys.size(), Store::kNotFound); + store.find_batch(keys.data(), keys.size(), found.data()); + for (size_t k = 0; k < keys.size(); ++k) { + BOOST_REQUIRE(found[k] < store.size()); + BOOST_TEST(store.row(found[k]) == monos[k]); + } + + // An absent key must miss through either shape. + Monomial absent; + absent.set(6); + absent.set(7); + absent.set(120); + while (std::find(monos.begin(), monos.end(), absent) != monos.end()) { + absent.set(9); + } + const RowBuffer absent_row = row_of(absent); + BOOST_TEST(!store.find(Key{.row = absent_row.view()}).has_value()); + BOOST_TEST(!store.find(Key{.spilled = &absent}).has_value()); +} + +// The codes array's element width follows slots_per_row_, since a codes word only ever sets bits below +// 2 * slots_per_row_. The row array is the operator's largest, and rows are payload -- never a hash +// input, never serialized -- so a narrowing here changes no term and no energy and a baseline diff +// cannot see it. This is the footprint gate. +// memory_bytes() - slack_bytes() is the *used* part of the arrays, which makes the figure exact rather +// than allocator-dependent. +BOOST_AUTO_TEST_CASE(codes_width_follows_the_slot_count) { + constexpr size_t kNumModes = 128; + constexpr size_t kRows = 500; + // 2 bytes of codes at up to 8 slots, 4 up to 16, 8 above -- with the mode lanes constant per slot. + const std::array, 3> kCases{{{8, 2}, {16, 4}, {17, 8}}}; + + for (const auto &[slots, codes_bytes] : kCases) { + SparseRowStore store(slots); + for (size_t i = 0; i < kRows; ++i) { + // One occupied mode per row: any slot count <= slots works, but staying at one keeps every + // row off the overflow side-map, whose bytes are counted separately and would blur this. + Monomial mono; + mono.set(2 * (i % kNumModes)); + store.push_back(mono); + } + const size_t expected = kRows * ((slots * sizeof(RowMode)) + codes_bytes); + BOOST_TEST(store.memory_bytes() - store.slack_bytes() == expected); + } +} + +// The narrowed storage must be invisible above the seam: a store at each codes width has to hold and +// return the same rows, hash them the same way and find them the same way. +BOOST_AUTO_TEST_CASE(a_narrowed_codes_word_reads_back_unchanged) { + constexpr size_t kNumModes = 128; + std::mt19937_64 rng(20260828); + + std::vector> monos; + for (size_t i = 0; i < 200; ++i) { + // At most 8 occupied modes, so every row fits the narrowest store's slots and none spills. + Monomial mono; + for (size_t k = 0; k < 8; ++k) { + mono.set(rng() % (2 * kNumModes)); + } + if (std::find(monos.begin(), monos.end(), mono) == monos.end()) { + monos.push_back(mono); + } + } + + SparseRowStore narrow(8); + SparseRowStore wide(SparseRowStore::kMaxSlots); + for (const auto &mono : monos) { + narrow.push_back(mono); + wide.push_back(mono); + narrow.emplace(mono, narrow.size() - 1); + wide.emplace(mono, wide.size() - 1); + } + BOOST_REQUIRE(narrow.memory_bytes() < wide.memory_bytes()); + + for (size_t i = 0; i < monos.size(); ++i) { + BOOST_REQUIRE(!narrow.spilled(i)); + BOOST_TEST(narrow.codes(i) == wide.codes(i)); + BOOST_TEST(narrow.row(i) == monos[i]); + BOOST_TEST(narrow.popcount(i) == wide.popcount(i)); + BOOST_TEST(narrow.slot_count(i) == wide.slot_count(i)); + BOOST_TEST(sparse_row_hash(narrow.view(i)) == sparse_row_hash(wide.view(i))); + BOOST_TEST(narrow.find(monos[i]).value() == i); + BOOST_TEST(narrow.find(narrow.view(i)).value() == i); + } +} diff --git a/cpp/tests/store_interchange_tests.cpp b/cpp/tests/store_interchange_tests.cpp new file mode 100644 index 00000000..8f51fc9c --- /dev/null +++ b/cpp/tests/store_interchange_tests.cpp @@ -0,0 +1,218 @@ +// 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. + +// Which of the operator store's consumers actually care which store they are given. The answer decides +// how much the backend swap has to touch, so it is asserted rather than reasoned about: +// +// InvertedIndex -- no: it reads rows only through the RowAccess.h accessors, so both stores build the +// same columns, the same parity words and the same tiering. +// MonomialMap -- no: keyed by a dense monomial, which both backends accept as a find() key, so +// init_op_map keeps its key type. +// for_each -- no: SparseRowStore offers OperatorIndex's fn(monomial, row_index) signature. + +#include + +#include +#include +#include +#include + +#include "monoprop/TypeAliases.h" +#include "monoprop/core/Monomial.h" +#include "monoprop/detail/operator/InvertedIndex.h" +#include "monoprop/detail/operator/OperatorIndex.h" +#include "monoprop/detail/operator/SparseRowStore.h" + +#include "RandomMonomial.h" + +using namespace monoprop; +using namespace monoprop::detail; + +namespace { + +// Every observable of a built index, so a difference cannot hide in a field the test forgot. +template +auto columns_agree(const InvertedIndex &a, const InvertedIndex &b) -> bool { + if (a.rows() != b.rows() || a.words() != b.words()) { + return false; + } + for (size_t c = 0; c < InvertedIndex::kNumColumns; ++c) { + if (a.column_is_dense(c) != b.column_is_dense(c)) { + return false; + } + if (a.column_is_dense(c)) { + for (size_t w = 0; w < a.words(); ++w) { + if (a.dense_column_data(c)[w] != b.dense_column_data(c)[w]) { + return false; + } + } + } + else if (a.sparse_column_rows(c) != b.sparse_column_rows(c)) { + return false; + } + } + for (size_t w = 0; w < a.words(); ++w) { + if (a.row_parity_words()[w] != b.row_parity_words()[w]) { + return false; + } + } + return a.tier_memory_bytes() == b.tier_memory_bytes(); +} + +template +auto check_inverted_index_is_store_agnostic(std::mt19937_64 &rng) -> void { + OperatorIndex packed; + // Capacity 4 on purpose, so a good share of rows spill and the index reads them through the + // accessors' overflow path rather than off a codes word. + SparseRowStore sparse(4); + size_t spilled = 0; + for (size_t t = 0; t < 400; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + packed.push_back(mono); + sparse.push_back(mono); + spilled += sparse.spilled(sparse.size() - 1) ? 1 : 0; + } + BOOST_TEST(spilled > 0U); + + InvertedIndex from_packed; + InvertedIndex from_sparse; + from_packed.rebuild(packed); + from_sparse.rebuild(sparse); + BOOST_TEST(columns_agree(from_packed, from_sparse)); + + // append_rows is the incremental path evolution actually takes, so it needs its own comparison -- + // and a like-for-like one. An appended index and a rebuilt index legitimately differ in tiering: + // rebuild() counts every column's postings up front and pre-promotes, while append_rows can only + // promote as rows arrive. That is InvertedIndex's own behaviour and says nothing about the store, + // so what is compared here is append-vs-append. + InvertedIndex appended_sparse; + InvertedIndex appended_packed; + appended_sparse.rebuild(sparse); + appended_packed.rebuild(packed); + for (size_t t = 0; t < 50; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + const size_t base = sparse.size(); + packed.push_back(mono); + sparse.push_back(mono); + appended_sparse.append_rows(sparse, base, 1); + appended_packed.append_rows(packed, base, 1); + } + BOOST_TEST(columns_agree(appended_sparse, appended_packed)); +} + +} // namespace + +BOOST_AUTO_TEST_CASE(store_interchange_inverted_index_is_store_agnostic) { + std::mt19937_64 rng(20260812U); + check_inverted_index_is_store_agnostic<32>(rng); + check_inverted_index_is_store_agnostic<64>(rng); +} + +// MonomialMap is keyed by a dense monomial, and both backends accept one as a find() key, so the map +// needs no change at all when the store is swapped -- which is what is asserted here. +BOOST_AUTO_TEST_CASE(store_interchange_monomial_map_keys_still_resolve) { + std::mt19937_64 rng(1234U); + constexpr size_t kNumModes = 32; + OperatorIndex packed; + SparseRowStore sparse(4); + + MonomialMap pending; + std::vector> stored; + for (size_t t = 0; t < 200; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + // One row per distinct monomial: emplace is insert-or-no-op, so a duplicate would resolve to the + // first row holding it and the coefficient below would name the wrong index. + if (sparse.find(mono).has_value()) { + continue; + } + const size_t i = sparse.size(); + packed.push_back(mono); + sparse.push_back(mono); + packed.emplace(mono, i); + sparse.emplace(mono, i); + stored.push_back(mono); + pending[mono] = static_cast(i); + } + // Terms the store does not hold, which MPOperator::get_operator must leave pending. + size_t absent = 0; + for (size_t t = 0; t < 200; ++t) { + Monomial mono; + for (size_t k = 0; k < 1 + (rng() % 5); ++k) { + mono.set(rng() % (2 * kNumModes)); + } + // emplace, not [], and count only what it inserted: the same absent monomial can be drawn twice, + // and assigning would leave the map smaller than the count. + if (!sparse.find(mono).has_value() && pending.emplace(mono, -1.0).second) { + ++absent; + } + } + BOOST_TEST(absent > 0U); + + // get_operator's loop, run against both stores: a key present in the store drains to its row, one + // absent stays. Both stores must agree on which is which, and on the row. + size_t drained = 0; + size_t retained = 0; + for (const auto &[mono, coeff] : pending) { + const auto in_sparse = sparse.find(mono); + const auto in_packed = packed.find(mono); + BOOST_TEST(in_sparse.has_value() == in_packed.has_value()); + if (in_sparse) { + BOOST_TEST(*in_sparse == *in_packed); + BOOST_TEST(coeff == static_cast(*in_sparse)); + ++drained; + } + else { + ++retained; + } + } + BOOST_TEST(drained == stored.size()); + BOOST_TEST(retained == absent); +} + +// evolved_operator_terms iterates the index, so both stores must offer the same signature and -- since +// the shared RowHashTable fixes slot order for a given insertion sequence -- visit the same row indices +// in the same order. The monomials come out equal even though the two stores hash differently, because +// that order is the *table's*, and both tables saw the same sequence of (index, hash) pairs from their +// own hash. +BOOST_AUTO_TEST_CASE(store_interchange_for_each_visits_every_row) { + std::mt19937_64 rng(4321U); + constexpr size_t kNumModes = 32; + OperatorIndex packed; + SparseRowStore sparse(4); + std::vector> stored; + for (size_t t = 0; t < 200; ++t) { + const auto mono = test_utils::random_monomial(rng, 8); + if (sparse.find(mono).has_value()) { + continue; // one row per distinct monomial, so an index maps to one term + } + const size_t i = sparse.size(); + packed.push_back(mono); + sparse.push_back(mono); + packed.emplace(mono, i); + sparse.emplace(mono, i); + stored.push_back(mono); + } + + std::map> from_packed; + std::map> from_sparse; + packed.for_each([&](const auto &mono, size_t idx) { from_packed.emplace(idx, mono); }); + sparse.for_each([&](const auto &mono, size_t idx) { from_sparse.emplace(idx, mono); }); + + BOOST_REQUIRE(from_packed.size() == stored.size()); + BOOST_REQUIRE(from_sparse.size() == stored.size()); + for (size_t i = 0; i < stored.size(); ++i) { + BOOST_TEST((from_packed.at(i) == stored[i])); + BOOST_TEST((from_sparse.at(i) == stored[i])); + } +} diff --git a/docs/content/docs/building.mdx b/docs/content/docs/building.mdx index 87b80b83..bf4b20e7 100644 --- a/docs/content/docs/building.mdx +++ b/docs/content/docs/building.mdx @@ -19,6 +19,51 @@ mechanism differs by build: The prebuilt wheels published to PyPI (`pip install monoprop`) are also built without MPI, so a from-source build is required for multi-rank runs. +### Other build-time options + +| Option | Default | Effect | +| --- | --- | --- | +| `monoprop_ENABLE_ARCH_FLAGS` | `ON` (`OFF` for the published wheels) | Compile with `-march=native` / `-xHost`. | +| `monoprop_WIDE_TERM_INDEX` | `OFF` | 64-bit term indices, for partitions holding more than ~2^32 terms. | +| `monoprop_ENABLE_MPI` | `OFF` | Multi-rank support, as above. | + +`monoprop_ENABLE_ARCH_FLAGS` also moves the sparse-row crossover, which is derived +from it rather than being an option of its own: `768` modes when the build actually +emits architecture flags, `256` when it does not — which covers the published wheels +and any `Debug` build, where the flags are suppressed. Dense monomials cost one pass +per storage word while sparse rows are flat in the width, so a target without a vector +`popcount` reaches the crossover sooner. The number is build-time because the ISA is, +and there is no crossover you could pin that would be right independently of the flags +you compiled with. + +Expect about one 32-mode block of machine dependence either way — near the crossing the +two backends are within a few percent of each other, so a threshold that is one block +off costs very little. If you are running a width close to it and care, measure your own +workload both ways with `monoprop_ROW_STORE`, which forces a backend per run and needs +no rebuild. + +### Choosing the row backend at run time + +A propagator stores its terms either as dense monomials or as sparse rows (a list of +the modes it occupies plus two bits per mode). The choice is made once, from the mode +count against the crossover above. Set `monoprop_ROW_STORE` to override it for every +propagator in the process: + +| Value | Effect | +| --- | --- | +| `auto` (default, or unset) | Sparse at or above the compiled-in crossover above, dense below. | +| `dense` | Always dense monomials. | +| `sparse` | Always sparse rows. | + +Anything else is an error rather than a silent fallback to `auto`, because the point +of setting it is to know which backend ran. + +Both backends compute the same terms and the same expectation value. They hash rows +differently, so they differ in term *order* — and therefore in floating-point +accumulation order, which is why a cross-backend comparison is a tolerance check and +not a byte diff. Overriding is mainly a testing and benchmarking tool; leave it unset +in production. + ## Prerequisites - a C++23-compliant compiler; on Linux the minimum supported versions are **GCC 14** and **Clang 18** diff --git a/docs/content/docs/testing.mdx b/docs/content/docs/testing.mdx index 04ea8da5..55ea072c 100644 --- a/docs/content/docs/testing.mdx +++ b/docs/content/docs/testing.mdx @@ -87,8 +87,9 @@ just test-wide # rebuilds via uv sync with monoprop_WIDE_TERM_INDEX=ON, then r CTest runs each Boost case as its own process, so an MPI build's `MPI_Init` probes every fabric device per case, whether or not the test sends anything. `monoprop_TEST_EXCLUDE_MPI_FABRIC=ON` -(default) skips that probe for the single-process `serial` variants only; multi-rank variants -keep the full component set, since they exchange real messages. Turn it off by adding +(default) skips that probe for the single-process `serial` variants only — including the +`sparse-rows` half of them, which is a second world-size-1 launch of the same case; multi-rank +variants keep the full component set, since they exchange real messages. Turn it off by adding `-Dmonoprop_TEST_EXCLUDE_MPI_FABRIC=OFF` to the `SKBUILD_CMAKE_ARGS` MPI build above. CTest registers every Boost case individually as a `serial` variant. When the @@ -98,6 +99,31 @@ one entry per rank count rather than per case, because the ranks have to reach the same collectives. Select either group with `ctest --test-dir build/editable/Release -L serial` (or `-L mpi-2`). +### The sparse-rows variants + +Every case is also registered a second time with `monoprop_ROW_STORE=sparse`, +labelled `sparse-rows`, so the whole suite runs against both row backends. This is +not redundancy: every fixture and oracle here is well below the automatic crossover +(`monoprop_SPARSE_ROW_MIN_MODES`), so without forcing it the support-form backend +would be compiled and never executed — and it is the one that ships for wide systems. +`cpp/tests/row_store_selection_tests.cpp` asserts the variable actually reached the +propagator, so the extra pass cannot silently be a second dense run. Run just that +configuration with `-L sparse-rows`, and its Python counterpart with +`just test-sparse-rows`. + +The MPI `sparse-rows` variants use a separate rank list, +`monoprop_MPI_SPARSE_ROWS_TEST_PROCS` (default `2`), so growing +`monoprop_MPI_TEST_PROCS` for dense coverage does not also multiply how many +sparse-row `mpiexec` launches run. + +Expect it to be *slower*: at fixture widths a monomial is one word, which is where +dense rows win. That is the crossover working, not a regression. + +The two backends agree on term sets and expectation values but not on term *order*, +so cross-backend checks use `tools/capture-baseline.py --compare` with a tolerance +(`just diff-baseline-sparse`), never the byte diff that `just diff-baseline` applies +within one backend. + ### Golden baselines A refactor that is meant to change no result is checked against a captured @@ -126,6 +152,12 @@ The capture is rank-local, so under `mpiexec` each rank writes its own `manifest-rank.json`. Output lands in `.baseline-capture/`, which is gitignored. +`just diff-baseline-sparse` is that tolerant comparison applied to the row +backends: it captures with `monoprop_ROW_STORE=sparse` and checks the result +against the golden capture as term sets plus a relative tolerance. The sparse rows +hash differently, so they accumulate in a different order by design — a byte diff +there would fail on every run. + ## Adding tests diff --git a/justfile b/justfile index 737a844b..02a40006 100644 --- a/justfile +++ b/justfile @@ -44,6 +44,20 @@ diff-baseline AGAINST='golden': just capture-baseline candidate diff -rq "{{ baseline_dir }}/{{ AGAINST }}" "{{ baseline_dir }}/candidate" +# Capture with the support-form row backend forced (monoprop_ROW_STORE=sparse) and check it against a +# stored capture as term *sets* plus a relative tolerance. Not a byte diff: the sparse rows hash +# differently, so they accumulate in a different order on purpose. This is how the two backends are +# held equivalent on the fixtures, every one of which is below the automatic crossover. +diff-baseline-sparse AGAINST='golden' TOL='1e-10': + rm -rf "{{ baseline_dir }}/sparse" + monoprop_ROW_STORE=sparse uv run --no-sync python tools/capture-baseline.py --out "{{ baseline_dir }}/sparse" + uv run --no-sync python tools/capture-baseline.py --compare "{{ baseline_dir }}/{{ AGAINST }}" "{{ baseline_dir }}/sparse" --tol "{{ TOL }}" + +# The Python suite with the support-form row backend forced, the counterpart of ctest's +# `-L sparse-rows` variants. +test-sparse-rows: + monoprop_ROW_STORE=sparse uv run --no-sync python -m pytest -m "not mpi" + # MPI is off by default in source builds, so build an MPI-enabled editable install # first, then run the suite under mpiexec with --no-sync (avoids a per-rank resync). # Pass RANKS as either a single integer or a semicolon-separated list (e.g. "1;2;4"). diff --git a/packages/monoprop-bench-tools/src/monoprop_bench_tools/report.py b/packages/monoprop-bench-tools/src/monoprop_bench_tools/report.py index 30410117..75c977bd 100644 --- a/packages/monoprop-bench-tools/src/monoprop_bench_tools/report.py +++ b/packages/monoprop-bench-tools/src/monoprop_bench_tools/report.py @@ -193,6 +193,15 @@ def _fmt_cpus(meta: dict) -> str: return f"{logical}/{physical}" +def _fmt_row_store(meta: dict) -> str: + """Render the row backend as ``asked → ran`` (the two differ whenever the setting is ``auto``).""" + asked = meta.get("monoprop_row_store", "—") + ran = meta.get("row_store_effective") + if ran is None: + return str(asked) + return str(asked) if asked == ran else f"{asked} → {ran}" + + def _config_table(labels: list[str], results: dict[str, dict]) -> list[str]: """Render the run-configuration table (one row per run label).""" metas = {lbl: results.get(lbl, {}).get("meta", {}) for lbl in labels} @@ -208,6 +217,7 @@ def _config_table(labels: list[str], results: dict[str, dict]) -> list[str]: "Ranks/node", "Partitions (requested)", "monoprop threads", + "Row store", "CPUs (logical/physical)", "Host", ] @@ -222,6 +232,7 @@ def _config_table(labels: list[str], results: dict[str, dict]) -> list[str]: str(metas[label].get("ranks_per_node", "—")), str(metas[label].get("partitions_env", "—")), str(metas[label].get("monoprop_threads", "default")), + _fmt_row_store(metas[label]), _fmt_cpus(metas[label]), str(metas[label].get("hostname", "—")), ] diff --git a/src/monoprop/bindings/binder.h b/src/monoprop/bindings/binder.h index 4076c7ff..4b63a5a3 100644 --- a/src/monoprop/bindings/binder.h +++ b/src/monoprop/bindings/binder.h @@ -223,6 +223,10 @@ auto bind_monomial_propagator(nb::module_ &mod) -> void { &MonomialPropagator::logical_num_modes, "Number of modes the operator actually uses"); + cls.def_prop_ro("rows_are_sparse", + &MonomialPropagator::rows_are_sparse, + "Whether this rank's rows live in the support-form backend (see monoprop_ROW_STORE)"); + cls.def_prop_ro_static( "storage_num_modes", [](nb::handle /*unused*/) { return MonomialPropagator::storage_num_modes; },