From 3c226611adfe13b452dc4ed15571249c1ab4d066 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sat, 29 Aug 2026 08:35:21 +0000 Subject: [PATCH 1/2] =?UTF-8?q?feat(cpp):=20=E2=9C=A8=20add=20the=20suppor?= =?UTF-8?q?t-form=20row=20store=20and=20its=20codes=20algebra?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit SparseRowStore is a third row representation: each row is a fixed-width list of the modes it occupies plus one word carrying two bits per occupied mode. CodesAlgebra.h is the structural algebra over that word -- 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. Not yet selected by anything. MPOperator still holds a single OperatorIndex, and neither RowAccess nor the layer-build scan knows this store exists; monoprop_SPARSE_ROW_MIN_MODES is defined and preferred_for_modes() reads it, but no caller asks. Wiring it in is the next change. Two pieces of shared machinery move to RowHashTable.h beside the index both stores key rows through: geometric_row_capacity, the 1.5x growth rule, and spilled_rows_bytes, the side-map accounting. A single number in two stores would otherwise be correctable in one and not the other, which would skew operator_memory_breakdown() for one backend only. splitmix_finalize comes out of SplitmixHash>::mix, whose body never depended on the width: the row hash accumulates over (mode, code) slots and so needs the mix with no Bitset in hand. That value routes MPI ownership, so having one definition rather than two matters more than the tidiness. Verified out of tree against the compile-time-width engine: 3000 random monomials round-trip through row(), popcount(), find(), find_batch(), for_each_position(), clone() and resized(); and over 5000 monomials the codes forms of cutoff_sums, is_paired, pauli_y_count, pauli_anticommutes and pauli_rotation_sign agree with the dense ones exactly. OperatorIndex comes out inert after the two helpers move -- identical iteration order, find results, clone order and memory_bytes. Assisted-by: ClaudeCode:claude-opus-5 --- CMakeLists.txt | 15 + cpp/monoprop/Bitset.h | 20 +- cpp/monoprop/CMakeLists.txt | 2 + cpp/monoprop/algebra/CMakeLists.txt | 1 + cpp/monoprop/algebra/CodesAlgebra.h | 317 ++++++++ cpp/monoprop/detail/operator/CMakeLists.txt | 1 + cpp/monoprop/detail/operator/OperatorIndex.h | 7 +- cpp/monoprop/detail/operator/RowHashTable.h | 18 + cpp/monoprop/detail/operator/SparseRowStore.h | 766 ++++++++++++++++++ 9 files changed, 1134 insertions(+), 13 deletions(-) create mode 100644 cpp/monoprop/algebra/CodesAlgebra.h create mode 100644 cpp/monoprop/detail/operator/SparseRowStore.h diff --git a/CMakeLists.txt b/CMakeLists.txt index c98b529e..6ee35ca3 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,7 @@ 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/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/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/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/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..2d08e32c --- /dev/null +++ b/cpp/monoprop/detail/operator/SparseRowStore.h @@ -0,0 +1,766 @@ +// 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 From 9a6284c25202be4362ba7b743b54cf0f7823e008 Mon Sep 17 00:00:00 2001 From: Ubuntu Date: Sat, 29 Aug 2026 15:36:25 +0000 Subject: [PATCH 2/2] =?UTF-8?q?feat(cpp):=20=E2=9C=A8=20select=20the=20sup?= =?UTF-8?q?port-form=20row=20store=20per=20propagator?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Wires SparseRowStore and CodesAlgebra into the engine: a propagator now picks a row backend once at construction and everything from the scan to the inserts is templated on it. 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, since the scan asks the store for a row per anticommuting term. LayerBuildEngine, fused_find_and_collect and probe_incoming_queries take the store as a template parameter; off that path the forwarding accessors pay one well-predicted branch. There is no accessor handing out a store, because there is no one type to hand out, so MonomialPropagator::indexing() is replaced by for_each_term() and num_local_terms(). The choice comes from the mode count against monoprop_SPARSE_ROW_MIN_MODES; monoprop_ROW_STORE=dense|sparse forces it process-wide, and an unrecognized value throws rather than falling back to auto -- the point of setting it is to know which backend ran. Every C++ case is registered a second time with the sparse backend forced (the sparse-rows ctest label, 541 cases where there were 270): every fixture is below the crossover, so without that the support-form backend would ship untested. row_store_selection_tests.cpp is what fails if the variable stops reaching the propagator. Gates: the dense capture is byte-identical to main's golden baseline, and the sparse capture agrees with it as term sets plus rtol 1e-10 (just diff-baseline-sparse). Assisted-by: ClaudeCode:claude-opus-5 --- AGENTS.md | 42 +- CMakeLists.txt | 5 +- README.md | 4 + benches/conftest.py | 24 + cpp/include/monoprop/MonomialPropagator.h | 41 +- cpp/monoprop/detail/EnvConfig.h | 33 +- .../detail/evolution/layer_build/Engine.h | 161 +++--- .../detail/evolution/layer_build/Resolve.h | 26 +- .../detail/evolution/layer_build/Scan.h | 23 +- .../MonomialPropagator.inl | 52 +- cpp/monoprop/detail/operator/MPOperator.h | 125 +++- cpp/monoprop/detail/operator/RowAccess.h | 47 +- cpp/monoprop/detail/operator/SparseRowStore.h | 3 +- cpp/tests/RandomMonomial.h | 49 ++ cpp/tests/TestOperator.h | 47 ++ cpp/tests/boost-test.cmake | 7 + cpp/tests/boostAddTests.cmake | 87 +++ cpp/tests/codes_algebra_tests.cpp | 288 ++++++++++ cpp/tests/codes_product_tests.cpp | 267 +++++++++ cpp/tests/evolution_detail_tests.cpp | 26 +- cpp/tests/mp_operator_tests.cpp | 36 +- cpp/tests/partition_equivalence_tests.cpp | 2 +- cpp/tests/pauli_build_layer_tests.cpp | 4 +- cpp/tests/row_accessor_tests.cpp | 36 +- cpp/tests/row_store_selection_tests.cpp | 92 +++ cpp/tests/simulator_copy_tests.cpp | 9 +- cpp/tests/sparse_index_tests.cpp | 314 +++++++++++ cpp/tests/sparse_row_store_tests.cpp | 532 ++++++++++++++++++ cpp/tests/store_interchange_tests.cpp | 218 +++++++ docs/content/docs/building.mdx | 45 ++ docs/content/docs/testing.mdx | 36 +- justfile | 14 + .../src/monoprop_bench_tools/report.py | 11 + src/monoprop/bindings/binder.h | 4 + 34 files changed, 2504 insertions(+), 206 deletions(-) create mode 100644 cpp/tests/RandomMonomial.h create mode 100644 cpp/tests/TestOperator.h create mode 100644 cpp/tests/codes_algebra_tests.cpp create mode 100644 cpp/tests/codes_product_tests.cpp create mode 100644 cpp/tests/row_store_selection_tests.cpp create mode 100644 cpp/tests/sparse_index_tests.cpp create mode 100644 cpp/tests/sparse_row_store_tests.cpp create mode 100644 cpp/tests/store_interchange_tests.cpp 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 6ee35ca3..980deaa6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -106,7 +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 + " 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/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/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/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/SparseRowStore.h b/cpp/monoprop/detail/operator/SparseRowStore.h index 2d08e32c..bf4c42e2 100644 --- a/cpp/monoprop/detail/operator/SparseRowStore.h +++ b/cpp/monoprop/detail/operator/SparseRowStore.h @@ -334,8 +334,7 @@ class SparseRowStore { // 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) { + [[gnu::always_inline]] auto with_codes(this Self &&self, F &&f) -> decltype(auto) { switch (self.codes_width_) { case CodesWidth::Narrow: return f(self.codes16_); 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; },