diff --git a/.gitignore b/.gitignore index a728f93..0034639 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,6 @@ academic/**/*.blg academic/**/*.log academic/**/*.out academic/**/*.run.xml + +# AI +.kilo/plans/* diff --git a/CMakeLists.txt b/CMakeLists.txt index 01444ac..e5eeaec 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -34,10 +34,15 @@ endif() # Build options # --------------------------------------------------------------------------- option(PIXIE_TESTS "Build unit tests" ON) -option(PIXIE_BENCHMARKS "Build benchmarks (includes comparison benchmarks against third-party libraries)" OFF) +option(PIXIE_BENCHMARKS "Build benchmarks" OFF) +option(PIXIE_THIRD_PARTY_BACKENDS "Build optional third-party backend integrations" OFF) option(PIXIE_DIAGNOSTICS "Include diagnostic logs" OFF) option(PIXIE_DOCS "Build Doxygen documentation" OFF) +if(PIXIE_THIRD_PARTY_BACKENDS) + add_compile_definitions(SDSL_SUPPORT) +endif() + if(PIXIE_DIAGNOSTICS) add_compile_definitions(PIXIE_DIAGNOSTICS) set(PIXIE_DIAGNOSTICS_LIBS spdlog::spdlog_header_only) @@ -67,10 +72,16 @@ if(PIXIE_BENCHMARKS) GIT_REPOSITORY https://github.com/google/benchmark.git GIT_TAG v1.9.4 ) - set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "Disable Google Benchmark tests") + set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "Disable Google Benchmark tests" FORCE) + set(BENCHMARK_ENABLE_GTEST_TESTS OFF CACHE BOOL "Disable Google Benchmark tests" FORCE) + set(BENCHMARK_ENABLE_INSTALL OFF CACHE BOOL "Disable Google Benchmark install targets" FORCE) + set(BENCHMARK_INSTALL_DOCS OFF CACHE BOOL "Disable Google Benchmark install docs" FORCE) FetchContent_MakeAvailable(googlebenchmark) +endif() - # Third-party libraries for comparison benchmarks +if(PIXIE_THIRD_PARTY_BACKENDS) + set(PASTA_BIT_VECTOR_BUILD_TESTS OFF CACHE BOOL "Disable pasta::bit_vector tests" FORCE) + set(PASTA_BIT_VECTOR_BUILD_BENCHMARKS OFF CACHE BOOL "Disable pasta::bit_vector benchmarks" FORCE) FetchContent_Declare( pasta_bit_vector GIT_REPOSITORY https://github.com/pasta-toolbox/bit_vector.git @@ -96,6 +107,8 @@ endif() if(PIXIE_TESTS) if(NOT TARGET gtest_main) + set(BUILD_GMOCK OFF CACHE BOOL "Disable GoogleMock" FORCE) + set(INSTALL_GTEST OFF CACHE BOOL "Disable GoogleTest install targets" FORCE) FetchContent_Declare( googletest GIT_REPOSITORY https://github.com/google/googletest.git @@ -139,6 +152,10 @@ if(PIXIE_TESTS) gtest gtest_main ${PIXIE_DIAGNOSTICS_LIBS}) + if(PIXIE_THIRD_PARTY_BACKENDS) + target_include_directories(test_rmm + PRIVATE ${sdsl_lite_SOURCE_DIR}/include) + endif() add_executable(louds_tree_tests src/tests/louds_tree_tests.cpp) @@ -170,9 +187,11 @@ if(PIXIE_BENCHMARKS) target_link_libraries(benchmarks benchmark benchmark_main - pasta_bit_vector ${PIXIE_DIAGNOSTICS_LIBS}) - target_compile_definitions(benchmarks PRIVATE PIXIE_THIRD_PARTY_BENCHMARKS) + if(PIXIE_THIRD_PARTY_BACKENDS) + target_link_libraries(benchmarks pasta_bit_vector) + target_compile_definitions(benchmarks PRIVATE PIXIE_THIRD_PARTY_BENCHMARKS) + endif() add_executable(bench_rmm src/benchmarks/bench_rmm.cpp) @@ -182,15 +201,17 @@ if(PIXIE_BENCHMARKS) benchmark ${PIXIE_DIAGNOSTICS_LIBS}) - add_executable(bench_rmm_sdsl - src/benchmarks/bench_rmm_sdsl.cpp) - target_include_directories(bench_rmm_sdsl - PUBLIC include - PRIVATE ${sdsl_lite_SOURCE_DIR}/include) - target_link_libraries(bench_rmm_sdsl - PRIVATE - benchmark - ${PIXIE_DIAGNOSTICS_LIBS}) + if(PIXIE_THIRD_PARTY_BACKENDS) + add_executable(bench_rmm_sdsl + src/benchmarks/bench_rmm_sdsl.cpp) + target_include_directories(bench_rmm_sdsl + PUBLIC include + PRIVATE ${sdsl_lite_SOURCE_DIR}/include) + target_link_libraries(bench_rmm_sdsl + PRIVATE + benchmark + ${PIXIE_DIAGNOSTICS_LIBS}) + endif() add_executable(louds_tree_benchmarks src/benchmarks/louds_tree_benchmarks.cpp) @@ -234,6 +255,16 @@ if(PIXIE_DOCS) FetchContent_GetProperties(doxygen-awesome-css SOURCE_DIR AWESOME_CSS_DIR) + set(PIXIE_DOC_INPUT + "${CMAKE_CURRENT_SOURCE_DIR}/include" + "${CMAKE_CURRENT_SOURCE_DIR}/README.md") + if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/src/docs/benchmark_results.md") + list(APPEND PIXIE_DOC_INPUT + "${CMAKE_CURRENT_SOURCE_DIR}/src/docs/benchmark_results.md") + endif() + string(JOIN " \\\n " PIXIE_DOC_INPUT + ${PIXIE_DOC_INPUT}) + set(DOXYFILE_IN ${CMAKE_CURRENT_SOURCE_DIR}/src/docs/Doxyfile.in) set(DOXYFILE_OUT ${CMAKE_CURRENT_BINARY_DIR}/docs/Doxyfile) configure_file(${DOXYFILE_IN} ${DOXYFILE_OUT} @ONLY) diff --git a/CMakePresets.json b/CMakePresets.json index cae9c15..67438e4 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -19,7 +19,9 @@ "inherits": "base", "binaryDir": "${sourceDir}/build/debug", "cacheVariables": { - "CMAKE_BUILD_TYPE": "Debug" + "CMAKE_BUILD_TYPE": "Debug", + "PIXIE_BENCHMARKS": "OFF", + "PIXIE_TESTS": "ON" } }, { @@ -28,7 +30,9 @@ "inherits": "base", "binaryDir": "${sourceDir}/build/release", "cacheVariables": { - "CMAKE_BUILD_TYPE": "Release" + "CMAKE_BUILD_TYPE": "Release", + "PIXIE_BENCHMARKS": "OFF", + "PIXIE_TESTS": "ON" } }, { @@ -38,9 +42,22 @@ "binaryDir": "${sourceDir}/build/release", "cacheVariables": { "CMAKE_BUILD_TYPE": "Release", + "PIXIE_TESTS": "OFF", "PIXIE_BENCHMARKS": "ON" } }, + { + "name": "benchmarks-third-party", + "displayName": "Benchmarks with third-party backends", + "inherits": "base", + "binaryDir": "${sourceDir}/build/release-third-party", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Release", + "PIXIE_BENCHMARKS": "ON", + "PIXIE_TESTS": "OFF", + "PIXIE_THIRD_PARTY_BACKENDS": "ON" + } + }, { "name": "benchmarks-diagnostic", "displayName": "Benchmarks diagnostic build", @@ -50,6 +67,7 @@ "BENCHMARK_ENABLE_LIBPFM": "ON", "CMAKE_BUILD_TYPE": "RelWithDebInfo", "PIXIE_DIAGNOSTICS": "ON", + "PIXIE_TESTS": "OFF", "PIXIE_BENCHMARKS": "ON" } }, @@ -90,17 +108,63 @@ { "name": "debug", "displayName": "Build Debug", - "configurePreset": "debug" + "configurePreset": "debug", + "targets": [ + "unittests", + "benchmark_tests", + "test_rmm", + "louds_tree_tests", + "excess_positions_tests" + ] }, { "name": "release", "displayName": "Build Release", - "configurePreset": "release" + "configurePreset": "release", + "targets": [ + "unittests", + "benchmark_tests", + "test_rmm", + "louds_tree_tests", + "excess_positions_tests" + ] + }, + { + "name": "benchmarks-all", + "displayName": "Build Benchmarks", + "configurePreset": "benchmarks-all", + "targets": [ + "benchmarks", + "bench_rmm", + "louds_tree_benchmarks", + "alignment_comparison", + "excess_positions_benchmarks" + ] + }, + { + "name": "benchmarks-third-party", + "displayName": "Build Benchmarks with third-party backends", + "configurePreset": "benchmarks-third-party", + "targets": [ + "benchmarks", + "bench_rmm", + "bench_rmm_sdsl", + "louds_tree_benchmarks", + "alignment_comparison", + "excess_positions_benchmarks" + ] }, { "name": "benchmarks-diagnostic", "displayName": "Benchmarks diagnostic", - "configurePreset": "benchmarks-diagnostic" + "configurePreset": "benchmarks-diagnostic", + "targets": [ + "benchmarks", + "bench_rmm", + "louds_tree_benchmarks", + "alignment_comparison", + "excess_positions_benchmarks" + ] }, { "name": "docs", @@ -113,12 +177,26 @@ { "name": "coverage", "displayName": "Build Coverage", - "configurePreset": "coverage" + "configurePreset": "coverage", + "targets": [ + "unittests", + "benchmark_tests", + "test_rmm", + "louds_tree_tests", + "excess_positions_tests" + ] }, { "name": "asan", "displayName": "Build AddressSanitizer", - "configurePreset": "asan" + "configurePreset": "asan", + "targets": [ + "unittests", + "benchmark_tests", + "test_rmm", + "louds_tree_tests", + "excess_positions_tests" + ] } ] } diff --git a/README.md b/README.md index 40c5a63..9537bfe 100644 --- a/README.md +++ b/README.md @@ -51,7 +51,7 @@ cmake -B build/release -DCMAKE_BUILD_TYPE=Release cmake --build build/release -j ``` -Tests are enabled by default (`PIXIE_TESTS=ON`). Benchmarks are opt-in; enable with `-DPIXIE_BENCHMARKS=ON` or configure with the `benchmarks-all` preset, you can use `benchmark-diagnostic` preset for performance diagnostics (Release with debug info + performance counters support). +Tests are enabled by default (`PIXIE_TESTS=ON`). Benchmarks are opt-in; enable with `-DPIXIE_BENCHMARKS=ON` or configure with the `benchmarks-all` preset. Use `benchmarks-third-party` for comparison backends such as sdsl-lite, and `benchmarks-diagnostic` for performance diagnostics (Release with debug info + performance counters support). --- @@ -114,20 +114,57 @@ Benchmarks are random 50/50 0-1 bitvectors up to $2^{34}$ bits. ./build/release/benchmarks ``` +Write JSON and plot size-scaled benchmark curves with a log-scaled x-axis: + +```sh +./build/release/benchmarks --benchmark_out=bitvector_bench.json --benchmark_out_format=json +python3 scripts/plot_size_benchmarks.py bitvector_bench.json -o graphs/bitvector_size.png --size-key n +``` + +### Excess Positions + +```sh +./build/release/excess_positions_benchmarks --benchmark_out=excess_positions.json --benchmark_out_format=json +python3 scripts/excess_benchmark_table.py excess_positions.json -o src/docs/excess_positions_benchmark_results.md +``` + +Generated benchmark documentation can be written to `src/docs/benchmark_results.md`; +the documentation pipeline does not run benchmarks. + ### RmM Tree ```sh ./build/release/bench_rmm ``` -For comparison with range min-max tree implementation from [sdsl-lite](https://github.com/simongog/sdsl-lite) (Release build required; use the release preset or `-DCMAKE_BUILD_TYPE=Release`): +For focused runs, `bench_rmm` accepts `--ops` with a comma-separated operation list. The benchmark harness only builds the query pools needed by the selected operations, so subset runs avoid much of the setup cost: + +```sh +./build/release/bench_rmm --ops=rank1,select1 --benchmark_out=rmm_rank_select.json --benchmark_out_format=json +``` + +By default, RmM benchmarks step through sizes by powers of two. Use `--per_octave=` for finer sampling between adjacent powers of two, or `--explicit_sizes=` for an exact size list. + +Google Benchmark filters are also used to limit RmM setup when `--ops` is not provided: + +```sh +./build/release/bench_rmm --benchmark_filter='^rank1$' --benchmark_out=rmm_rank1.json --benchmark_out_format=json +``` + +For comparison with range min-max tree implementation from [sdsl-lite](https://github.com/simongog/sdsl-lite), use the third-party benchmark preset. This defines `SDSL_SUPPORT` and builds `bench_rmm_sdsl`: ```bash +cmake --preset benchmarks-third-party +cmake --build --preset benchmarks-third-party sudo cpupower frequency-set --governor performance -./build/release/bench_rmm_sdsl --benchmark_out=rmm_bench_sdsl.json +./build/release-third-party/bench_rmm_sdsl --benchmark_out=rmm_bench_sdsl.json ``` -For visualization, write the JSON output to a file using `--benchmark_out=` (e.g. `./build/release/bench_rmm --benchmark_out=rmm_bench.json`) and plot it with `scripts/plot_rmm.py` (add `--sdsl-json rmm_bench_sdsl.json` for comparison). +For visualization, write the JSON output to a file using `--benchmark_out=` (e.g. `./build/release/bench_rmm --benchmark_out=rmm_bench.json`) and plot it with `scripts/plot_rmm.py` (add `--sdsl-json rmm_bench_sdsl.json` for per-operation sdsl-lite comparison plots). For size-scaled tree plots, use: + +```sh +python3 scripts/plot_size_benchmarks.py rmm_bench.json -o graphs/rmm_size.png --size-key N +``` --- @@ -152,8 +189,10 @@ int main() { ```cpp #include -#include +#include #include +#include +#include using namespace pixie; @@ -166,7 +205,15 @@ int main() { // └─ C // └─ c1 std::string bits = "11101001011000"; - RmMTree t(bits); + std::vector words((bits.size() + 63) / 64); + for (std::size_t i = 0; i < bits.size(); ++i) { + if (bits[i] == '1') { + words[i / 64] |= std::uint64_t{1} << (i % 64); + } + } + + // RmMTree is non-owning: keep words alive and immutable while using t. + RmMTree t(words, bits.size()); std::cout << "close(1): " << t.close(1) << "\n"; // expected 6 (A) std::cout << "open(3): " << t.open(3) << "\n"; // expected 2 (a1) diff --git a/include/pixie/rmm_base.h b/include/pixie/rmm_base.h new file mode 100644 index 0000000..3899a50 --- /dev/null +++ b/include/pixie/rmm_base.h @@ -0,0 +1,200 @@ +#pragma once + +#include +#include + +namespace pixie { + +/** + * @brief CRTP facade for rank/select and range min-max tree operations. + * + * Positions are zero-based unless explicitly documented otherwise. Operations + * returning positions use npos when the requested value does not exist or the + * query is outside the valid domain. Balanced-parentheses navigation follows + * SDSL-style zero-based parenthesis indexing: close/open/enclose accept and + * return bit positions, not prefix-boundary positions. + */ +template +class RmMBase { + public: + /** + * @brief Sentinel returned by position queries when no valid answer exists. + */ + static constexpr std::size_t npos = std::numeric_limits::max(); + + /** + * @brief Number of bits in the represented sequence. + */ + std::size_t size() const { return impl().size_impl(); } + + /** + * @brief Count 1 bits in the prefix [0, end_position). + */ + std::size_t rank1(std::size_t end_position) const { + return impl().rank1_impl(end_position); + } + + /** + * @brief Count 0 bits in the prefix [0, end_position). + */ + std::size_t rank0(std::size_t end_position) const { + return impl().rank0_impl(end_position); + } + + /** + * @brief Return the zero-based position of the rank-th 1 bit. + * @param rank One-based rank of the requested 1 bit. + */ + std::size_t select1(std::size_t rank) const { + return impl().select1_impl(rank); + } + + /** + * @brief Return the zero-based position of the rank-th 0 bit. + * @param rank One-based rank of the requested 0 bit. + */ + std::size_t select0(std::size_t rank) const { + return impl().select0_impl(rank); + } + + /** + * @brief Count "10" bit patterns fully contained in [0, end_position). + * + * The returned count is the number of start positions p such that + * p + 1 < end_position, bit[p] == 1, and bit[p + 1] == 0. + */ + std::size_t rank10(std::size_t end_position) const { + return impl().rank10_impl(end_position); + } + + /** + * @brief Return the zero-based start position of the rank-th "10" pattern. + * @param rank One-based rank of the requested pattern. + */ + std::size_t select10(std::size_t rank) const { + return impl().select10_impl(rank); + } + + /** + * @brief Prefix excess on [0, end_position): +1 for 1 bits, -1 for 0 bits. + */ + int excess(std::size_t end_position) const { + return impl().excess_impl(end_position); + } + + /** + * @brief Forward excess search from a prefix boundary. + * + * Finds the first bit position p >= start_position such that + * excess(p + 1) == excess(start_position) + delta. + */ + std::size_t fwdsearch(std::size_t start_position, int delta) const { + return impl().fwdsearch_impl(start_position, delta); + } + + /** + * @brief Backward excess search from a prefix boundary. + * + * Finds the greatest prefix boundary p < start_position such that + * excess(p) == excess(start_position) + delta. + */ + std::size_t bwdsearch(std::size_t start_position, int delta) const { + return impl().bwdsearch_impl(start_position, delta); + } + + /** + * @brief Position of the first minimum relative excess in [range_begin, + * range_end]. + * + * Relative excess is accumulated from zero over the inclusive bit range. + */ + std::size_t range_min_query_pos(std::size_t range_begin, + std::size_t range_end) const { + return impl().range_min_query_pos_impl(range_begin, range_end); + } + + /** + * @brief Minimum relative excess value over [range_begin, range_end]. + * + * Relative excess is accumulated from zero over the inclusive bit range. + */ + int range_min_query_val(std::size_t range_begin, + std::size_t range_end) const { + return impl().range_min_query_val_impl(range_begin, range_end); + } + + /** + * @brief Position of the first maximum relative excess in [range_begin, + * range_end]. + * + * Relative excess is accumulated from zero over the inclusive bit range. + */ + std::size_t range_max_query_pos(std::size_t range_begin, + std::size_t range_end) const { + return impl().range_max_query_pos_impl(range_begin, range_end); + } + + /** + * @brief Maximum relative excess value over [range_begin, range_end]. + * + * Relative excess is accumulated from zero over the inclusive bit range. + */ + int range_max_query_val(std::size_t range_begin, + std::size_t range_end) const { + return impl().range_max_query_val_impl(range_begin, range_end); + } + + /** + * @brief Count positions attaining the minimum relative excess in + * [range_begin, range_end]. + */ + std::size_t mincount(std::size_t range_begin, std::size_t range_end) const { + return impl().mincount_impl(range_begin, range_end); + } + + /** + * @brief Select a position attaining the minimum relative excess in + * [range_begin, range_end]. + * @param rank One-based rank among positions attaining the range minimum. + */ + std::size_t minselect(std::size_t range_begin, + std::size_t range_end, + std::size_t rank) const { + return impl().minselect_impl(range_begin, range_end, rank); + } + + /** + * @brief Matching closing parenthesis for the parenthesis at open_position. + * + * Uses SDSL-style zero-based indexing. If open_position refers to a closing + * parenthesis, the expected result is open_position. + */ + std::size_t close(std::size_t open_position) const { + return impl().close_impl(open_position); + } + + /** + * @brief Matching opening parenthesis for the parenthesis at close_position. + * + * Uses SDSL-style zero-based indexing. If close_position refers to an opening + * parenthesis, the expected result is close_position. + */ + std::size_t open(std::size_t close_position) const { + return impl().open_impl(close_position); + } + + /** + * @brief Closest opening parenthesis strictly enclosing position. + * + * Uses SDSL-style zero-based indexing. If position refers to a closing + * parenthesis, this is equivalent to open(position). + */ + std::size_t enclose(std::size_t open_position) const { + return impl().enclose_impl(open_position); + } + + private: + const Impl& impl() const { return static_cast(*this); } +}; + +} // namespace pixie diff --git a/include/pixie/rmm_tree.h b/include/pixie/rmm_tree.h index fe8bbe7..e26a5b1 100644 --- a/include/pixie/rmm_tree.h +++ b/include/pixie/rmm_tree.h @@ -1,6 +1,7 @@ #pragma once #include #include +#include #include #include @@ -9,7 +10,8 @@ #include #include #include -#include +#include +#include #include namespace pixie { @@ -29,12 +31,14 @@ namespace pixie { * - mincount / minselect (count/selection of minima) * - close / open / enclose (BP navigation) * - * The bitvector is LSB-first inside each 64-bit word. + * The bitvector is LSB-first inside each 64-bit word. RmMTree stores a + * non-owning view of those words; callers must keep the backing storage alive + * and immutable for the lifetime of the tree. */ -class RmMTree { +class RmMTree : public RmMBase { // ------------ bitvector ------------ - std::vector bits; // LSB-first - size_t num_bits = 0; // number of bits + std::span bits; // LSB-first, externally owned + size_t num_bits = 0; // number of bits // ------------ blocking ------------ size_t block_bits = 64; // block size (bits), leaf covers <= block_bits bits @@ -93,44 +97,34 @@ class RmMTree { RmMTree() = default; /** - * @brief Build from a '0'/'1' string. - * @param bit_string Bitstring (characters '0' and '1'). - * @param leaf_block_bits Desired leaf size (power of two, 0 = auto). - * @param max_overhead Max allowed overhead fraction (<0 to disable - * constraint). - * @details Block size priority: (1) respect @p max_overhead, (2) explicit @p - * leaf_block_bits, (3) set to ceil_pow2(log2(num_bits)). - */ - explicit RmMTree(const std::string& bit_string, - const size_t& leaf_block_bits /*0=auto*/ = 0, - const float& max_overhead /*<0=off*/ = -1.0) { - build_from_string(bit_string, leaf_block_bits, max_overhead); - } - - /** - * @brief Build from 64-bit words (LSB-first). - * @param words Array of words holding bits LSB-first. + * @brief Build from a non-owning view of 64-bit words (LSB-first). + * @param words Non-owning view of words holding bits LSB-first. * @param bit_count Number of valid bits. * @param leaf_block_bits Desired leaf size (power of two, 0 = auto). * @param max_overhead Max allowed overhead fraction (<0 to disable * constraint). * @details Block size priority: (1) respect @p max_overhead, (2) explicit @p * leaf_block_bits, (3) set to ceil_pow2(log2(num_bits)). + * + * The caller must keep @p words alive and immutable for the lifetime of this + * tree. */ - explicit RmMTree(const std::vector& words, + explicit RmMTree(std::span words, size_t bit_count, const size_t& leaf_block_bits /*0=auto*/ = 0, const float& max_overhead /*<0=off*/ = -1.0) { build_from_words(words, bit_count, leaf_block_bits, max_overhead); } + size_t size_impl() const { return num_bits; } + // --------- queries: rank/select/excess ---------- /** * @brief Number of ones in prefix [0, @p end_position). * @details Returns 0 for @p end_position == 0. */ - size_t rank1(const size_t& end_position) const { + size_t rank1_impl(const size_t& end_position) const { if (end_position == 0) { return 0; } @@ -153,17 +147,17 @@ class RmMTree { /** * @brief Number of zeros in prefix [0, @p end_position). - * @details Computed as @p end_position - rank1(@p end_position). + * @details Computed as @p end_position - rank1_impl(@p end_position). */ - size_t rank0(const size_t& end_position) const { - return end_position - rank1(end_position); + size_t rank0_impl(const size_t& end_position) const { + return end_position - rank1_impl(end_position); } /** * @brief 1-based select of the @p target_one_rank-th one. * @return Position of @p target_one_rank-th '1' or npos if not found. */ - size_t select1(size_t target_one_rank) const { + size_t select1_impl(size_t target_one_rank) const { if (target_one_rank == 0 || num_bits == 0) { return npos; } @@ -194,7 +188,7 @@ class RmMTree { * @brief 1-based select of the @p target_zero_rank-th zero. * @return Position of @p target_zero_rank-th '0' or npos if not found. */ - size_t select0(size_t target_zero_rank) const { + size_t select0_impl(size_t target_zero_rank) const { if (target_zero_rank == 0 || num_bits == 0) { return npos; } @@ -229,7 +223,7 @@ class RmMTree { * @details Counts p where bit[p]==1 and bit[p+1]==0 with p+1<@p * end_position. */ - size_t rank10(const size_t& end_position) const { + size_t rank10_impl(const size_t& end_position) const { if (end_position <= 1) { return 0; } @@ -262,7 +256,7 @@ class RmMTree { * @brief 1-based select of the @p target_pattern_rank-th "10" start. * @return Position p such that bits[p..p+1]=="10", or npos if not found. */ - size_t select10(size_t target_pattern_rank) const { + size_t select10_impl(size_t target_pattern_rank) const { if (target_pattern_rank == 0 || num_bits == 0) { return npos; } @@ -317,17 +311,17 @@ class RmMTree { /** * @brief Prefix excess on [0, @p end_position): +1 for '1', −1 for '0'. */ - inline int excess(const size_t& end_position) const { - return int64_t(rank1(end_position)) * 2 - int64_t(end_position); + inline int excess_impl(const size_t& end_position) const { + return int64_t(rank1_impl(end_position)) * 2 - int64_t(end_position); } /** * @brief Forward search: first position p ≥ @p start_position where - * excess(p) = excess(@p start_position) + @p delta. + * excess_impl(p) = excess_impl(@p start_position) + @p delta. * @details Scans remainder of current leaf, then descends using precomputed * bounds. Returns npos if no such position exists. */ - size_t fwdsearch(const size_t& start_position, const int& delta) const { + size_t fwdsearch_impl(const size_t& start_position, const int& delta) const { if (start_position >= num_bits) { return npos; } @@ -387,11 +381,11 @@ class RmMTree { /** * @brief Backward search: last position p ≤ @p start_position where - * excess(p) = excess(@p start_position) + @p delta. + * excess_impl(p) = excess_impl(@p start_position) + @p delta. * @details Scans inside the leaf to the left, then climbs to examine left * siblings. Returns npos if no such position exists. */ - size_t bwdsearch(const size_t& start_position, const int& delta) const { + size_t bwdsearch_impl(const size_t& start_position, const int& delta) const { if (start_position > num_bits || start_position == 0) { return npos; } @@ -399,15 +393,16 @@ class RmMTree { // 1) scan inside the block const size_t leaf_block_index = block_of(start_position - 1); const size_t block_begin = leaf_block_index * block_bits; - int leaf_delta = 0; // excess(start_position) - excess(block_begin) + int leaf_delta = + 0; // excess_impl(start_position) - excess_impl(block_begin) const size_t leaf_result = leaf_bwd_bp_simd( leaf_block_index, block_begin, start_position, delta, leaf_delta); if (leaf_result != npos) { return leaf_result; } - // need = target - excess(block_begin) = excess(start_position) + delta - - // excess(block_begin) = leaf_delta + delta + // need = target - excess_impl(block_begin) = excess_impl(start_position) + + // delta - excess_impl(block_begin) = leaf_delta + delta int remaining_delta = leaf_delta + delta; size_t node_index = leaf_index_of(block_begin); size_t segment_base = block_begin; @@ -455,8 +450,8 @@ class RmMTree { * (inclusive). * @return Position of first occurrence of minimum, or npos on invalid range. */ - size_t range_min_query_pos(const size_t& range_begin, - const size_t& range_end) const { + size_t range_min_query_pos_impl(const size_t& range_begin, + const size_t& range_end) const { if (range_begin > range_end || range_end >= num_bits) { return npos; } @@ -555,18 +550,18 @@ class RmMTree { * @brief Value of the minimum prefix excess on [@p range_begin, @p range_end] * relative to @p range_begin. * @details Equivalent to min_{t in [@p range_begin..@p range_end]} - * (excess(t+1) - excess(@p range_begin)). + * (excess_impl(t+1) - excess_impl(@p range_begin)). */ - int range_min_query_val(const size_t& range_begin, - const size_t& range_end) const { + int range_min_query_val_impl(const size_t& range_begin, + const size_t& range_end) const { if (range_begin > range_end || range_end >= num_bits) { return 0; } - size_t min_position = range_min_query_pos(range_begin, range_end); + size_t min_position = range_min_query_pos_impl(range_begin, range_end); if (min_position == npos) { return 0; } - return excess(min_position + 1) - excess(range_begin); + return excess_impl(min_position + 1) - excess_impl(range_begin); } /** @@ -575,8 +570,8 @@ class RmMTree { * (inclusive). * @return Position of first occurrence of maximum, or npos on invalid range. */ - size_t range_max_query_pos(const size_t& range_begin, - const size_t& range_end) const { + size_t range_max_query_pos_impl(const size_t& range_begin, + const size_t& range_end) const { if (range_begin > range_end || range_end >= num_bits) { return npos; } @@ -675,23 +670,24 @@ class RmMTree { * @brief Value of the maximum prefix excess on [@p range_begin, @p * range_end] relative to @p range_begin. */ - int range_max_query_val(const size_t& range_begin, - const size_t& range_end) const { + int range_max_query_val_impl(const size_t& range_begin, + const size_t& range_end) const { if (range_begin > range_end || range_end >= num_bits) { return 0; } - size_t max_position = range_max_query_pos(range_begin, range_end); + size_t max_position = range_max_query_pos_impl(range_begin, range_end); if (max_position == npos) { return 0; } - return excess(max_position + 1) - excess(range_begin); + return excess_impl(max_position + 1) - excess_impl(range_begin); } /** * @brief How many times the minimum prefix excess occurs on [@p range_begin, * @p range_end]. */ - size_t mincount(const size_t& range_begin, const size_t& range_end) const { + size_t mincount_impl(const size_t& range_begin, + const size_t& range_end) const { if (range_begin > range_end || range_end >= num_bits) { return 0; } @@ -774,9 +770,9 @@ class RmMTree { * @return Position or npos if @p target_min_rank exceeds the number of * minima. */ - size_t minselect(const size_t& range_begin, - const size_t& range_end, - size_t target_min_rank) const { + size_t minselect_impl(const size_t& range_begin, + const size_t& range_end, + size_t target_min_rank) const { if (range_begin > range_end || range_end >= num_bits || target_min_rank == 0) { return npos; @@ -930,39 +926,47 @@ class RmMTree { // ----- parentheses navigation (BP) ----- /** - * @brief close(@p open_position): matching ')' for '(' at @p open_position. + * @brief close_impl(@p open_position): matching ')' for '(' at @p + * open_position. + * @todo This method still uses the older boundary-adjusted search convention; + * align it with SDSL-style zero-based parenthesis indexing. * @return Position of matching ')', or npos. */ - inline size_t close(const size_t& open_position) const { + inline size_t close_impl(const size_t& open_position) const { if (open_position >= num_bits) { return npos; } - return fwdsearch(open_position, -1); + return fwdsearch_impl(open_position, -1); } /** - * @brief open(@p close_position): matching '(' for ')' at @p close_position. + * @brief open_impl(@p close_position): matching '(' for ')' at @p + * close_position. + * @todo This method still uses the older boundary-adjusted search convention; + * align it with SDSL-style zero-based parenthesis indexing. * @return Position of matching '(', or npos. */ - inline size_t open(const size_t& close_position) const { + inline size_t open_impl(const size_t& close_position) const { // bwdsearch allows i in [1..num_bits] if (close_position == 0 || close_position > num_bits) { return npos; } - const size_t result = bwdsearch(close_position, 0); + const size_t result = bwdsearch_impl(close_position, 0); return (result == npos ? npos : result + 1); } /** - * @brief enclose(@p position): opening '(' that strictly encloses @p + * @brief enclose_impl(@p position): opening '(' that strictly encloses @p * position. + * @todo This method still uses the older boundary-adjusted search convention; + * align it with SDSL-style zero-based parenthesis indexing. * @return Position of enclosing '(', or npos. */ - inline size_t enclose(const size_t& position) const { + inline size_t enclose_impl(const size_t& position) const { if (position == 0 || position > num_bits) { return npos; } - const size_t result = bwdsearch(position, -2); + const size_t result = bwdsearch_impl(position, -2); return (result == npos ? npos : result + 1); } @@ -2178,24 +2182,6 @@ class RmMTree { return candidate_block_bits; } - /** - * @brief Build internal structures from a 0/1 string. - * @param leaf_block_bits Desired leaf size (0 = auto). - * @param max_overhead Overhead cap (<0 = no cap). - */ - void build_from_string(const std::string& bit_string_input, - const size_t& leaf_block_bits = 0, - const float& max_overhead = -1.0) { - num_bits = bit_string_input.size(); - bits.assign((num_bits + 63) / 64, 0); - for (size_t i = 0; i < num_bits; ++i) { - if (bit_string_input[i] == '1') { - set1(i); - } - } - build(leaf_block_bits, max_overhead); - } - /** * @brief Build internal structures from 64-bit words. * @param words Words with LSB-first bits. @@ -2203,14 +2189,15 @@ class RmMTree { * @param leaf_block_bits Desired leaf size (0 = auto). * @param max_overhead Overhead cap (<0 = no cap). */ - void build_from_words(const std::vector& words, + void build_from_words(std::span words, const size_t& bit_count, const size_t& leaf_block_bits = 0, const float& max_overhead = -1.0) { bits = words; num_bits = bit_count; if (bits.size() * 64 < num_bits) { - bits.resize((num_bits + 63) / 64); + throw std::invalid_argument( + "RmMTree bit_count exceeds the provided word span"); } build(leaf_block_bits, max_overhead); } @@ -2222,13 +2209,6 @@ class RmMTree { return (bits[position >> 6] >> (position & 63)) & 1u; } - /** - * @brief Set bit at position @p position to 1. - */ - inline void set1(const size_t& position) noexcept { - bits[position >> 6] |= (std::uint64_t(1) << (position & 63)); - } - /** * @brief Number of ones in node @p node_index computed from size and total * excess. @@ -2494,8 +2474,8 @@ class RmMTree { return npos; } - // leaf_delta = excess(start_position) - excess(leaf_block_begin) - // = 2*rank1([leaf_begin, start)) - (start - leaf_begin) + // leaf_delta = excess_impl(start_position) - excess_impl(leaf_block_begin) + // = 2*rank1_impl([leaf_begin, start)) - (start - leaf_begin) const int len = int(start_position - leaf_block_begin); const int ones = int(rank1_in_block(leaf_block_begin, start_position)); leaf_delta = ones * 2 - len; diff --git a/include/pixie/rmm_tree_sdsl.h b/include/pixie/rmm_tree_sdsl.h new file mode 100644 index 0000000..3245857 --- /dev/null +++ b/include/pixie/rmm_tree_sdsl.h @@ -0,0 +1,204 @@ +#pragma once + +#ifndef SDSL_SUPPORT +#error "pixie/rmm_tree_sdsl.h requires SDSL_SUPPORT" +#endif + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace pixie { + +class SdslRmMTree : public RmMBase { + public: + using BpSupport = sdsl::bp_support_sada<>; + static constexpr std::size_t npos = RmMBase::npos; + + SdslRmMTree() = default; + + SdslRmMTree(std::span words, + std::size_t bit_count, + std::size_t) { + size_ = bit_count; + const std::size_t valid_words = (bit_count + 63) / 64; + for (std::size_t i = 0; i < valid_words; ++i) { + std::uint64_t word = words[i]; + if (i + 1 == valid_words && (bit_count & 63) != 0) { + word &= (std::uint64_t{1} << (bit_count & 63)) - 1; + } + ones_ += std::popcount(word); + } + zeros_ = size_ - ones_; + + bits_ = sdsl::bit_vector(size_); + for (std::size_t i = 0; i < size_; ++i) { + bits_[i] = (words[i >> 6] >> (i & 63)) & 1u; + } + tree_ = BpSupport(&bits_); + } + + SdslRmMTree(const SdslRmMTree& other) + : size_(other.size_), + ones_(other.ones_), + zeros_(other.zeros_), + bits_(other.bits_) { + reset_support(); + } + + SdslRmMTree& operator=(const SdslRmMTree& other) { + if (this == &other) { + return *this; + } + size_ = other.size_; + ones_ = other.ones_; + zeros_ = other.zeros_; + bits_ = other.bits_; + reset_support(); + return *this; + } + + SdslRmMTree(SdslRmMTree&& other) noexcept + : size_(other.size_), + ones_(other.ones_), + zeros_(other.zeros_), + bits_(std::move(other.bits_)) { + reset_support(); + } + + SdslRmMTree& operator=(SdslRmMTree&& other) noexcept { + if (this == &other) { + return *this; + } + size_ = other.size_; + ones_ = other.ones_; + zeros_ = other.zeros_; + bits_ = std::move(other.bits_); + reset_support(); + return *this; + } + + std::size_t size_impl() const { return size_; } + + std::size_t rank1_impl(std::size_t end_position) const { + if (size_ == 0 || end_position == 0) { + return 0; + } + return tree_.rank(std::min(end_position, size_) - 1); + } + + std::size_t rank0_impl(std::size_t end_position) const { + if (size_ == 0 || end_position == 0) { + return 0; + } + if (end_position >= size_) { + return zeros_; + } + return tree_.preceding_closing_parentheses(end_position); + } + + std::size_t select1_impl(std::size_t rank) const { + if (rank == 0 || rank > ones_) { + return npos; + } + return tree_.select(rank); + } + + std::size_t select0_impl(std::size_t) const { return npos; } + std::size_t rank10_impl(std::size_t) const { return 0; } + std::size_t select10_impl(std::size_t) const { return npos; } + + int excess_impl(std::size_t end_position) const { + if (size_ == 0 || end_position == 0) { + return 0; + } + return tree_.excess(std::min(end_position, size_) - 1); + } + + std::size_t fwdsearch_impl(std::size_t start_position, int delta) const { + if (start_position >= size_) { + return npos; + } + std::int64_t current = excess_impl(start_position); + const std::int64_t target = current + delta; + for (std::size_t position = start_position; position < size_; ++position) { + current += bits_[position] ? 1 : -1; + if (current == target) { + return position; + } + } + return npos; + } + + std::size_t bwdsearch_impl(std::size_t, int) const { return npos; } + + std::size_t range_min_query_pos_impl(std::size_t range_begin, + std::size_t range_end) const { + if (size_ == 0) { + return 0; + } + return tree_.rmq(range_begin, range_end); + } + + int range_min_query_val_impl(std::size_t range_begin, + std::size_t range_end) const { + if (size_ == 0) { + return 0; + } + const auto min_position = tree_.rmq(range_begin, range_end); + if (min_position >= size_) { + return 0; + } + const auto base_excess = + (range_begin == 0 ? 0 : tree_.excess(range_begin - 1)); + return tree_.excess(min_position) - base_excess; + } + + std::size_t range_max_query_pos_impl(std::size_t, std::size_t) const { + return npos; + } + int range_max_query_val_impl(std::size_t, std::size_t) const { return 0; } + std::size_t mincount_impl(std::size_t, std::size_t) const { return 0; } + std::size_t minselect_impl(std::size_t, std::size_t, std::size_t) const { + return npos; + } + + std::size_t close_impl(std::size_t open_position) const { + if (size_ == 0) { + return 0; + } + return tree_.find_close(open_position); + } + + std::size_t open_impl(std::size_t close_position) const { + if (size_ == 0) { + return 0; + } + return tree_.find_open(close_position); + } + + std::size_t enclose_impl(std::size_t open_position) const { + if (size_ == 0) { + return 0; + } + return tree_.enclose(open_position); + } + + private: + void reset_support() { tree_ = BpSupport(&bits_); } + + std::size_t size_{}; + std::size_t ones_{}; + std::size_t zeros_{}; + sdsl::bit_vector bits_; + BpSupport tree_; +}; + +} // namespace pixie diff --git a/scripts/excess_benchmark_table.py b/scripts/excess_benchmark_table.py new file mode 100644 index 0000000..0915341 --- /dev/null +++ b/scripts/excess_benchmark_table.py @@ -0,0 +1,109 @@ +#!/usr/bin/env python3 +"""Render excess-position benchmark JSON as a compact Markdown table.""" + +from __future__ import annotations + +import argparse +import json +from pathlib import Path +from typing import Dict, Iterable, Tuple + + +def _target_from_benchmark(bench: dict) -> str: + if "params" in bench and "X" in bench["params"]: + return str(bench["params"]["X"]) + name = bench.get("name", "") + for suffix in ("_mean", "_median", "_stddev", "_cv"): + name = name.removesuffix(suffix) + marker = "/X:" + if marker in name: + return name.split(marker, 1)[1].split("/", 1)[0] + args = bench.get("args", []) + if args: + return str(args[0]) + return "" + + +def _method_from_benchmark(bench: dict) -> str: + name = bench.get("name", "unknown") + name = name.split("/", 1)[0] + return name.removeprefix("BM_ExcessPositions512").lstrip("_") or "Current" + + +def _is_data_row(bench: dict) -> bool: + name = bench.get("name", "") + if not name.startswith("BM_ExcessPositions512"): + return False + if any(name.endswith(suffix) for suffix in ("_median", "_stddev", "_cv")): + return False + run_type = bench.get("run_type") + return run_type in (None, "iteration", "aggregate") + + +def _collect(benchmarks: Iterable[dict]) -> Dict[Tuple[str, str], Tuple[float, float]]: + rows: Dict[Tuple[str, str], Tuple[float, float]] = {} + for bench in benchmarks: + if not _is_data_row(bench): + continue + name = bench.get("name", "") + if bench.get("run_type") == "aggregate" and not name.endswith("_mean"): + continue + method = _method_from_benchmark(bench).removesuffix("_mean") + target = _target_from_benchmark(bench) + if not target: + continue + real_time = float(bench.get("real_time", 0.0)) + cpu_time = float(bench.get("cpu_time", 0.0)) + rows[(method, target)] = (real_time, cpu_time) + return rows + + +def _sort_target(value: str) -> int: + try: + return int(value) + except ValueError: + return 0 + + +def render_markdown(rows: Dict[Tuple[str, str], Tuple[float, float]]) -> str: + methods = sorted({method for method, _ in rows}) + targets = sorted({target for _, target in rows}, key=_sort_target) + if not methods or not targets: + raise SystemExit("No excess benchmark rows found.") + + lines = ["| Method | " + " | ".join(f"X={x}" for x in targets) + " |"] + lines.append("|---|" + "|".join("---:" for _ in targets) + "|") + for method in methods: + cells = [] + for target in targets: + value = rows.get((method, target)) + cells.append("" if value is None else f"{value[1]:.2f} ns") + lines.append(f"| {method} | " + " | ".join(cells) + " |") + return "\n".join(lines) + + +def main() -> int: + parser = argparse.ArgumentParser( + description="Create a Markdown CPU-time table from excess benchmark JSON." + ) + parser.add_argument("json", type=Path, help="Google Benchmark JSON output") + parser.add_argument( + "-o", + "--output", + type=Path, + help="Optional output Markdown file. Defaults to stdout.", + ) + args = parser.parse_args() + + data = json.loads(args.json.read_text(encoding="utf-8")) + table = render_markdown(_collect(data.get("benchmarks", []))) + if args.output: + args.output.parent.mkdir(parents=True, exist_ok=True) + args.output.write_text(table + "\n", encoding="utf-8") + else: + print(table) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/plot_rmm.py b/scripts/plot_rmm.py index ce88c99..b543d4b 100644 --- a/scripts/plot_rmm.py +++ b/scripts/plot_rmm.py @@ -17,7 +17,17 @@ import argparse import json import os -import pandas as pd +import tempfile +from collections import defaultdict +from statistics import median + +os.environ.setdefault( + "MPLCONFIGDIR", os.path.join(tempfile.gettempdir(), "pixie-matplotlib") +) + +import matplotlib + +matplotlib.use("Agg") import matplotlib.pyplot as plt OPS_ORDER = [ @@ -42,20 +52,59 @@ ] -def load_bench_json(path: str) -> pd.DataFrame: +def clean_name(name: str) -> str: + for suffix in ("_mean", "_median", "_stddev", "_cv"): + name = name.removesuffix(suffix) + return name + + +def load_bench_json(path: str) -> list[dict]: with open(path, "r", encoding="utf-8") as f: data = json.load(f) if isinstance(data, dict) and "benchmarks" in data: data = data["benchmarks"] - df = pd.DataFrame(data) - required = {"name", "cpu_time", "N"} - if not required.issubset(df.columns): - missing = ", ".join(sorted(required - set(df.columns))) - raise ValueError( - f"{path!r}: missing required columns in benchmark JSON: {missing}" - ) - df = df.dropna(subset=["cpu_time", "N"]) - return df + rows = [] + for bench in data: + if not {"name", "cpu_time", "N"}.issubset(bench): + continue + if bench.get("run_type") == "aggregate" and bench.get("aggregate_name") != "mean": + continue + try: + n = int(float(bench["N"])) + cpu_time = float(bench["cpu_time"]) + except (TypeError, ValueError): + continue + if n <= 0 or cpu_time <= 0: + continue + rows.append({"name": clean_name(str(bench["name"])), "N": n, "cpu_time": cpu_time}) + if not rows: + raise ValueError(f"{path!r}: no usable benchmark rows found") + return rows + + +def ops_in(rows: list[dict]) -> set[str]: + return {row["name"] for row in rows} + + +def op_points(rows: list[dict], op: str) -> list[tuple[int, float]]: + grouped = defaultdict(list) + for row in rows: + if row["name"] == op: + grouped[row["N"]].append(row["cpu_time"]) + return sorted((n, median(values)) for n, values in grouped.items()) + + +def smooth(points: list[tuple[int, float]], window: int) -> list[float]: + if window <= 1: + return [time for _, time in points] + values = [time for _, time in points] + radius = window // 2 + smoothed = [] + for index in range(len(values)): + lo = max(0, index - radius) + hi = min(len(values), index + radius + 1) + smoothed.append(median(values[lo:hi])) + return smoothed def main(): @@ -129,9 +178,9 @@ def main(): if args.sdsl_json is not None: df_sdsl = load_bench_json(args.sdsl_json) - rmm_ops = set(df_rmm["name"].dropna().astype(str).unique()) + rmm_ops = ops_in(df_rmm) if df_sdsl is not None: - sdsl_ops = set(df_sdsl["name"].dropna().astype(str).unique()) + sdsl_ops = ops_in(df_sdsl) common_ops = rmm_ops & sdsl_ops ops_to_plot = [op for op in OPS_ORDER if op in common_ops] extra = sorted(common_ops - set(OPS_ORDER)) @@ -140,66 +189,52 @@ def main(): ops_to_plot = OPS_ORDER for op in ops_to_plot: - d_rmm = df_rmm[df_rmm["name"] == op].copy() - d_sdsl = None + d_rmm = op_points(df_rmm, op) + d_sdsl = [] if df_sdsl is not None: - d_sdsl = df_sdsl[df_sdsl["name"] == op].copy() - if d_rmm.empty or (df_sdsl is not None and (d_sdsl is None or d_sdsl.empty)): + d_sdsl = op_points(df_sdsl, op) + if not d_rmm or (df_sdsl is not None and not d_sdsl): continue plt.figure() - if not d_rmm.empty: - d_rmm = ( - d_rmm.groupby("N", as_index=False)["cpu_time"].median().sort_values("N") - ) - y_rmm = d_rmm["cpu_time"] - if args.smooth and args.smooth > 1: - d_rmm["ns_smooth"] = ( - d_rmm["cpu_time"] - .rolling(window=args.smooth, center=True, min_periods=1) - .median() - ) - y_rmm = d_rmm["ns_smooth"] + if d_rmm: + x_rmm = [n for n, _ in d_rmm] + raw_rmm = [time for _, time in d_rmm] + y_rmm = smooth(d_rmm, args.smooth) plt.scatter( - d_rmm["N"], - d_rmm["cpu_time"], - s=8, - alpha=0.3, + x_rmm, + raw_rmm, + s=18, + alpha=0.7, linewidths=0, ) plt.plot( - d_rmm["N"], + x_rmm, y_rmm, - linewidth=1.5, - label="RmMTree", + marker="o", + markersize=3, + linewidth=0.85, + label="Pixie RmMTree", ) - if d_sdsl is not None and not d_sdsl.empty: - d_sdsl = ( - d_sdsl.groupby("N", as_index=False)["cpu_time"] - .median() - .sort_values("N") - ) - y_sdsl = d_sdsl["cpu_time"] - if args.smooth and args.smooth > 1: - d_sdsl["ns_smooth"] = ( - d_sdsl["cpu_time"] - .rolling(window=args.smooth, center=True, min_periods=1) - .median() - ) - y_sdsl = d_sdsl["ns_smooth"] + if d_sdsl: + x_sdsl = [n for n, _ in d_sdsl] + raw_sdsl = [time for _, time in d_sdsl] + y_sdsl = smooth(d_sdsl, args.smooth) plt.scatter( - d_sdsl["N"], - d_sdsl["cpu_time"], - s=8, - alpha=0.3, - linewidths=0.9, + x_sdsl, + raw_sdsl, + s=18, + alpha=0.7, + linewidths=0.85, marker="x", ) plt.plot( - d_sdsl["N"], + x_sdsl, y_sdsl, - linewidth=1.5, + marker="x", + markersize=3, + linewidth=0.85, linestyle="--", label="sdsl-lite", ) @@ -208,14 +243,14 @@ def main(): plt.xscale("log", base=2) plt.xlabel("Sequence size N (bits)") plt.ylabel("Time per operation, ns") - if d_sdsl is not None and not d_sdsl.empty: - plt.title(f"RmMTree vs sdsl-lite: {op}") + if d_sdsl: + plt.title(f"RmMTree comparison: {op}") else: plt.title(f"RmMTree: {op}") plt.grid(True, which="both", linestyle="--", alpha=0.4) plt.legend(loc="best") out = os.path.join(args.save_dir, f"{op}.png") - plt.savefig(out, bbox_inches="tight", dpi=160) + plt.savefig(out, bbox_inches="tight", dpi=180) if not args.show: plt.close() print(f"[saved] {out}") diff --git a/scripts/plot_size_benchmarks.py b/scripts/plot_size_benchmarks.py new file mode 100644 index 0000000..d7fab7c --- /dev/null +++ b/scripts/plot_size_benchmarks.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python3 +"""Plot Google Benchmark timing against evaluated bitvector/tree sizes.""" + +from __future__ import annotations + +import argparse +import json +import os +import tempfile +from pathlib import Path +from typing import Dict, Iterable, List, Tuple + +os.environ.setdefault( + "MPLCONFIGDIR", str(Path(tempfile.gettempdir()) / "pixie-matplotlib") +) + +import matplotlib + +matplotlib.use("Agg") +import matplotlib.pyplot as plt + + +Point = Tuple[int, float] + + +def _clean_name(name: str) -> str: + for suffix in ("_mean", "_median", "_stddev", "_cv"): + if name.endswith(suffix): + return name[: -len(suffix)] + return name + + +def _is_timing_row(bench: dict) -> bool: + name = bench.get("name", "") + if any(name.endswith(suffix) for suffix in ("_stddev", "_cv")): + return False + run_type = bench.get("run_type") + if run_type == "aggregate": + return name.endswith("_mean") + return run_type in (None, "iteration") + + +def _guess_size_key(benchmarks: Iterable[dict]) -> str: + candidates = ("n", "N", "tree_size", "size") + for bench in benchmarks: + params = bench.get("params", {}) + for key in candidates: + if key in params: + return key + for key in bench.get("ArgNames", []): + if key in candidates: + return key + for key in candidates: + if key in bench: + return key + return "n" + + +def _extract_size_from_name(name: str, size_key: str) -> int: + marker = f"{size_key}:" + if marker not in name: + return 0 + try: + return int(name.split(marker, 1)[1].split("/", 1)[0]) + except ValueError: + return 0 + + +def _extract_size(bench: dict, size_key: str) -> int: + params = bench.get("params", {}) + if size_key in params: + return int(params[size_key]) + if size_key in bench: + return int(bench[size_key]) + args = bench.get("args", []) + if args: + return int(args[0]) + return _extract_size_from_name(bench.get("name", ""), size_key) + + +def _load_series(path: Path, size_key: str | None) -> Dict[str, List[Point]]: + data = json.loads(path.read_text(encoding="utf-8")) + benchmarks = data.get("benchmarks", data if isinstance(data, list) else []) + key = size_key or _guess_size_key(benchmarks) + series: Dict[str, List[Point]] = {} + for bench in benchmarks: + if not _is_timing_row(bench): + continue + name = _clean_name(bench.get("name", "unknown")).split("/", 1)[0] + size = _extract_size(bench, key) + cpu_time = float(bench.get("cpu_time", bench.get("real_time", 0.0))) + if size <= 0 or cpu_time <= 0: + continue + series.setdefault(name, []).append((size, cpu_time)) + + collapsed: Dict[str, List[Point]] = {} + for name, points in series.items(): + by_size: Dict[int, List[float]] = {} + for size, time in points: + by_size.setdefault(size, []).append(time) + collapsed[name] = sorted( + (size, sorted(times)[len(times) // 2]) for size, times in by_size.items() + ) + return collapsed + + +def _plot(series_by_label: Dict[str, List[Point]], output: Path, title: str) -> None: + plt.figure(figsize=(10, 6)) + for label, points in sorted(series_by_label.items()): + xs = [size for size, _ in points] + ys = [time for _, time in points] + plt.plot( + xs, + ys, + marker="o", + markersize=3, + linewidth=0.85, + label=label, + ) + plt.xscale("log", base=2) + plt.xlabel("Evaluated size") + plt.ylabel("CPU time (ns)") + plt.title(title) + plt.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.35) + plt.legend(fontsize=7) + plt.tight_layout() + output.parent.mkdir(parents=True, exist_ok=True) + plt.savefig(output, dpi=180) + + +def main() -> int: + parser = argparse.ArgumentParser( + description=( + "Plot bitvector/tree size benchmarks with x-log scale, thin lines, " + "and markers at evaluated sizes." + ) + ) + parser.add_argument("json", type=Path, nargs="+", help="Benchmark JSON file(s)") + parser.add_argument( + "-o", "--output", type=Path, default=Path("graphs/size_benchmarks.png") + ) + parser.add_argument("--title", default="Benchmark time vs size") + parser.add_argument("--size-key", help="Override size key, e.g. n, N, tree_size") + args = parser.parse_args() + + combined: Dict[str, List[Point]] = {} + multi_file = len(args.json) > 1 + for path in args.json: + for name, points in _load_series(path, args.size_key).items(): + label = f"{path.stem}:{name}" if multi_file else name + combined[label] = points + if not combined: + raise SystemExit("No valid benchmark series found.") + _plot(combined, args.output, args.title) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/benchmarks/bench_rmm.cpp b/src/benchmarks/bench_rmm.cpp index 5478b8e..38e0102 100644 --- a/src/benchmarks/bench_rmm.cpp +++ b/src/benchmarks/bench_rmm.cpp @@ -1,605 +1,8 @@ -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using pixie::RmMTree; - -struct Args { - int min_exp = 14; - int max_exp = 22; - std::size_t Q = 200000; - double p1 = 0.5; - std::uint64_t seed = 42; - std::size_t block_bits = 64; - int per_octave = 10; - std::vector explicit_sizes; -}; - -static Args args; - -static void parse_args_and_strip(int& argc, char**& argv) { - auto is_value_token = [&](int candidate_index) -> bool { - if (candidate_index >= argc) { - return false; - } - std::string candidate = argv[candidate_index]; - return !(candidate.size() >= 2 && candidate[0] == '-' && - candidate[1] == '-'); - }; - - auto get_value = [&](const std::string& key) -> std::optional { - std::string prefix = "--" + key; - std::string prefix_with_equals = prefix + "="; - for (int argument_index = 1; argument_index < argc; ++argument_index) { - std::string argument = argv[argument_index]; - if (argument.rfind(prefix_with_equals, 0) == 0) { - return argument.substr(prefix_with_equals.size()); - } - if (argument == prefix && is_value_token(argument_index + 1)) { - return std::string(argv[argument_index + 1]); - } - } - return std::nullopt; - }; - - auto strip = [&](const std::string& key) { - std::string prefix = "--" + key; - std::string prefix_with_equals = prefix + "="; - int write_index = 0; - for (int argument_index = 0; argument_index < argc; ++argument_index) { - std::string argument = argv[argument_index]; - if (argument_index > 0 && argument.rfind(prefix_with_equals, 0) == 0) { - continue; - } - if (argument_index > 0 && argument == prefix && - is_value_token(argument_index + 1)) { - ++argument_index; - continue; - } - argv[write_index++] = argv[argument_index]; - } - argc = write_index; - }; - - if (auto argument_value = get_value("min_exp")) { - args.min_exp = std::stoi(*argument_value); - strip("min_exp"); - } - if (auto argument_value = get_value("max_exp")) { - args.max_exp = std::stoi(*argument_value); - strip("max_exp"); - } - if (auto argument_value = get_value("Q")) { - args.Q = std::stoull(*argument_value); - strip("Q"); - } - if (auto argument_value = get_value("p1")) { - args.p1 = std::stod(*argument_value); - strip("p1"); - } - if (auto argument_value = get_value("seed")) { - args.seed = std::stoull(*argument_value); - strip("seed"); - } - if (auto argument_value = get_value("block_bits")) { - args.block_bits = std::stoull(*argument_value); - strip("block_bits"); - } - if (auto argument_value = get_value("per_octave")) { - args.per_octave = std::stoi(*argument_value); - strip("per_octave"); - } - - if (auto argument_value = get_value("explicit_sizes")) { - std::string sizes_text = *argument_value; - strip("explicit_sizes"); - std::size_t position = 0; - while (position < sizes_text.size()) { - while (position < sizes_text.size() && - (sizes_text[position] == ',' || - std::isspace(static_cast(sizes_text[position])))) { - ++position; - } - std::size_t start_index = position; - while (position < sizes_text.size() && - std::isdigit(static_cast(sizes_text[position]))) { - ++position; - } - if (start_index < position) { - args.explicit_sizes.push_back(std::stoull( - sizes_text.substr(start_index, position - start_index))); - } - } - std::sort(args.explicit_sizes.begin(), args.explicit_sizes.end()); - args.explicit_sizes.erase( - std::unique(args.explicit_sizes.begin(), args.explicit_sizes.end()), - args.explicit_sizes.end()); - } -} - -static std::string random_dyck_bits(std::size_t N, - double p1, - std::mt19937_64& rng) { - std::size_t pairs = N >> 1; - std::size_t L = pairs << 1; - std::string s; - s.resize(L); - if (L == 0) { - return s; - } - p1 = std::clamp(p1, 0., 1.); - std::size_t opens_left = pairs; - std::size_t closes_left = pairs; - int balance = 0; - std::uniform_real_distribution U(0., 1.); - for (std::size_t i = 0; i < L; ++i) { - if (opens_left == 0) { - s[i] = '0'; - --closes_left; - --balance; - continue; - } - - if (closes_left == 0) { - s[i] = '1'; - --opens_left; - ++balance; - continue; - } - - if (balance == 0) { - s[i] = '1'; - --opens_left; - ++balance; - continue; - } - if (U(rng) < p1) { - s[i] = '1'; - --opens_left; - ++balance; - } else { - s[i] = '0'; - --closes_left; - --balance; - } - } - return s; -} - -static std::vector build_size_grid() { - if (!args.explicit_sizes.empty()) { - std::vector v = args.explicit_sizes; - v.erase(std::remove(v.begin(), v.end(), 0), v.end()); - return v; - } - - const int lo = args.min_exp, hi = args.max_exp; - if (args.per_octave <= 0) { - std::vector v; - for (int e = lo; e <= hi; ++e) { - v.push_back(std::size_t(1) << e); - } - return v; - } - - int per_octave_steps = args.per_octave; - std::set N_cands; - std::size_t min_N = std::size_t(1) << lo, max_N = std::size_t(1) << hi; - - for (int e = lo; e <= hi; ++e) { - for (int t = 0; t <= per_octave_steps; ++t) { - long double x = e + (long double)t / (long double)per_octave_steps; - std::size_t N = (std::size_t)std::llround(std::pow(2.0L, x)); - if (N < min_N) { - N = min_N; - } - if (N > max_N) { - N = max_N; - } - if (N) { - N_cands.insert(N); - } - } - } - - std::vector v(N_cands.begin(), N_cands.end()); - if (v.front() != min_N) { - v.insert(v.begin(), min_N); - } - if (v.back() != max_N) { - v.push_back(max_N); - } - return v; -} - -struct Pools { - std::vector inds_any; - std::vector rank10_end_positions; - std::vector open_positions_zero_based; - std::vector open_positions_one_based; - std::vector close_positions_one_based; - std::vector inds; - std::vector inds_1N; - std::vector deltas; - std::vector> segs; - - std::vector ks1, ks0, ks10; - std::vector minselect_q; -}; - -struct Dataset { - std::size_t N{}; - std::string bits; - RmMTree t; - - std::size_t cnt1{}, cnt0{}, cnt10{}; - Pools pool; -}; - -static std::vector> keepalive; - -static std::size_t count10(const std::string& s) { - std::size_t c = 0; - if (s.size() < 2) { - return 0; - } - for (std::size_t i = 0; i + 1 < s.size(); ++i) { - if (s[i] == '1' && s[i + 1] == '0') { - ++c; - } - } - return c; -} - -static Dataset build_dataset(std::size_t N) { - std::mt19937_64 rng(args.seed ^ - static_cast(N) * 0x9E3779B185EBCA87ull); - - Dataset d; - d.bits = random_dyck_bits(N, args.p1, rng); - d.N = d.bits.size(); - N = d.N; - - if (args.block_bits == 0) { - d.t = RmMTree(d.bits, 0, -1.0f); - } else { - d.t = RmMTree(d.bits, args.block_bits); - } - - d.cnt1 = std::count(d.bits.begin(), d.bits.end(), '1'); - d.cnt0 = N - d.cnt1; - d.cnt10 = count10(d.bits); - - std::vector open_positions_zero_based; - std::vector open_positions_one_based; - std::vector close_positions_one_based; - open_positions_zero_based.reserve(d.N >> 1); - open_positions_one_based.reserve(d.N >> 1); - close_positions_one_based.reserve(d.N >> 1); - for (std::size_t i = 0; i < N; ++i) { - if (d.bits[i] == '1') { - open_positions_zero_based.push_back(i); - open_positions_one_based.push_back(i + 1); - } else { - close_positions_one_based.push_back(i + 1); - } - } - - std::size_t sample_count = std::max( - 1, std::min(args.Q, std::size_t(32768))); - - std::uniform_int_distribution nonedge_index_distribution( - N > 1 ? 1 : 0, N > 1 ? (N - 1) : 0); - std::uniform_int_distribution ind_dist(0, N ? (N - 1) : 0); - std::uniform_int_distribution ind_dist_1N(1, N ? N : 1); - std::uniform_int_distribution rank10_end_position_distribution( - N >= 2 ? 2 : 0, N >= 2 ? N : 0); - std::uniform_int_distribution delta_distribution(-8, +8); - - d.pool.inds_any.resize(sample_count); - d.pool.rank10_end_positions.resize(sample_count); - d.pool.inds.resize(sample_count); - d.pool.inds_1N.resize(sample_count); - d.pool.deltas.resize(sample_count); - d.pool.segs.resize(sample_count); - - auto rand_ij = [&]() -> std::pair { - if (N == 0) { - return {0, 0}; - } - std::size_t a = ind_dist(rng), b = ind_dist(rng); - if (a > b) { - std::swap(a, b); - } - return {a, b}; - }; - - for (std::size_t i = 0; i < sample_count; ++i) { - d.pool.inds_any[i] = nonedge_index_distribution(rng); - d.pool.rank10_end_positions[i] = rank10_end_position_distribution(rng); - d.pool.inds[i] = (N ? ind_dist(rng) : 0); - d.pool.inds_1N[i] = (N ? ind_dist_1N(rng) : 0); - d.pool.deltas[i] = delta_distribution(rng); - d.pool.segs[i] = rand_ij(); - } - - auto fill_from_candidates = [&](const std::vector& candidates, - std::vector& destination, - std::size_t fallback_value) { - destination.resize(sample_count); - if (candidates.empty()) { - std::fill(destination.begin(), destination.end(), fallback_value); - return; - } - std::uniform_int_distribution index_distribution( - 0, candidates.size() - 1); - for (std::size_t sample_index = 0; sample_index < sample_count; - ++sample_index) { - destination[sample_index] = candidates[index_distribution(rng)]; - } - }; - const std::size_t one_based_fallback = (N > 0 ? 1 : 0); - fill_from_candidates(open_positions_zero_based, - d.pool.open_positions_zero_based, 0); - fill_from_candidates(open_positions_one_based, - d.pool.open_positions_one_based, one_based_fallback); - fill_from_candidates(close_positions_one_based, - d.pool.close_positions_one_based, one_based_fallback); - - auto fill_ks = [&](std::size_t total, std::vector& out) { - out.resize(sample_count); - if (total == 0) { - std::fill(out.begin(), out.end(), 1); - return; - } - std::uniform_int_distribution dist(1, total); - for (std::size_t i = 0; i < sample_count; ++i) { - out[i] = dist(rng); - } - }; - fill_ks(d.cnt1, d.pool.ks1); - fill_ks(d.cnt0, d.pool.ks0); - fill_ks(d.cnt10, d.pool.ks10); - - d.pool.minselect_q.resize(sample_count); - for (std::size_t i = 0; i < sample_count; ++i) { - auto [l, r] = d.pool.segs[i]; - std::size_t c = d.t.mincount(l, r); - if (c == 0) { - c = 1; - } - std::uniform_int_distribution uq(1, c); - d.pool.minselect_q[i] = uq(rng); - } - - return d; -} - -template -static void register_op(const std::string& op, - std::shared_ptr data, - Fn&& body) { - auto idx_ptr = std::make_shared(0); - - auto* b = benchmark::RegisterBenchmark( - op.c_str(), [data, idx_ptr, body](benchmark::State& state) { - const Dataset& D = *data; - for (auto _ : state) { - std::size_t i = (*idx_ptr)++; - body(state, D, i); - } - state.counters["N"] = static_cast(D.N); - state.counters["seed"] = static_cast(args.seed); - state.counters["block_bits"] = static_cast(args.block_bits); - }); - - b->Unit(benchmark::kNanosecond); -} - -static void register_all() { - auto Ns = build_size_grid(); - for (std::size_t N : Ns) { - auto data = std::make_shared(build_dataset(N)); - keepalive.push_back(data); - - const auto& P = data->pool; - - register_op("rank1", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.rank1(P.inds_any[k % P.inds_any.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("rank0", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.rank0(P.inds_any[k % P.inds_any.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("select1", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.select1(P.ks1[k % P.ks1.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("select0", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.select0(P.ks0[k % P.ks0.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op( - "rank10", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - benchmark::DoNotOptimize(dataset.t.rank10( - dataset.pool.rank10_end_positions - [sample_index % dataset.pool.rank10_end_positions.size()])); - }); - - register_op("select10", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.select10(P.ks10[k % P.ks10.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("excess", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.excess(P.inds_any[k % P.inds_any.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("fwdsearch", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.fwdsearch(P.inds[k % P.inds.size()], - P.deltas[k % P.deltas.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("bwdsearch", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto r = D.t.bwdsearch(P.inds_1N[k % P.inds_1N.size()], - P.deltas[k % P.deltas.size()]); - benchmark::DoNotOptimize(r); - }); - - register_op("range_min_query_pos", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - const std::size_t segment_index = k % P.segs.size(); - auto [range_begin, range_end] = P.segs[segment_index]; - auto min_position = - D.t.range_min_query_pos(range_begin, range_end); - benchmark::DoNotOptimize(min_position); - }); - - register_op("range_min_query_val", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - const std::size_t segment_index = k % P.segs.size(); - auto [range_begin, range_end] = P.segs[segment_index]; - auto min_value = - D.t.range_min_query_val(range_begin, range_end); - benchmark::DoNotOptimize(min_value); - }); - - register_op("range_max_query_pos", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto [i, j] = P.segs[k % P.segs.size()]; - auto r = D.t.range_max_query_pos(i, j); - benchmark::DoNotOptimize(r); - }); - - register_op("range_max_query_val", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto [i, j] = P.segs[k % P.segs.size()]; - auto r = D.t.range_max_query_val(i, j); - benchmark::DoNotOptimize(r); - }); - - register_op("mincount", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - auto [i, j] = P.segs[k % P.segs.size()]; - auto r = D.t.mincount(i, j); - benchmark::DoNotOptimize(r); - }); - - register_op("minselect", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - const std::size_t idx = k % P.segs.size(); - auto [i, j] = P.segs[idx]; - auto q = P.minselect_q[idx]; - auto r = D.t.minselect(i, j, q); - benchmark::DoNotOptimize(r); - }); - - register_op("close", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - if (dataset.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t open_position = - dataset.pool.open_positions_zero_based - [sample_index % - dataset.pool.open_positions_zero_based.size()]; - benchmark::DoNotOptimize(dataset.t.close(open_position)); - }); - - register_op( - "open", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - const std::size_t close_position_one_based = - dataset.pool.close_positions_one_based - [sample_index % - dataset.pool.close_positions_one_based.size()]; - benchmark::DoNotOptimize(dataset.t.open(close_position_one_based)); - }); - - register_op( - "enclose", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - const std::size_t open_position_one_based = - dataset.pool.open_positions_one_based - [sample_index % dataset.pool.open_positions_one_based.size()]; - benchmark::DoNotOptimize(dataset.t.enclose(open_position_one_based)); - }); - } -} +#include "rmm_benchmark_base.h" int main(int argc, char** argv) { - benchmark::MaybeReenterWithoutASLR(argc, argv); - parse_args_and_strip(argc, argv); - - auto has = [&](const char* key) { - std::string p1 = std::string(key) + "="; - for (int i = 1; i < argc; ++i) { - std::string s = argv[i]; - if (s == key || s.rfind(p1, 0) == 0) { - return true; - } - } - return false; - }; - - static std::vector extra; - if (!has("--benchmark_out_format")) { - extra.emplace_back("--benchmark_out_format=json"); - } - if (!has("--benchmark_counters_tabular")) { - extra.emplace_back("--benchmark_counters_tabular=true"); - } - if (!has("--benchmark_time_unit")) { - extra.emplace_back("--benchmark_time_unit=ns"); - } - - static std::vector argv_vec; - argv_vec.assign(argv, argv + argc); - for (auto& s : extra) { - argv_vec.push_back(s.data()); - } - argv_vec.push_back(nullptr); - argc = (int)argv_vec.size() - 1; - argv = argv_vec.data(); - - benchmark::Initialize(&argc, argv); - register_all(); - benchmark::RunSpecifiedBenchmarks(); - benchmark::Shutdown(); - return 0; + pixie_bench::RmMBenchmark benchmark; + return benchmark.Run(argc, argv); } diff --git a/src/benchmarks/bench_rmm_sdsl.cpp b/src/benchmarks/bench_rmm_sdsl.cpp index d7e450e..0dd9eac 100644 --- a/src/benchmarks/bench_rmm_sdsl.cpp +++ b/src/benchmarks/bench_rmm_sdsl.cpp @@ -1,594 +1,26 @@ -#include +#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include -using sdsl::bit_vector; -using sdsl::bp_support_sada; -using BpSupport = bp_support_sada<>; +#include "rmm_benchmark_base.h" -struct Args { - int min_exp = 14; - int max_exp = 22; - std::size_t Q = 200000; - double p1 = 0.5; - std::uint64_t seed = 42; - std::size_t block_bits = 64; - int per_octave = 10; - std::vector explicit_sizes; -}; - -static Args args; - -static void parse_args_and_strip(int& argc, char**& argv) { - auto is_value_token = [&](int candidate_index) -> bool { - if (candidate_index >= argc) { - return false; - } - std::string candidate = argv[candidate_index]; - return !(candidate.size() >= 2 && candidate[0] == '-' && - candidate[1] == '-'); - }; +namespace pixie_bench { - auto get_value = [&](const std::string& key) -> std::optional { - std::string prefix = "--" + key; - std::string prefix_with_equals = prefix + "="; - for (int argument_index = 1; argument_index < argc; ++argument_index) { - std::string argument = argv[argument_index]; - if (argument.rfind(prefix_with_equals, 0) == 0) { - return argument.substr(prefix_with_equals.size()); - } - if (argument == prefix && is_value_token(argument_index + 1)) { - return std::string(argv[argument_index + 1]); - } - } - return std::nullopt; - }; - - auto strip = [&](const std::string& key) { - std::string prefix = "--" + key; - std::string prefix_with_equals = prefix + "="; - int write_index = 0; - for (int argument_index = 0; argument_index < argc; ++argument_index) { - std::string argument = argv[argument_index]; - if (argument_index > 0 && argument.rfind(prefix_with_equals, 0) == 0) { - continue; - } - if (argument_index > 0 && argument == prefix && - is_value_token(argument_index + 1)) { - ++argument_index; - continue; - } - argv[write_index++] = argv[argument_index]; - } - argc = write_index; - }; - - if (auto argument_value = get_value("min_exp")) { - args.min_exp = std::stoi(*argument_value); - strip("min_exp"); - } - if (auto argument_value = get_value("max_exp")) { - args.max_exp = std::stoi(*argument_value); - strip("max_exp"); - } - if (auto argument_value = get_value("Q")) { - args.Q = std::stoull(*argument_value); - strip("Q"); - } - if (auto argument_value = get_value("p1")) { - args.p1 = std::stod(*argument_value); - strip("p1"); - } - if (auto argument_value = get_value("seed")) { - args.seed = std::stoull(*argument_value); - strip("seed"); - } - if (auto argument_value = get_value("block_bits")) { - args.block_bits = std::stoull(*argument_value); - strip("block_bits"); - } - if (auto argument_value = get_value("per_octave")) { - args.per_octave = std::stoi(*argument_value); - strip("per_octave"); - } - - if (auto argument_value = get_value("explicit_sizes")) { - std::string sizes_text = *argument_value; - strip("explicit_sizes"); - std::size_t position = 0; - while (position < sizes_text.size()) { - while (position < sizes_text.size() && - (sizes_text[position] == ',' || - std::isspace(static_cast(sizes_text[position])))) { - ++position; - } - std::size_t start_index = position; - while (position < sizes_text.size() && - std::isdigit(static_cast(sizes_text[position]))) { - ++position; - } - if (start_index < position) { - args.explicit_sizes.push_back(std::stoull( - sizes_text.substr(start_index, position - start_index))); - } - } - std::sort(args.explicit_sizes.begin(), args.explicit_sizes.end()); - args.explicit_sizes.erase( - std::unique(args.explicit_sizes.begin(), args.explicit_sizes.end()), - args.explicit_sizes.end()); - } -} - -static std::string random_dyck_bits(std::size_t N, - double p1, - std::mt19937_64& rng) { - std::size_t pairs = N >> 1; - std::size_t L = pairs << 1; - std::string s; - s.resize(L); - if (L == 0) { - return s; - } - p1 = std::clamp(p1, 0., 1.); - std::size_t opens_left = pairs; - std::size_t closes_left = pairs; - int balance = 0; - std::uniform_real_distribution U(0., 1.); - for (std::size_t i = 0; i < L; ++i) { - if (opens_left == 0) { - s[i] = '0'; - --closes_left; - --balance; - continue; - } - - if (closes_left == 0) { - s[i] = '1'; - --opens_left; - ++balance; - continue; - } - - if (balance == 0) { - s[i] = '1'; - --opens_left; - ++balance; - continue; - } - if (U(rng) < p1) { - s[i] = '1'; - --opens_left; - ++balance; - } else { - s[i] = '0'; - --closes_left; - --balance; - } - } - return s; -} +template <> +struct RmMBenchmarkTraits { + static constexpr std::size_t DefaultBlockBits = 64; -static std::vector build_size_grid() { - if (!args.explicit_sizes.empty()) { - std::vector v = args.explicit_sizes; - v.erase(std::remove(v.begin(), v.end(), 0), v.end()); - return v; + static bool SupportsOp(std::string_view op) { + return op == "rank1" || op == "rank0" || op == "select1" || + op == "excess" || op == "fwdsearch" || op == "range_min_query_pos" || + op == "range_min_query_val" || op == "close" || op == "open" || + op == "enclose"; } - - const int lo = args.min_exp, hi = args.max_exp; - if (args.per_octave <= 0) { - std::vector v; - for (int e = lo; e <= hi; ++e) { - v.push_back(std::size_t(1) << e); - } - return v; - } - - int per_octave_steps = args.per_octave; - std::set N_cands; - std::size_t min_N = std::size_t(1) << lo, max_N = std::size_t(1) << hi; - - for (int e = lo; e <= hi; ++e) { - for (int t = 0; t <= per_octave_steps; ++t) { - long double x = e + (long double)t / (long double)per_octave_steps; - std::size_t N = (std::size_t)std::llround(std::pow(2.0L, x)); - if (N < min_N) { - N = min_N; - } - if (N > max_N) { - N = max_N; - } - if (N) { - N_cands.insert(N); - } - } - } - - std::vector v(N_cands.begin(), N_cands.end()); - if (v.front() != min_N) { - v.insert(v.begin(), min_N); - } - if (v.back() != max_N) { - v.push_back(max_N); - } - return v; -} - -struct Pools { - std::vector inds_any; - std::vector rank10_end_positions; - std::vector open_positions_zero_based; - std::vector open_positions_one_based; - std::vector close_positions_one_based; - std::vector inds; - std::vector inds_1N; - std::vector deltas; - std::vector> segs; - - std::vector ks1, ks0, ks10; - std::vector minselect_q; }; -struct Dataset { - std::size_t N{}; - std::string bits; - bit_vector bv; - BpSupport t; - - std::size_t cnt1{}, cnt0{}, cnt10{}; - Pools pool; -}; - -static std::vector> keepalive; - -static std::size_t count10(const std::string& s) { - std::size_t c = 0; - if (s.size() < 2) { - return 0; - } - for (std::size_t i = 0; i + 1 < s.size(); ++i) { - if (s[i] == '1' && s[i + 1] == '0') { - ++c; - } - } - return c; -} - -static std::shared_ptr build_dataset(std::size_t N) { - std::mt19937_64 rng(args.seed ^ - static_cast(N) * 0x9E3779B185EBCA87ull); - - auto d = std::make_shared(); - d->bits = random_dyck_bits(N, args.p1, rng); - d->N = d->bits.size(); - N = d->N; - - d->bv = bit_vector(N); - for (std::size_t i = 0; i < N; ++i) { - d->bv[i] = (d->bits[i] == '1'); - } - d->t = BpSupport(&d->bv); - - d->cnt1 = std::count(d->bits.begin(), d->bits.end(), '1'); - d->cnt0 = N - d->cnt1; - d->cnt10 = count10(d->bits); - - std::vector open_positions_zero_based; - std::vector open_positions_one_based; - std::vector close_positions_one_based; - open_positions_zero_based.reserve(d->N >> 1); - open_positions_one_based.reserve(d->N >> 1); - close_positions_one_based.reserve(d->N >> 1); - for (std::size_t i = 0; i < N; ++i) { - if (d->bits[i] == '1') { - open_positions_zero_based.push_back(i); - open_positions_one_based.push_back(i + 1); - } else { - close_positions_one_based.push_back(i + 1); - } - } - - std::size_t sample_count = std::max( - 1, std::min(args.Q, std::size_t(32768))); - - std::uniform_int_distribution nonedge_index_distribution( - N > 1 ? 1 : 0, N > 1 ? (N - 1) : 0); - std::uniform_int_distribution ind_dist(0, N ? (N - 1) : 0); - std::uniform_int_distribution ind_dist_1N(1, N ? N : 1); - std::uniform_int_distribution rank10_end_position_distribution( - N >= 2 ? 2 : 0, N >= 2 ? N : 0); - std::uniform_int_distribution delta_distribution(-8, +8); - - d->pool.inds_any.resize(sample_count); - d->pool.rank10_end_positions.resize(sample_count); - d->pool.inds.resize(sample_count); - d->pool.inds_1N.resize(sample_count); - d->pool.deltas.resize(sample_count); - d->pool.segs.resize(sample_count); - - auto rand_ij = [&]() -> std::pair { - if (N == 0) { - return {0, 0}; - } - std::size_t a = ind_dist(rng), b = ind_dist(rng); - if (a > b) { - std::swap(a, b); - } - return {a, b}; - }; - - for (std::size_t i = 0; i < sample_count; ++i) { - d->pool.inds_any[i] = nonedge_index_distribution(rng); - d->pool.rank10_end_positions[i] = rank10_end_position_distribution(rng); - d->pool.inds[i] = (N ? ind_dist(rng) : 0); - d->pool.inds_1N[i] = (N ? ind_dist_1N(rng) : 0); - d->pool.deltas[i] = delta_distribution(rng); - d->pool.segs[i] = rand_ij(); - } - - auto fill_from_candidates = [&](const std::vector& candidates, - std::vector& destination, - std::size_t fallback_value) { - destination.resize(sample_count); - if (candidates.empty()) { - std::fill(destination.begin(), destination.end(), fallback_value); - return; - } - std::uniform_int_distribution index_distribution( - 0, candidates.size() - 1); - for (std::size_t sample_index = 0; sample_index < sample_count; - ++sample_index) { - destination[sample_index] = candidates[index_distribution(rng)]; - } - }; - const std::size_t one_based_fallback = (N > 0 ? 1 : 0); - fill_from_candidates(open_positions_zero_based, - d->pool.open_positions_zero_based, 0); - fill_from_candidates(open_positions_one_based, - d->pool.open_positions_one_based, one_based_fallback); - fill_from_candidates(close_positions_one_based, - d->pool.close_positions_one_based, one_based_fallback); - - auto fill_ks = [&](std::size_t total, std::vector& out) { - out.resize(sample_count); - if (total == 0) { - std::fill(out.begin(), out.end(), 1); - return; - } - std::uniform_int_distribution dist(1, total); - for (std::size_t i = 0; i < sample_count; ++i) { - out[i] = dist(rng); - } - }; - fill_ks(d->cnt1, d->pool.ks1); - fill_ks(d->cnt0, d->pool.ks0); - fill_ks(d->cnt10, d->pool.ks10); - - d->pool.minselect_q.clear(); - - return d; -} - -template -static void register_op(const std::string& op, - std::shared_ptr data, - Fn&& body) { - auto idx_ptr = std::make_shared(0); - - auto* b = benchmark::RegisterBenchmark( - op.c_str(), [data, idx_ptr, body](benchmark::State& state) { - const Dataset& D = *data; - for (auto _ : state) { - std::size_t i = (*idx_ptr)++; - body(state, D, i); - } - state.counters["N"] = static_cast(D.N); - state.counters["seed"] = static_cast(args.seed); - state.counters["block_bits"] = static_cast(args.block_bits); - }); - - b->Unit(benchmark::kNanosecond); -} - -static void register_all() { - auto Ns = build_size_grid(); - for (std::size_t N : Ns) { - auto data = build_dataset(N); - keepalive.push_back(data); - - const auto& P = data->pool; - - register_op("rank1", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - std::size_t pos = P.inds_any[k % P.inds_any.size()]; - std::size_t r = 0; - if (pos > 0) { - r = D.t.rank(pos - 1); - } - benchmark::DoNotOptimize(r); - }); - - register_op("rank0", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - std::size_t pos = P.inds_any[k % P.inds_any.size()]; - std::size_t r = 0; - if (pos == 0) { - r = 0; - } else if (pos >= D.N) { - r = D.cnt0; - } else { - r = D.t.preceding_closing_parentheses(pos); - } - benchmark::DoNotOptimize(r); - }); - - register_op("select1", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.cnt1 == 0) { - benchmark::DoNotOptimize(0); - return; - } - std::size_t q = P.ks1[k % P.ks1.size()]; - auto r = D.t.select(q); - benchmark::DoNotOptimize(r); - }); - - register_op("excess", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - std::size_t pos = P.inds_any[k % P.inds_any.size()]; - BpSupport::difference_type r = 0; - if (pos > 0) { - r = D.t.excess(pos - 1); - } - benchmark::DoNotOptimize(r); - }); - - register_op("range_min_query_pos", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t segment_index = k % P.segs.size(); - auto [range_begin, range_end] = P.segs[segment_index]; - auto min_position = D.t.rmq(range_begin, range_end); - benchmark::DoNotOptimize(min_position); - }); - - register_op("range_min_query_val", data, - [&](benchmark::State&, const Dataset& D, std::size_t k) { - if (D.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t segment_index = k % P.segs.size(); - auto [range_begin, range_end] = P.segs[segment_index]; - auto min_position = D.t.rmq(range_begin, range_end); - BpSupport::difference_type min_value = 0; - if (min_position < D.N) { - const auto base_excess = - (range_begin == 0 ? 0 : D.t.excess(range_begin - 1)); - min_value = D.t.excess(min_position) - base_excess; - } - benchmark::DoNotOptimize(min_value); - }); - - register_op("close", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - if (dataset.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t open_position = - dataset.pool.open_positions_zero_based - [sample_index % - dataset.pool.open_positions_zero_based.size()]; - benchmark::DoNotOptimize(dataset.t.find_close(open_position)); - }); - - register_op( - "open", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - if (dataset.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t close_position_one_based = - dataset.pool.close_positions_one_based - [sample_index % - dataset.pool.close_positions_one_based.size()]; - const std::size_t close_position_zero_based = - (close_position_one_based > 0 ? close_position_one_based - 1 : 0); - benchmark::DoNotOptimize( - dataset.t.find_open(close_position_zero_based)); - }); - - register_op( - "enclose", data, - [&](benchmark::State&, const Dataset& dataset, - std::size_t sample_index) { - if (dataset.N == 0) { - benchmark::DoNotOptimize(0); - return; - } - const std::size_t open_position_one_based = - dataset.pool.open_positions_one_based - [sample_index % dataset.pool.open_positions_one_based.size()]; - const std::size_t open_position_zero_based = - (open_position_one_based > 0 ? open_position_one_based - 1 : 0); - benchmark::DoNotOptimize(dataset.t.enclose(open_position_zero_based)); - }); - } -} +} // namespace pixie_bench int main(int argc, char** argv) { - benchmark::MaybeReenterWithoutASLR(argc, argv); - parse_args_and_strip(argc, argv); - - auto has = [&](const char* key) { - std::string p1 = std::string(key) + "="; - for (int i = 1; i < argc; ++i) { - std::string s = argv[i]; - if (s == key || s.rfind(p1, 0) == 0) { - return true; - } - } - return false; - }; - - static std::vector extra; - if (!has("--benchmark_out_format")) { - extra.emplace_back("--benchmark_out_format=json"); - } - if (!has("--benchmark_counters_tabular")) { - extra.emplace_back("--benchmark_counters_tabular=true"); - } - if (!has("--benchmark_time_unit")) { - extra.emplace_back("--benchmark_time_unit=ns"); - } - - static std::vector argv_vec; - argv_vec.assign(argv, argv + argc); - for (auto& s : extra) { - argv_vec.push_back(s.data()); - } - argv_vec.push_back(nullptr); - argc = (int)argv_vec.size() - 1; - argv = argv_vec.data(); - - benchmark::Initialize(&argc, argv); - register_all(); - benchmark::RunSpecifiedBenchmarks(); - benchmark::Shutdown(); - return 0; + pixie_bench::RmMBenchmark benchmark; + return benchmark.Run(argc, argv); } diff --git a/src/benchmarks/rmm_benchmark_base.h b/src/benchmarks/rmm_benchmark_base.h new file mode 100644 index 0000000..6f2b601 --- /dev/null +++ b/src/benchmarks/rmm_benchmark_base.h @@ -0,0 +1,788 @@ +#pragma once + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace pixie_bench { + +struct RmMBenchmarkArgs { + int min_exp = 14; + int max_exp = 22; + std::size_t Q = 200000; + double p1 = 0.5; + std::uint64_t seed = 42; + std::size_t block_bits = 64; + bool block_bits_specified = false; + int per_octave = 0; + std::vector explicit_sizes; + std::set ops; + std::string benchmark_filter; +}; + +struct RmMBenchmarkPools { + std::vector inds_any; + std::vector rank10_end_positions; + std::vector open_positions_zero_based; + std::vector open_positions_one_based; + std::vector close_positions_zero_based; + std::vector close_positions_one_based; + std::vector inds; + std::vector inds_1N; + std::vector deltas; + std::vector> segs; + + std::vector ks1; + std::vector ks0; + std::vector ks10; + std::vector minselect_q; +}; + +inline const std::vector& AllRmMOps() { + static const std::vector ops = { + "rank1", + "rank0", + "select1", + "select0", + "rank10", + "select10", + "excess", + "fwdsearch", + "bwdsearch", + "range_min_query_pos", + "range_min_query_val", + "range_max_query_pos", + "range_max_query_val", + "mincount", + "minselect", + "close", + "open", + "enclose", + }; + return ops; +} + +inline std::vector SplitCsv(const std::string& text) { + std::vector out; + std::size_t position = 0; + while (position < text.size()) { + while (position < text.size() && + (text[position] == ',' || + std::isspace(static_cast(text[position])))) { + ++position; + } + const std::size_t start = position; + while (position < text.size() && text[position] != ',') { + ++position; + } + std::string value = text.substr(start, position - start); + while (!value.empty() && + std::isspace(static_cast(value.back()))) { + value.pop_back(); + } + if (!value.empty()) { + out.push_back(value); + } + } + return out; +} + +inline std::string RandomDyckBits(std::size_t N, + double p1, + std::mt19937_64& rng) { + const std::size_t pairs = N >> 1; + const std::size_t L = pairs << 1; + std::string s; + s.resize(L); + if (L == 0) { + return s; + } + p1 = std::clamp(p1, 0., 1.); + std::size_t opens_left = pairs; + std::size_t closes_left = pairs; + int balance = 0; + std::uniform_real_distribution U(0., 1.); + for (std::size_t i = 0; i < L; ++i) { + if (opens_left == 0) { + s[i] = '0'; + --closes_left; + --balance; + continue; + } + + if (closes_left == 0) { + s[i] = '1'; + --opens_left; + ++balance; + continue; + } + + if (balance == 0 || U(rng) < p1) { + s[i] = '1'; + --opens_left; + ++balance; + } else { + s[i] = '0'; + --closes_left; + --balance; + } + } + return s; +} + +inline std::size_t Count10(const std::string& s) { + std::size_t count = 0; + if (s.size() < 2) { + return 0; + } + for (std::size_t i = 0; i + 1 < s.size(); ++i) { + if (s[i] == '1' && s[i + 1] == '0') { + ++count; + } + } + return count; +} + +inline std::vector PackWordsLsbFirst(const std::string& bits) { + std::vector words((bits.size() + 63) / 64, 0); + for (std::size_t i = 0; i < bits.size(); ++i) { + if (bits[i] == '1') { + words[i >> 6] |= (std::uint64_t(1) << (i & 63)); + } + } + return words; +} + +template +struct RmMBenchmarkTraits { + static constexpr std::size_t DefaultBlockBits = 64; + + static bool SupportsOp(std::string_view) { return true; } +}; + +template +struct RmMBenchmarkDataset { + std::size_t N{}; + std::string bits; + std::vector words; + Tree tree; + + std::size_t cnt1{}; + std::size_t cnt0{}; + std::size_t cnt10{}; + RmMBenchmarkPools pool; +}; + +template +class RmMBenchmark { + public: + int Run(int argc, char** argv) { + benchmark::MaybeReenterWithoutASLR(argc, argv); + ParseArgsAndStrip(argc, argv); + AddDefaultBenchmarkArgs(argc, argv); + + benchmark::Initialize(&argc, argv); + RegisterAll(); + benchmark::RunSpecifiedBenchmarks(); + benchmark::Shutdown(); + return 0; + } + + const RmMBenchmarkArgs& args() const { return args_; } + + protected: + bool OpEnabled(std::string_view op) const { + const std::string op_string(op); + if (!args_.ops.empty()) { + return args_.ops.contains(op_string); + } + if (!args_.benchmark_filter.empty() && args_.benchmark_filter != "all" && + args_.benchmark_filter != ".*") { + try { + return std::regex_search(op_string, std::regex(args_.benchmark_filter)); + } catch (const std::regex_error&) { + return true; + } + } + return true; + } + + private: + using Dataset = RmMBenchmarkDataset; + + bool SupportsOp(std::string_view op) const { + return RmMBenchmarkTraits::SupportsOp(op); + } + + bool ActiveOp(std::string_view op) const { + return OpEnabled(op) && SupportsOp(op); + } + + void ParseArgsAndStrip(int& argc, char**& argv) { + auto is_value_token = [&](int candidate_index) -> bool { + if (candidate_index >= argc) { + return false; + } + std::string candidate = argv[candidate_index]; + return !(candidate.size() >= 2 && candidate[0] == '-' && + candidate[1] == '-'); + }; + + auto get_value = [&](const std::string& key) -> std::optional { + const std::string prefix = "--" + key; + const std::string prefix_with_equals = prefix + "="; + for (int argument_index = 1; argument_index < argc; ++argument_index) { + std::string argument = argv[argument_index]; + if (argument.rfind(prefix_with_equals, 0) == 0) { + return argument.substr(prefix_with_equals.size()); + } + if (argument == prefix && is_value_token(argument_index + 1)) { + return std::string(argv[argument_index + 1]); + } + } + return std::nullopt; + }; + + auto strip = [&](const std::string& key) { + const std::string prefix = "--" + key; + const std::string prefix_with_equals = prefix + "="; + int write_index = 0; + for (int argument_index = 0; argument_index < argc; ++argument_index) { + std::string argument = argv[argument_index]; + if (argument_index > 0 && argument.rfind(prefix_with_equals, 0) == 0) { + continue; + } + if (argument_index > 0 && argument == prefix && + is_value_token(argument_index + 1)) { + ++argument_index; + continue; + } + argv[write_index++] = argv[argument_index]; + } + argc = write_index; + }; + + if (auto argument_value = get_value("min_exp")) { + args_.min_exp = std::stoi(*argument_value); + strip("min_exp"); + } + if (auto argument_value = get_value("max_exp")) { + args_.max_exp = std::stoi(*argument_value); + strip("max_exp"); + } + if (auto argument_value = get_value("Q")) { + args_.Q = std::stoull(*argument_value); + strip("Q"); + } + if (auto argument_value = get_value("p1")) { + args_.p1 = std::stod(*argument_value); + strip("p1"); + } + if (auto argument_value = get_value("seed")) { + args_.seed = std::stoull(*argument_value); + strip("seed"); + } + if (auto argument_value = get_value("block_bits")) { + args_.block_bits = std::stoull(*argument_value); + args_.block_bits_specified = true; + strip("block_bits"); + } + if (auto argument_value = get_value("per_octave")) { + args_.per_octave = std::stoi(*argument_value); + strip("per_octave"); + } + if (auto argument_value = get_value("ops")) { + for (const auto& op : SplitCsv(*argument_value)) { + args_.ops.insert(op); + } + strip("ops"); + } + if (auto argument_value = get_value("benchmark_filter")) { + args_.benchmark_filter = *argument_value; + } + + if (auto argument_value = get_value("explicit_sizes")) { + strip("explicit_sizes"); + for (const auto& size_text : SplitCsv(*argument_value)) { + if (std::all_of(size_text.begin(), size_text.end(), [](char ch) { + return std::isdigit(static_cast(ch)); + })) { + args_.explicit_sizes.push_back(std::stoull(size_text)); + } + } + std::sort(args_.explicit_sizes.begin(), args_.explicit_sizes.end()); + args_.explicit_sizes.erase( + std::unique(args_.explicit_sizes.begin(), args_.explicit_sizes.end()), + args_.explicit_sizes.end()); + } + } + + void AddDefaultBenchmarkArgs(int& argc, char**& argv) { + auto has = [&](const char* key) { + const std::string prefix = std::string(key) + "="; + for (int i = 1; i < argc; ++i) { + const std::string argument = argv[i]; + if (argument == key || argument.rfind(prefix, 0) == 0) { + return true; + } + } + return false; + }; + + if (!has("--benchmark_out_format")) { + extra_args_.emplace_back("--benchmark_out_format=json"); + } + if (!has("--benchmark_counters_tabular")) { + extra_args_.emplace_back("--benchmark_counters_tabular=true"); + } + if (!has("--benchmark_time_unit")) { + extra_args_.emplace_back("--benchmark_time_unit=ns"); + } + + argv_vec_.assign(argv, argv + argc); + for (auto& argument : extra_args_) { + argv_vec_.push_back(argument.data()); + } + argv_vec_.push_back(nullptr); + argc = static_cast(argv_vec_.size()) - 1; + argv = argv_vec_.data(); + } + + std::vector BuildSizeGrid() const { + if (!args_.explicit_sizes.empty()) { + std::vector sizes = args_.explicit_sizes; + sizes.erase(std::remove(sizes.begin(), sizes.end(), 0), sizes.end()); + return sizes; + } + + const int low_exp = args_.min_exp; + const int high_exp = args_.max_exp; + if (args_.per_octave <= 0) { + std::vector sizes; + for (int exp = low_exp; exp <= high_exp; ++exp) { + sizes.push_back(std::size_t(1) << exp); + } + return sizes; + } + + std::set candidates; + const std::size_t min_size = std::size_t(1) << low_exp; + const std::size_t max_size = std::size_t(1) << high_exp; + + for (int exp = low_exp; exp <= high_exp; ++exp) { + for (int step = 0; step <= args_.per_octave; ++step) { + const long double x = + exp + static_cast(step) / args_.per_octave; + std::size_t size = + static_cast(std::llround(std::pow(2.0L, x))); + size = std::clamp(size, min_size, max_size); + if (size != 0) { + candidates.insert(size); + } + } + } + + std::vector sizes(candidates.begin(), candidates.end()); + if (sizes.front() != min_size) { + sizes.insert(sizes.begin(), min_size); + } + if (sizes.back() != max_size) { + sizes.push_back(max_size); + } + return sizes; + } + + std::shared_ptr BuildDataset(std::size_t N) { + std::mt19937_64 rng(args_.seed ^ + static_cast(N) * 0x9E3779B185EBCA87ull); + + auto data = std::make_shared(); + data->bits = RandomDyckBits(N, args_.p1, rng); + data->N = data->bits.size(); + N = data->N; + + data->words = PackWordsLsbFirst(data->bits); + const std::size_t block_bits = + args_.block_bits_specified ? args_.block_bits + : RmMBenchmarkTraits::DefaultBlockBits; + data->tree = + Tree(std::span(data->words), N, block_bits); + FillPools(*data, rng); + return data; + } + + void FillPools(Dataset& data, std::mt19937_64& rng) { + const std::size_t N = data.N; + const bool need_rank10 = ActiveOp("rank10"); + const bool need_select1 = ActiveOp("select1"); + const bool need_select0 = ActiveOp("select0"); + const bool need_select10 = ActiveOp("select10"); + const bool need_fwdsearch = ActiveOp("fwdsearch"); + const bool need_bwdsearch = ActiveOp("bwdsearch"); + const bool need_range_queries = + ActiveOp("range_min_query_pos") || ActiveOp("range_min_query_val") || + ActiveOp("range_max_query_pos") || ActiveOp("range_max_query_val") || + ActiveOp("mincount") || ActiveOp("minselect"); + const bool need_minselect = ActiveOp("minselect"); + const bool need_open_positions = ActiveOp("close") || ActiveOp("enclose"); + const bool need_close_positions = ActiveOp("open"); + + if (need_select1 || need_select0 || ActiveOp("rank0")) { + data.cnt1 = std::count(data.bits.begin(), data.bits.end(), '1'); + data.cnt0 = N - data.cnt1; + } + if (need_select10) { + data.cnt10 = Count10(data.bits); + } + + std::vector open_positions_zero_based; + std::vector open_positions_one_based; + std::vector close_positions_zero_based; + std::vector close_positions_one_based; + if (need_open_positions || need_close_positions) { + open_positions_zero_based.reserve(N >> 1); + open_positions_one_based.reserve(N >> 1); + close_positions_zero_based.reserve(N >> 1); + close_positions_one_based.reserve(N >> 1); + for (std::size_t i = 0; i < N; ++i) { + if (data.bits[i] == '1') { + open_positions_zero_based.push_back(i); + open_positions_one_based.push_back(i + 1); + } else { + close_positions_zero_based.push_back(i); + close_positions_one_based.push_back(i + 1); + } + } + } + + const std::size_t sample_count = + std::max(1, std::min(args_.Q, 32768)); + + std::uniform_int_distribution nonedge_index_distribution( + N > 1 ? 1 : 0, N > 1 ? (N - 1) : 0); + std::uniform_int_distribution ind_dist(0, N ? (N - 1) : 0); + std::uniform_int_distribution ind_dist_1N(1, N ? N : 1); + std::uniform_int_distribution rank10_end_position_distribution( + N >= 2 ? 2 : 0, N >= 2 ? N : 0); + std::uniform_int_distribution delta_distribution(-8, +8); + + if (ActiveOp("rank1") || ActiveOp("rank0") || ActiveOp("excess")) { + data.pool.inds_any.resize(sample_count); + } + if (need_rank10) { + data.pool.rank10_end_positions.resize(sample_count); + } + if (need_fwdsearch) { + data.pool.inds.resize(sample_count); + } + if (need_bwdsearch) { + data.pool.inds_1N.resize(sample_count); + } + if (need_fwdsearch || need_bwdsearch) { + data.pool.deltas.resize(sample_count); + } + if (need_range_queries) { + data.pool.segs.resize(sample_count); + } + + auto rand_ij = [&]() -> std::pair { + if (N == 0) { + return {0, 0}; + } + std::size_t left = ind_dist(rng); + std::size_t right = ind_dist(rng); + if (left > right) { + std::swap(left, right); + } + return {left, right}; + }; + + for (std::size_t i = 0; i < sample_count; ++i) { + if (!data.pool.inds_any.empty()) { + data.pool.inds_any[i] = nonedge_index_distribution(rng); + } + if (!data.pool.rank10_end_positions.empty()) { + data.pool.rank10_end_positions[i] = + rank10_end_position_distribution(rng); + } + if (!data.pool.inds.empty()) { + data.pool.inds[i] = (N ? ind_dist(rng) : 0); + } + if (!data.pool.inds_1N.empty()) { + data.pool.inds_1N[i] = (N ? ind_dist_1N(rng) : 0); + } + if (!data.pool.deltas.empty()) { + data.pool.deltas[i] = delta_distribution(rng); + } + if (!data.pool.segs.empty()) { + data.pool.segs[i] = rand_ij(); + } + } + + auto fill_from_candidates = [&](const std::vector& candidates, + std::vector& destination, + std::size_t fallback_value) { + destination.resize(sample_count); + if (candidates.empty()) { + std::fill(destination.begin(), destination.end(), fallback_value); + return; + } + std::uniform_int_distribution index_distribution( + 0, candidates.size() - 1); + for (std::size_t sample_index = 0; sample_index < sample_count; + ++sample_index) { + destination[sample_index] = candidates[index_distribution(rng)]; + } + }; + + const std::size_t one_based_fallback = (N > 0 ? 1 : 0); + if (ActiveOp("close")) { + fill_from_candidates(open_positions_zero_based, + data.pool.open_positions_zero_based, 0); + } + if (ActiveOp("enclose")) { + fill_from_candidates(open_positions_one_based, + data.pool.open_positions_one_based, + one_based_fallback); + } + if (ActiveOp("open")) { + fill_from_candidates(close_positions_zero_based, + data.pool.close_positions_zero_based, + one_based_fallback); + } + + auto fill_ks = [&](std::size_t total, std::vector& out) { + out.resize(sample_count); + if (total == 0) { + std::fill(out.begin(), out.end(), 1); + return; + } + std::uniform_int_distribution dist(1, total); + for (std::size_t i = 0; i < sample_count; ++i) { + out[i] = dist(rng); + } + }; + + if (need_select1) { + fill_ks(data.cnt1, data.pool.ks1); + } + if (need_select0) { + fill_ks(data.cnt0, data.pool.ks0); + } + if (need_select10) { + fill_ks(data.cnt10, data.pool.ks10); + } + + if (need_minselect) { + data.pool.minselect_q.resize(sample_count); + for (std::size_t i = 0; i < sample_count; ++i) { + auto [left, right] = data.pool.segs[i]; + std::size_t count = data.tree.mincount(left, right); + if (count == 0) { + count = 1; + } + std::uniform_int_distribution query_distribution(1, count); + data.pool.minselect_q[i] = query_distribution(rng); + } + } + } + + void RegisterAll() { + const auto sizes = BuildSizeGrid(); + bool any_op = false; + for (const auto& op : AllRmMOps()) { + any_op = any_op || ActiveOp(op); + } + if (!any_op) { + std::cerr << "No RmM benchmark operations match the requested filter.\n"; + return; + } + + for (std::size_t size : sizes) { + auto data = BuildDataset(size); + keepalive_.push_back(data); + + RegisterOp("rank1", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.rank1( + pool.inds_any[sample_index % pool.inds_any.size()]); + }); + RegisterOp("rank0", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.rank0( + pool.inds_any[sample_index % pool.inds_any.size()]); + }); + RegisterOp("select1", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.select1( + pool.ks1[sample_index % pool.ks1.size()]); + }); + RegisterOp("select0", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.select0( + pool.ks0[sample_index % pool.ks0.size()]); + }); + RegisterOp( + "rank10", data, [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.rank10( + pool.rank10_end_positions[sample_index % + pool.rank10_end_positions.size()]); + }); + RegisterOp("select10", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.select10( + pool.ks10[sample_index % pool.ks10.size()]); + }); + RegisterOp("excess", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.excess( + pool.inds_any[sample_index % pool.inds_any.size()]); + }); + RegisterOp("fwdsearch", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.fwdsearch( + pool.inds[sample_index % pool.inds.size()], + pool.deltas[sample_index % pool.deltas.size()]); + }); + RegisterOp("bwdsearch", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + return dataset.tree.bwdsearch( + pool.inds_1N[sample_index % pool.inds_1N.size()], + pool.deltas[sample_index % pool.deltas.size()]); + }); + RegisterOp("range_min_query_pos", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const auto [left, right] = + pool.segs[sample_index % pool.segs.size()]; + return dataset.tree.range_min_query_pos(left, right); + }); + RegisterOp("range_min_query_val", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const auto [left, right] = + pool.segs[sample_index % pool.segs.size()]; + return dataset.tree.range_min_query_val(left, right); + }); + RegisterOp("range_max_query_pos", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const auto [left, right] = + pool.segs[sample_index % pool.segs.size()]; + return dataset.tree.range_max_query_pos(left, right); + }); + RegisterOp("range_max_query_val", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const auto [left, right] = + pool.segs[sample_index % pool.segs.size()]; + return dataset.tree.range_max_query_val(left, right); + }); + RegisterOp("mincount", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const auto [left, right] = + pool.segs[sample_index % pool.segs.size()]; + return dataset.tree.mincount(left, right); + }); + RegisterOp("minselect", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const std::size_t index = sample_index % pool.segs.size(); + const auto [left, right] = pool.segs[index]; + return dataset.tree.minselect(left, right, + pool.minselect_q[index]); + }); + RegisterOp( + "close", data, [](const Dataset& dataset, std::size_t sample_index) { + if (dataset.N == 0) { + return std::size_t{0}; + } + const auto& pool = dataset.pool; + const std::size_t open_position = + pool.open_positions_zero_based + [sample_index % pool.open_positions_zero_based.size()]; + return dataset.tree.close(open_position); + }); + RegisterOp( + "open", data, [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const std::size_t close_position = + pool.close_positions_zero_based + [sample_index % pool.close_positions_zero_based.size()]; + return dataset.tree.open(close_position); + }); + RegisterOp( + "enclose", data, + [](const Dataset& dataset, std::size_t sample_index) { + const auto& pool = dataset.pool; + const std::size_t open_position = + pool.open_positions_zero_based + [sample_index % pool.open_positions_zero_based.size()]; + return dataset.tree.enclose(open_position); + }); + } + } + + template + void RegisterOp(const std::string& op, + std::shared_ptr data, + Fn&& body) { + if (!ActiveOp(op)) { + return; + } + + auto idx_ptr = std::make_shared(0); + + auto* benchmark = benchmark::RegisterBenchmark( + op.c_str(), [this, data, idx_ptr, + body = std::forward(body)](benchmark::State& state) { + const Dataset& dataset = *data; + for (auto _ : state) { + const std::size_t sample_index = (*idx_ptr)++; + auto result = body(dataset, sample_index); + benchmark::DoNotOptimize(result); + } + state.counters["N"] = static_cast(dataset.N); + state.counters["seed"] = static_cast(args_.seed); + const std::size_t block_bits = + args_.block_bits_specified + ? args_.block_bits + : RmMBenchmarkTraits::DefaultBlockBits; + state.counters["block_bits"] = static_cast(block_bits); + }); + + benchmark->Unit(benchmark::kNanosecond); + } + + RmMBenchmarkArgs args_; + std::vector> keepalive_; + std::vector extra_args_; + std::vector argv_vec_; +}; + +} // namespace pixie_bench diff --git a/src/docs/Doxyfile.in b/src/docs/Doxyfile.in index bffa501..297e6c7 100644 --- a/src/docs/Doxyfile.in +++ b/src/docs/Doxyfile.in @@ -1005,8 +1005,7 @@ WARN_LOGFILE = # spaces. See also FILE_PATTERNS and EXTENSION_MAPPING # Note: If this tag is empty the current directory is searched. -INPUT = @CMAKE_CURRENT_SOURCE_DIR@/include \ - @CMAKE_CURRENT_SOURCE_DIR@/README.md +INPUT = @PIXIE_DOC_INPUT@ # This tag can be used to specify the character encoding of the source files # that Doxygen parses. Internally Doxygen uses the UTF-8 encoding. Doxygen uses diff --git a/src/tests/test_rmm.cpp b/src/tests/test_rmm.cpp index 569a6f7..c8b4bfa 100644 --- a/src/tests/test_rmm.cpp +++ b/src/tests/test_rmm.cpp @@ -4,13 +4,12 @@ #include #include -#include #include -#include -#include #include #include +#include #include +#include #include #include @@ -101,37 +100,70 @@ static std::string random_dyck_bits(std::mt19937_64& rng, size_t m) { return s; } -struct Limits { - size_t CASES = 20; - size_t OPS_PER_CASE = 600; - size_t MAX_N = 20000; -}; +static size_t sdsl_style_close_ref(const std::string& bits, size_t position) { + if (position >= bits.size()) { + return NaiveRmM::npos; + } + if (bits[position] == '0') { + return position; + } + int balance = 1; + for (size_t i = position + 1; i < bits.size(); ++i) { + balance += bits[i] == '1' ? 1 : -1; + if (balance == 0) { + return i; + } + } + return NaiveRmM::npos; +} -static Limits load_limits_from_env() { - Limits L; - if (const char* s = std::getenv("RMM_CASES")) { - L.CASES = std::max(1, std::strtoull(s, nullptr, 10)); +static size_t sdsl_style_open_ref(const std::string& bits, size_t position) { + if (position >= bits.size()) { + return NaiveRmM::npos; } - if (const char* s = std::getenv("RMM_OPS")) { - L.OPS_PER_CASE = std::max(1, std::strtoull(s, nullptr, 10)); + if (bits[position] == '1') { + return position; } - if (const char* s = std::getenv("RMM_MAX_N")) { - L.MAX_N = std::max(1, std::strtoull(s, nullptr, 10)); + int balance = 0; + for (size_t i = position + 1; i > 0;) { + --i; + balance += bits[i] == '0' ? 1 : -1; + if (balance == 0 && bits[i] == '1') { + return i; + } } - return L; + return NaiveRmM::npos; } -static uint64_t choose_seed() { - if (const char* s = std::getenv("RMM_SEED")) { - return std::strtoull(s, nullptr, 10); +static size_t sdsl_style_enclose_ref(const std::string& bits, size_t position) { + if (position >= bits.size()) { + return NaiveRmM::npos; + } + if (bits[position] == '0') { + return sdsl_style_open_ref(bits, position); } - std::random_device rd; - return ((uint64_t)rd() << 32) ^ (uint64_t)rd() ^ - (uint64_t)std::chrono::high_resolution_clock::now() - .time_since_epoch() - .count(); + std::vector stack; + stack.reserve(position + 1); + for (size_t i = 0; i <= position; ++i) { + if (bits[i] == '1') { + if (i == position) { + return stack.empty() ? NaiveRmM::npos : stack.back(); + } + stack.push_back(i); + } else if (!stack.empty()) { + stack.pop_back(); + } + } + return NaiveRmM::npos; } +static constexpr uint64_t kSeed = 42; +static constexpr size_t kRandomCases = 20; +static constexpr size_t kOpsPerCase = 600; +static constexpr size_t kMaxBits = 20000; +static constexpr size_t kLongOpsPerCase = 2000; +static constexpr size_t kLongMaxBits = 65536; + static void run_case_and_compare(const std::string& bits, std::mt19937_64& rng, size_t ops_per_case, @@ -139,12 +171,12 @@ static void run_case_and_compare(const std::string& bits, const bool use_words = std::uniform_int_distribution(0, 1)(rng); pixie::RmMTree rm; NaiveRmM nv; + auto words = pack_words_lsb_first(bits); if (use_words) { - auto words = pack_words_lsb_first(bits); - rm = pixie::RmMTree(words, bits.size()); + rm = pixie::RmMTree(std::span(words), bits.size()); nv = NaiveRmM(words, bits.size()); } else { - rm = pixie::RmMTree(bits); + rm = pixie::RmMTree(std::span(words), bits.size()); nv = NaiveRmM(bits); } @@ -333,36 +365,27 @@ static void run_case_and_compare(const std::string& bits, class RmMRandomTest : public ::testing::Test { protected: - uint64_t seed{}; - std::mt19937_64 rng; - Limits L; - void SetUp() override { - L = load_limits_from_env(); - seed = choose_seed(); - rng.seed(seed); - std::cerr << "[ RmMRandomTest ] seed=" << seed << " CASES=" << L.CASES - << " OPS=" << L.OPS_PER_CASE << " MAX_N=" << L.MAX_N << "\n"; - } + std::mt19937_64 rng{kSeed}; }; TEST_F(RmMRandomTest, RandomMixedBits) { - std::uniform_int_distribution len_u(1, (int)L.MAX_N); - for (size_t t = 0; t < L.CASES; ++t) { + std::uniform_int_distribution len_u(1, (int)kMaxBits); + for (size_t t = 0; t < kRandomCases; ++t) { const size_t n = (size_t)len_u(rng); const std::string bits = random_bits(rng, n); - run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); + run_case_and_compare(bits, rng, kOpsPerCase, kSeed); } } TEST_F(RmMRandomTest, RandomDyckBits) { - std::uniform_int_distribution len_even(0, (int)(L.MAX_N / 2)); - for (size_t t = 0; t < L.CASES; ++t) { + std::uniform_int_distribution len_even(0, (int)(kMaxBits / 2)); + for (size_t t = 0; t < kRandomCases; ++t) { size_t m = 1 + (size_t)len_even(rng); std::string bits = random_dyck_bits(rng, m); if (bits.empty()) { bits = "10"; } - run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); + run_case_and_compare(bits, rng, kOpsPerCase, kSeed); } } @@ -378,25 +401,24 @@ TEST_F(RmMRandomTest, ShortInputs) { { SCOPED_TRACE(::testing::Message() << "short-inputs:string bits=" << bits); - run_case_and_compare(bits, rng, /*ops_per_case=*/200, - /*seed=*/0xC0FFEEull); + run_case_and_compare(bits, rng, /*ops_per_case=*/200, kSeed); } { SCOPED_TRACE(::testing::Message() << "short-inputs:words bits=" << bits); - run_case_and_compare(bits, rng, /*ops_per_case=*/200, - /*seed=*/0xC0FFEEull); + run_case_and_compare(bits, rng, /*ops_per_case=*/200, kSeed); } } } } -static void run_case_string_vs_words(const std::string& bits, - std::mt19937_64& rng) { +static void run_case_repeated_span_construction(const std::string& bits, + std::mt19937_64& rng) { const size_t N = bits.size(); - pixie::RmMTree rm_s(bits); + auto span_words = pack_words_lsb_first(bits); + pixie::RmMTree rm_s(std::span(span_words), N); auto words = pack_words_lsb_first(bits); - pixie::RmMTree rm_w(words, N); + pixie::RmMTree rm_w(std::span(words), N); std::uniform_int_distribution pos_i(0, N); std::uniform_int_distribution pos_i_nz(0, N ? N - 1 : 0); @@ -488,17 +510,18 @@ static void run_case_string_vs_words(const std::string& bits, } } -TEST_F(RmMRandomTest, WordsVsString) { - std::uniform_int_distribution len_u(0, (int)L.MAX_N); - for (size_t t = 0; t < L.CASES; ++t) { +TEST_F(RmMRandomTest, RepeatedSpanConstruction) { + std::uniform_int_distribution len_u(0, (int)kMaxBits); + for (size_t t = 0; t < kRandomCases; ++t) { const size_t n = (size_t)len_u(rng); const std::string bits = random_bits(rng, n); - run_case_string_vs_words(bits, rng); + run_case_repeated_span_construction(bits, rng); } } TEST(RmMEdgeCases, EmptyInput) { - pixie::RmMTree rm(std::string{}); + std::vector words; + pixie::RmMTree rm(std::span(words), 0); NaiveRmM nv(std::string{}); EXPECT_EQ(rm.rank1(0), nv.rank1(0)); EXPECT_EQ(rm.rank0(0), nv.rank0(0)); @@ -577,7 +600,9 @@ TEST(RmMEdgeCases, MultiwordPattern10AcrossWordBoundaries) { bits[boundary + 1] = '0'; } - pixie::RmMTree rm(bits, /*leaf_block_bits=*/256); + auto words = pack_words_lsb_first(bits); + pixie::RmMTree rm(std::span(words), bits.size(), + /*leaf_block_bits=*/256); NaiveRmM nv(bits); expect_rank_select_equal(rm, nv, n); @@ -591,7 +616,10 @@ TEST(RmMEdgeCases, PartialLastLeafSelects) { for (size_t i = 576; i < n; ++i) { mostly_zero[i] = '1'; } - pixie::RmMTree rm_select1(mostly_zero, /*leaf_block_bits=*/256); + auto mostly_zero_words = pack_words_lsb_first(mostly_zero); + pixie::RmMTree rm_select1(std::span(mostly_zero_words), + mostly_zero.size(), + /*leaf_block_bits=*/256); NaiveRmM nv_select1(mostly_zero); expect_rank_select_equal(rm_select1, nv_select1, n); @@ -599,7 +627,10 @@ TEST(RmMEdgeCases, PartialLastLeafSelects) { for (size_t i = 576; i < n; ++i) { mostly_one[i] = '0'; } - pixie::RmMTree rm_select0(mostly_one, /*leaf_block_bits=*/256); + auto mostly_one_words = pack_words_lsb_first(mostly_one); + pixie::RmMTree rm_select0(std::span(mostly_one_words), + mostly_one.size(), + /*leaf_block_bits=*/256); NaiveRmM nv_select0(mostly_one); expect_rank_select_equal(rm_select0, nv_select0, n); } @@ -615,7 +646,9 @@ TEST(RmMEdgeCases, Select10OnIncompleteInternalNode) { bits[i + 1] = '0'; } - pixie::RmMTree rm(bits, leaf_block_bits); + auto words = pack_words_lsb_first(bits); + pixie::RmMTree rm(std::span(words), bits.size(), + leaf_block_bits); NaiveRmM nv(bits); const size_t pairs10 = nv.rank10(n); @@ -636,7 +669,9 @@ TEST(RmMEdgeCases, InvalidArgumentsGuards) { bits[i] = '0'; } - pixie::RmMTree rm(bits, /*leaf_block_bits=*/256); + auto words = pack_words_lsb_first(bits); + pixie::RmMTree rm(std::span(words), bits.size(), + /*leaf_block_bits=*/256); EXPECT_EQ(rm.select1(0), pixie::RmMTree::npos); EXPECT_EQ(rm.select0(0), pixie::RmMTree::npos); @@ -665,23 +700,52 @@ TEST(RmMEdgeCases, InvalidArgumentsGuards) { /** * bit_count is larger than the provided words buffer. - * Verifies that words beyond the provided buffer are treated as zeros after - * resize. + * Verifies that non-owning construction rejects spans that are too short. */ -TEST(RmMEdgeCases, WordsConstructorResizesInputStorage) { +TEST(RmMEdgeCases, SpanConstructorRejectsShortInputStorage) { std::vector words = {0xAAAAAAAAAAAAAAAAull}; const size_t bit_count = 300; - pixie::RmMTree rm(words, bit_count); - NaiveRmM nv(words, bit_count); + EXPECT_THROW( + (void)pixie::RmMTree(std::span(words), bit_count), + std::invalid_argument); +} + +#ifdef SDSL_SUPPORT +TEST(RmMSdslEdgeCases, IgnoresWordsBeyondBitCount) { + std::vector words = { + 0b101ull, std::numeric_limits::max()}; + pixie::SdslRmMTree rm(std::span(words), + /*bit_count=*/3, + /*unused=*/0); + + EXPECT_EQ(rm.size(), 3u); + EXPECT_EQ(rm.rank1(3), 2u); + EXPECT_EQ(rm.rank0(3), 1u); + EXPECT_EQ(rm.select1(1), 0u); + EXPECT_EQ(rm.select1(2), 2u); + EXPECT_EQ(rm.select1(3), pixie::SdslRmMTree::npos); +} + +TEST(RmMSdslEdgeCases, ParenthesesNavigationMatchesSdslStyleRefs) { + const std::string bits = "11101001011000"; + auto words = pack_words_lsb_first(bits); + pixie::SdslRmMTree rm(std::span(words), bits.size(), + /*unused=*/0); - expect_rank_select_equal(rm, nv, bit_count); - expect_range_ops_equal(rm, nv, bit_count); + for (size_t position = 0; position < bits.size(); ++position) { + SCOPED_TRACE(::testing::Message() + << "position=" << position << " bits=" << bits); + EXPECT_EQ(rm.close(position), sdsl_style_close_ref(bits, position)); + EXPECT_EQ(rm.open(position), sdsl_style_open_ref(bits, position)); + EXPECT_EQ(rm.enclose(position), sdsl_style_enclose_ref(bits, position)); + } } +#endif /** * Same bitvector built through different configuration paths (auto vs explicit - * leaf size, different overhead caps, and words-based constructor). Query + * leaf size and different overhead caps). Query * results must be identical. */ TEST(RmMEdgeCases, ExplicitBuildParametersAndOverheadCap) { @@ -690,57 +754,103 @@ TEST(RmMEdgeCases, ExplicitBuildParametersAndOverheadCap) { const std::string bits = random_bits(rng, n); NaiveRmM nv(bits); - pixie::RmMTree rm_auto(bits, /*leaf_block_bits=*/0, /*max_overhead=*/1.f); - pixie::RmMTree rm_explicit(bits, /*leaf_block_bits=*/512, - /*max_overhead=*/2.f); auto words = pack_words_lsb_first(bits); - pixie::RmMTree rm_words(words, n, /*leaf_block_bits=*/256, - /*max_overhead=*/1.f); + const std::span word_span(words); + pixie::RmMTree rm_auto(word_span, n, 0, 1.f); + pixie::RmMTree rm_explicit(word_span, n, 512, 2.f); + pixie::RmMTree rm_span(word_span, n, 256, 1.f); expect_rank_select_equal(rm_auto, nv, n); expect_range_ops_equal(rm_auto, nv, n); expect_rank_select_equal(rm_explicit, nv, n); expect_range_ops_equal(rm_explicit, nv, n); - expect_rank_select_equal(rm_words, nv, n); - expect_range_ops_equal(rm_words, nv, n); + expect_rank_select_equal(rm_span, nv, n); + expect_range_ops_equal(rm_span, nv, n); } -TEST(RmMTreeStress, LongRandom) { - Limits L; - L.OPS_PER_CASE = 2000; - L.MAX_N = 65536; - - if (const char* s = std::getenv("RMM_CASES")) { - L.CASES = std::max(1, std::strtoull(s, nullptr, 10)); +TEST(RmMEdgeCases, BoundaryHeavyQueries) { + constexpr size_t leaf_block_bits = 128; + const size_t n = 777; + std::string bits(n, '0'); + for (size_t i = 0; i < n; ++i) { + bits[i] = (((i * 37 + i / 7) % 11) < 5) ? '1' : '0'; } - if (const char* s = std::getenv("RMM_OPS")) { - L.OPS_PER_CASE = std::max(1, std::strtoull(s, nullptr, 10)); + for (size_t position : {63u, 64u, 65u, 127u, 128u, 129u, 255u, 256u}) { + bits[position] = (position & 1u) ? '1' : '0'; } - if (const char* s = std::getenv("RMM_MAX_N")) { - L.MAX_N = std::max(1, std::strtoull(s, nullptr, 10)); + + auto words = pack_words_lsb_first(bits); + pixie::RmMTree rm(std::span(words), n, leaf_block_bits); + NaiveRmM nv(bits); + + for (size_t end_position : std::array{0, 1, 63, 64, 65, 127, 128, + 129, 255, 256, 512, n}) { + EXPECT_EQ(rm.rank1(end_position), nv.rank1(end_position)); + EXPECT_EQ(rm.rank0(end_position), nv.rank0(end_position)); + EXPECT_EQ(rm.rank10(end_position), nv.rank10(end_position)); } - const uint64_t seed = choose_seed(); - std::mt19937_64 rng(seed); - size_t LOG_EVERY = 5; - if (const char* s = std::getenv("RMM_LOG_EVERY")) { - size_t v = std::strtoull(s, nullptr, 10); - if (v) { - LOG_EVERY = v; + const size_t ones = nv.rank1(n); + const size_t zeros = n - ones; + const size_t pairs10 = nv.rank10(n); + for (size_t rank : std::array{1, 2, 7, 31, 64, ones, ones + 1}) { + EXPECT_EQ(rm.select1(rank), nv.select1(rank)) << "select1 rank=" << rank; + } + for (size_t rank : std::array{1, 2, 7, 31, 64, zeros, zeros + 1}) { + EXPECT_EQ(rm.select0(rank), nv.select0(rank)) << "select0 rank=" << rank; + } + for (size_t rank : + std::array{1, 2, 7, 31, 64, pairs10, pairs10 + 1}) { + EXPECT_EQ(rm.select10(rank), nv.select10(rank)) << "select10 rank=" << rank; + } + + for (auto [left, right] : + std::array, 10>{{{0, 0}, + {57, 75}, + {63, 130}, + {120, 260}, + {128, 511}, + {129, 511}, + {250, 520}, + {512, n - 1}, + {0, n - 1}, + {n - 9, n - 1}}}) { + EXPECT_EQ(rm.range_min_query_pos(left, right), + nv.range_min_query_pos(left, right)); + EXPECT_EQ(rm.range_min_query_val(left, right), + nv.range_min_query_val(left, right)); + EXPECT_EQ(rm.range_max_query_pos(left, right), + nv.range_max_query_pos(left, right)); + EXPECT_EQ(rm.range_max_query_val(left, right), + nv.range_max_query_val(left, right)); + + const size_t count = nv.mincount(left, right); + EXPECT_EQ(rm.mincount(left, right), count); + EXPECT_EQ(rm.minselect(left, right, 1), nv.minselect(left, right, 1)); + EXPECT_EQ(rm.minselect(left, right, count), + nv.minselect(left, right, count)); + EXPECT_EQ(rm.minselect(left, right, count + 1), + nv.minselect(left, right, count + 1)); + } + + for (size_t position : std::array{0, 1, 63, 64, 65, 127, 128, 129, + 255, 256, 512, n - 1}) { + for (int delta : {-40, -9, -8, -2, -1, 0, 1, 2, 8, 9, 40}) { + EXPECT_EQ(rm.fwdsearch(position, delta), nv.fwdsearch(position, delta)) + << "fwdsearch position=" << position << " delta=" << delta; + EXPECT_EQ(rm.bwdsearch(position + 1, delta), + nv.bwdsearch(position + 1, delta)) + << "bwdsearch position=" << (position + 1) << " delta=" << delta; } } +} - std::cerr << "[ LongRandom ] seed=" << seed << " CASES=" << L.CASES - << " OPS=" << L.OPS_PER_CASE << " MAX_N=" << L.MAX_N - << " LOG_EVERY=" << LOG_EVERY << "\n"; - +TEST_F(RmMRandomTest, LongRandom) { std::uniform_int_distribution coin(0, 1); - std::uniform_int_distribution len_u(1, (int)L.MAX_N); - std::uniform_int_distribution len_even(0, (int)(L.MAX_N / 2)); - - size_t total_ops = 0; + std::uniform_int_distribution len_u(1, (int)kLongMaxBits); + std::uniform_int_distribution len_even(0, (int)(kLongMaxBits / 2)); - for (size_t iter = 1; iter <= L.CASES; ++iter) { + for (size_t iter = 0; iter < kRandomCases; ++iter) { std::string bits; if (coin(rng) == 0) { bits = random_bits(rng, (size_t)len_u(rng)); @@ -752,19 +862,13 @@ TEST(RmMTreeStress, LongRandom) { } } - run_case_and_compare(bits, rng, L.OPS_PER_CASE, seed); - total_ops += L.OPS_PER_CASE; - - if (iter % LOG_EVERY == 0) { - std::cerr << "[ LongRandom ] iter=" << iter << " total_ops=" << total_ops - << " last_N=" << bits.size() << " ok\n"; - } + run_case_and_compare(bits, rng, kLongOpsPerCase, kSeed); } } TEST(RmMTest, RankBasic) { std::vector bits = {0b10110}; - pixie::RmMTree rm(bits, 5); + pixie::RmMTree rm(std::span(bits), 5); EXPECT_EQ(rm.rank1(0), 0); // No bits EXPECT_EQ(rm.rank1(1), 0); // 0 @@ -776,7 +880,7 @@ TEST(RmMTest, RankBasic) { TEST(RmMTest, RankWithZeros) { std::vector bits = {0}; - pixie::RmMTree rm(bits, 5); + pixie::RmMTree rm(std::span(bits), 5); for (size_t i = 0; i <= 5; i++) { EXPECT_EQ(rm.rank1(i), 0); @@ -785,7 +889,7 @@ TEST(RmMTest, RankWithZeros) { TEST(RmMTest, SelectBasic) { std::vector bits = {0b1100010110010110}; - pixie::RmMTree rm(bits, 16); + pixie::RmMTree rm(std::span(bits), 16); EXPECT_EQ(rm.select1(1), 1); EXPECT_EQ(rm.select1(2), 2); @@ -805,7 +909,7 @@ TEST(RmMTest, MainRankTest) { } size_t rm_size = 65536 * 32 * 64; - pixie::RmMTree rm(bits, rm_size); + pixie::RmMTree rm(std::span(bits), rm_size); size_t rank = 0; for (size_t i = 0; i < rm_size; ++i) { ASSERT_EQ(rank, rm.rank1(i)); @@ -821,7 +925,7 @@ TEST(RmMTest, MainSelectTest) { } size_t rm_size = 65536 * 32 * 64; - pixie::RmMTree rm(bits, rm_size); + pixie::RmMTree rm(std::span(bits), rm_size); size_t rank = 0; for (size_t i = 0; i < rm_size; ++i) {