diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml index 247d974..2d9e248 100644 --- a/.github/workflows/ci.yaml +++ b/.github/workflows/ci.yaml @@ -1,83 +1,186 @@ name: ci -on: [ push, pull_request ] +on: + push: + branches: [ master ] + pull_request: env: CARGO_TERM_COLOR: always + RUSTUP_TOOLCHAIN: stable jobs: + msrv: + name: MSRV 1.59 + runs-on: ubuntu-latest + env: + RUSTUP_TOOLCHAIN: 1.59.0 + steps: + - name: Checkout sources + uses: actions/checkout@v7 + + - name: Install MSRV toolchain + uses: dtolnay/rust-toolchain@1.59.0 + with: + targets: aarch64-unknown-linux-gnu + + - name: Cache dependencies + uses: Swatinem/rust-cache@v2 + + - name: Test on the MSRV + run: cargo test --all-features --locked + + - name: Check AArch64 NEON on the MSRV + run: cargo check --lib --all-features --locked --target aarch64-unknown-linux-gnu + check: name: Check runs-on: ubuntu-latest steps: - name: Checkout sources - uses: actions/checkout@v2 + uses: actions/checkout@v7 - - name: Install pinned toolchain - uses: actions-rs/toolchain@v1 - with: - profile: minimal - toolchain: 1.95.0 - override: true + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable - name: Cache dependencies - uses: Swatinem/rust-cache@v1 + uses: Swatinem/rust-cache@v2 - name: Run cargo check - uses: actions-rs/cargo@v1 - with: - command: check - args: --all-features + run: cargo check --all-targets --all-features --locked lints: name: Lints runs-on: ubuntu-latest steps: - name: Checkout sources - uses: actions/checkout@v2 + uses: actions/checkout@v7 - - name: Install pinned toolchain - uses: actions-rs/toolchain@v1 + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable with: - profile: minimal - toolchain: 1.95.0 - override: true - components: rustfmt, clippy + components: rustfmt,clippy - name: Cache dependencies - uses: Swatinem/rust-cache@v1 + uses: Swatinem/rust-cache@v2 - name: Run cargo fmt - uses: actions-rs/cargo@v1 - with: - command: fmt - args: --all -- --check + run: cargo fmt --all -- --check - name: Run cargo clippy - uses: actions-rs/cargo@v1 + run: cargo clippy --all-targets --all-features --locked -- -D warnings + + native-avx2: + name: Native AVX2 + runs-on: ubuntu-latest + steps: + - name: Checkout sources + uses: actions/checkout@v7 + + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable + + - name: Cache dependencies + uses: Swatinem/rust-cache@v2 + + - name: Record and require AVX2 + run: | + lscpu + grep -qw avx2 /proc/cpuinfo + + - name: Exercise runtime dispatch and AVX2 masks + run: | + cargo test --all-features --locked simd::dispatch_tests::selected_backend_matches_the_host -- --exact --nocapture + cargo test --all-features --locked simd::avx2::tests::primitive_masks_match_scalar_when_avx2_is_available -- --exact --nocapture + + scalar-x86: + name: Scalar x86-64 + runs-on: ubuntu-latest + steps: + - name: Checkout sources + uses: actions/checkout@v7 + + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable + + - name: Cache dependencies + uses: Swatinem/rust-cache@v2 + + - name: Install QEMU + run: | + sudo apt-get update + sudo apt-get install --yes qemu-user + qemu-x86_64 --version + + - name: Exercise scalar runtime dispatch without AVX2 + env: + CARGO_TARGET_X86_64_UNKNOWN_LINUX_GNU_RUNNER: qemu-x86_64 -cpu Nehalem + run: | + cargo test --all-features --locked --target x86_64-unknown-linux-gnu simd::dispatch_tests::selected_backend_matches_the_host -- --exact --nocapture + cargo test --all-features --locked --target x86_64-unknown-linux-gnu + + scalar-targets: + name: Scalar fallback targets + runs-on: ubuntu-latest + steps: + - name: Checkout sources + uses: actions/checkout@v7 + + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable with: - command: clippy - args: --all-features -- -D warnings + targets: i686-unknown-linux-gnu,powerpc64le-unknown-linux-gnu,wasm32-unknown-unknown + + - name: Cache dependencies + uses: Swatinem/rust-cache@v2 + + - name: Check scalar-only architectures + run: | + cargo check --lib --all-features --locked --target i686-unknown-linux-gnu + cargo check --lib --all-features --locked --target powerpc64le-unknown-linux-gnu + cargo check --lib --all-features --locked --target wasm32-unknown-unknown test: - name: Test Suite + name: Test Suite (${{ matrix.name }}) runs-on: ${{ matrix.os }} strategy: matrix: - os: [ ubuntu-latest, macOS-latest, windows-latest ] + include: + - name: Linux, default features + os: ubuntu-latest + features: "" + - name: Linux, serde + os: ubuntu-latest + features: --features with_serde + - name: Linux, unstable sort + os: ubuntu-latest + features: --features sort_unstable + - name: Linux, all features + os: ubuntu-latest + features: --all-features + - name: macOS ARM64, all features + os: macos-15 + features: --all-features + neon: true + - name: Windows, all features + os: windows-latest + features: --all-features steps: - name: Checkout sources - uses: actions/checkout@v2 + uses: actions/checkout@v7 - - name: Install pinned toolchain - uses: actions-rs/toolchain@v1 - with: - profile: minimal - toolchain: 1.95.0 - override: true + - name: Install stable toolchain + uses: dtolnay/rust-toolchain@stable - name: Cache dependencies - uses: Swatinem/rust-cache@v1 + uses: Swatinem/rust-cache@v2 - name: Run tests - run: cargo test --verbose --all-features \ No newline at end of file + run: cargo test --verbose --locked ${{ matrix.features }} + + - name: Verify native NEON dispatch + if: matrix.neon + run: | + test "$(uname -m)" = arm64 + cargo test --all-features --locked simd::dispatch_tests::selected_backend_matches_the_host -- --exact --nocapture + cargo test --all-features --locked --test portable_index every_primitive_integer_type_and_block_tail_matches_forward_brute_force -- --exact diff --git a/Cargo.lock b/Cargo.lock index 802d636..cc0506e 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,12 @@ # It is not intended for manual editing. version = 3 +[[package]] +name = "anes" +version = "0.1.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "4b46cbb362ab8752921c97e041f5e366ee6297bd428a31275b9fcf1e380f7299" + [[package]] name = "atty" version = "0.2.14" @@ -15,9 +21,9 @@ dependencies = [ [[package]] name = "autocfg" -version = "1.0.1" +version = "1.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cdb031dd78e28731d87d56cc8ffef4a8f36ca26c38fe2de700543e627f8a464a" +checksum = "f2032f911046de80f0a198e0901378627c33f59ea0ac00e363d481118bd70a53" [[package]] name = "bincode" @@ -30,84 +36,88 @@ dependencies = [ [[package]] name = "bitflags" -version = "1.3.1" +version = "1.3.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2da1976d75adbe5fbc88130ecd119529cf1cc6a93ae1546d8696ee66f0d21af1" +checksum = "bef38d45163c2f1dde094a7dfd33ccf595c92905c8f8f4fdc18d06fb1037718a" [[package]] -name = "bstr" -version = "0.2.16" +name = "cast" +version = "0.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "90682c8d613ad3373e66de8c6411e0ae2ab2571e879d2efbf73558cc66f21279" -dependencies = [ - "lazy_static", - "memchr", - "regex-automata", - "serde", -] +checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" [[package]] -name = "bumpalo" -version = "3.7.0" +name = "cfg-if" +version = "1.0.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9c59e7af012c713f529e7a3ee57ce9b31ddd858d4b512923602f74608b009631" +checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] -name = "cast" -version = "0.2.7" +name = "ciborium" +version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4c24dab4283a142afa2fdca129b80ad2c6284e073930f964c3a1293c225ee39a" +checksum = "42e69ffd6f0917f5c029256a24d0161db17cea3997d185db0d35926308770f0e" dependencies = [ - "rustc_version", + "ciborium-io", + "ciborium-ll", + "serde", ] [[package]] -name = "cfg-if" -version = "1.0.0" +name = "ciborium-io" +version = "0.2.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" +checksum = "05afea1e0a06c9be33d539b876f1ce3692f4afea2cb41f740e7743225ed1c757" + +[[package]] +name = "ciborium-ll" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "57663b653d948a338bfb3eeba9bb2fd5fcfaecb9e199e87e1eda4d9e8b240fd9" +dependencies = [ + "ciborium-io", + "half", +] [[package]] name = "clap" -version = "2.33.3" +version = "3.2.25" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "37e58ac78573c40708d45522f0d80fa2f01cc4f9b4e2bf749807255454312002" +checksum = "4ea181bf566f71cb9a5d17a59e1871af638180a18fb0035c92ae62b705207123" dependencies = [ "bitflags", + "clap_lex", + "indexmap", "textwrap", - "unicode-width", ] [[package]] -name = "cpu-time" -version = "0.1.0" +name = "clap_lex" +version = "0.2.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ad9182963eedd274a445d6a43a50c4097537a238c8ad8980e400c3bfbc956426" +checksum = "2850f2f5a82cbf437dd5af4d49848fbdfc27c157c3d010345776f952765261c5" dependencies = [ - "libc", - "winapi", + "os_str_bytes", ] [[package]] name = "criterion" -version = "0.3.5" +version = "0.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1604dafd25fba2fe2d5895a9da139f8dc9b319a5fe5354ca137cbbce4e178d10" +checksum = "e7c76e09c1aae2bc52b3d2f29e13c6572553b30c4aa1b8a49fd70de6412654cb" dependencies = [ + "anes", "atty", "cast", + "ciborium", "clap", "criterion-plot", - "csv", "itertools", "lazy_static", "num-traits", "oorandom", - "plotters", - "rayon", "regex", "serde", - "serde_cbor", "serde_derive", "serde_json", "tinytemplate", @@ -116,267 +126,159 @@ dependencies = [ [[package]] name = "criterion-plot" -version = "0.4.4" +version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d00996de9f2f7559f7f4dc286073197f83e92256a59ed395f9aac01fe717da57" +checksum = "6b50826342786a51a89e2da3a28f1c32b06e387201bc2d19791f622c673706b1" dependencies = [ "cast", "itertools", ] [[package]] -name = "crossbeam-channel" -version = "0.5.1" +name = "crunchy" +version = "0.2.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "06ed27e177f16d65f0f0c22a213e17c696ace5dd64b14258b52f9417ccb52db4" -dependencies = [ - "cfg-if", - "crossbeam-utils", -] +checksum = "460fbee9c2c2f33933d720630a6a0bac33ba7053db5344fac858d4b8952d77d5" [[package]] -name = "crossbeam-deque" -version = "0.8.1" +name = "either" +version = "1.10.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6455c0ca19f0d2fbf751b908d5c55c1f5cbc65e03c4225427254b46890bdde1e" -dependencies = [ - "cfg-if", - "crossbeam-epoch", - "crossbeam-utils", -] +checksum = "11157ac094ffbdde99aa67b23417ebdd801842852b500e395a45a9c0aac03e4a" [[package]] -name = "crossbeam-epoch" -version = "0.9.5" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "4ec02e091aa634e2c3ada4a392989e7c3116673ef0ac5b72232439094d73b7fd" -dependencies = [ - "cfg-if", - "crossbeam-utils", - "lazy_static", - "memoffset", - "scopeguard", -] - -[[package]] -name = "crossbeam-utils" -version = "0.8.8" +name = "getrandom" +version = "0.2.15" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0bf124c720b7686e3c2663cf54062ab0f68a88af2fb6a030e87e30bf721fcb38" +checksum = "c4567c8db10ae91089c99af84c68c38da3ec2f087c3f82960bcdbf3656b6f4d7" dependencies = [ "cfg-if", - "lazy_static", -] - -[[package]] -name = "csv" -version = "1.1.6" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "22813a6dc45b335f9bade10bf7271dc477e81113e89eb251a0bc2a8a81c536e1" -dependencies = [ - "bstr", - "csv-core", - "itoa", - "ryu", - "serde", + "libc", + "wasi", ] [[package]] -name = "csv-core" -version = "0.1.10" +name = "half" +version = "2.2.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b2466559f260f48ad25fe6317b3c8dac77b5bdb5763ac7d9d6103530663bc90" +checksum = "02b4af3693f1b705df946e9fe5631932443781d0aabb423b62fcd4d73f6d2fd0" dependencies = [ - "memchr", + "crunchy", ] [[package]] -name = "either" -version = "1.6.1" +name = "hashbrown" +version = "0.12.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e78d4f1cc4ae33bbfc157ed5d5a5ef3bc29227303d595861deb238fcec4e9457" +checksum = "8a9ee70c43aaf417c914396645a0fa852624801b24ebb7ae78fe8272889ac888" [[package]] -name = "getrandom" -version = "0.1.16" +name = "hermit-abi" +version = "0.1.19" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8fc3cb4d91f53b50155bdcfd23f6a4c39ae1969c2ae85982b135750cccaf5fce" +checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" dependencies = [ - "cfg-if", "libc", - "wasi", ] [[package]] -name = "half" -version = "1.7.1" +name = "indexmap" +version = "1.9.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "62aca2aba2d62b4a7f5b33f3712cb1b0692779a56fb510499d5c0aa594daeaf3" - -[[package]] -name = "hermit-abi" -version = "0.1.19" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "62b467343b94ba476dcb2500d242dadbb39557df889310ac77c5d99100aaac33" +checksum = "bd070e393353796e801d209ad339e89596eb4c8d430d18ede6a1cced8fafbd99" dependencies = [ - "libc", + "autocfg", + "hashbrown", ] [[package]] name = "itertools" -version = "0.10.1" +version = "0.10.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "69ddb889f9d0d08a67338271fa9b62996bc788c7796a5c18cf057420aaed5eaf" +checksum = "b0fd2260e829bddf4cb6ea802289de2f86d6a7a690192fbe91b3f46e0f2c8473" dependencies = [ "either", ] [[package]] name = "itoa" -version = "0.4.7" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "dd25036021b0de88a0aff6b850051563c6516d0bf53f8638938edbb9de732736" - -[[package]] -name = "js-sys" -version = "0.3.52" +version = "1.0.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ce791b7ca6638aae45be056e068fc756d871eb3b3b10b8efa62d1c9cec616752" -dependencies = [ - "wasm-bindgen", -] +checksum = "49f1f14873335454500d59611f1cf4a4b0f786f9ac11f4312a78e4cf2566695b" [[package]] name = "lazy_static" -version = "1.4.0" +version = "1.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e2abad23fbc42b3700f2f279844dc832adb2b2eb069b2df918f455c4e18cc646" +checksum = "bbd2bcb4c963f2ddae06a2efc7e9f3591312473c50c6685e1f298068316e66fe" [[package]] name = "libc" -version = "0.2.99" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a7f823d141fe0a24df1e23b4af4e3c7ba9e5966ec514ea068c93024aa7deb765" - -[[package]] -name = "log" -version = "0.4.14" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "51b9bbe6c47d51fc3e1a9b945965946b4c44142ab8792c50835a980d362c2710" -dependencies = [ - "cfg-if", -] - -[[package]] -name = "memchr" -version = "2.4.0" +version = "0.2.163" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b16bd47d9e329435e309c58469fe0791c2d0d1ba96ec0954152a5ae2b04387dc" - -[[package]] -name = "memoffset" -version = "0.6.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "59accc507f1338036a0477ef61afdae33cde60840f4dfe481319ce3ad116ddf9" -dependencies = [ - "autocfg", -] +checksum = "1fdaeca4cf44ed4ac623e86ef41f056e848dbeab7ec043ecb7326ba300b36fd0" [[package]] name = "num-traits" -version = "0.2.14" +version = "0.2.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9a64b1ec5cda2586e284722486d802acf1f7dbdc623e2bfc57e65ca1cd099290" +checksum = "da0df0e5185db44f69b44f26786fe401b6c293d1907744beaa7fa62b2e5a517a" dependencies = [ "autocfg", ] -[[package]] -name = "num_cpus" -version = "1.13.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "05499f3756671c15885fee9034446956fff3f243d6077b91e5767df161f766b3" -dependencies = [ - "hermit-abi", - "libc", -] - [[package]] name = "oorandom" -version = "11.1.3" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0ab1bc2a289d34bd04a330323ac98a1b4bc82c9d9fcb1e66b63caa84da26b575" - -[[package]] -name = "plotters" -version = "0.3.1" +version = "11.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "32a3fd9ec30b9749ce28cd91f255d569591cdf937fe280c312143e3c4bad6f2a" -dependencies = [ - "num-traits", - "plotters-backend", - "plotters-svg", - "wasm-bindgen", - "web-sys", -] +checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" [[package]] -name = "plotters-backend" -version = "0.3.2" +name = "os_str_bytes" +version = "6.5.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d88417318da0eaf0fdcdb51a0ee6c3bed624333bff8f946733049380be67ac1c" - -[[package]] -name = "plotters-svg" -version = "0.3.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "521fa9638fa597e1dc53e9412a4f9cefb01187ee1f7413076f9e6749e2885ba9" -dependencies = [ - "plotters-backend", -] +checksum = "4d5d9eb14b174ee9aa2ef96dc2b94637a2d4b6e7cb873c7e171f0c20c6cf3eac" [[package]] name = "ppv-lite86" -version = "0.2.10" +version = "0.2.17" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ac74c624d6b2d21f425f752262f42188365d7b8ff1aff74c82e45136510a4857" +checksum = "5b40af805b3121feab8a3c29f04d8ad262fa8e0561883e7653e024ae4479e6de" [[package]] name = "proc-macro2" -version = "1.0.28" +version = "1.0.80" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5c7ed8b8c7b886ea3ed7dde405212185f423ab44682667c8c6dd14aa1d9f6612" +checksum = "a56dea16b0a29e94408b9aa5e2940a4eedbd128a1ba20e8f7ae60fd3d465af0e" dependencies = [ - "unicode-xid", + "unicode-ident", ] [[package]] name = "quote" -version = "1.0.9" +version = "1.0.35" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c3d0b9745dc2debf507c8422de05d7226cc1f0644216dfdfead988f9b1ab32a7" +checksum = "291ec9ab5efd934aaf503a6466c5d5251535d108ee747472c3977cc5acc868ef" dependencies = [ "proc-macro2", ] [[package]] name = "rand" -version = "0.7.3" +version = "0.8.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6a6b1679d49b24bbfe0c803429aa1874472f50d9b363131f0e89fc356b544d03" +checksum = "5ca0ecfa931c29007047d1bc58e623ab12e5590e8c7cc53200d5202b69266d8a" dependencies = [ - "getrandom", "libc", "rand_chacha", "rand_core", - "rand_hc", ] [[package]] name = "rand_chacha" -version = "0.2.2" +version = "0.3.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f4c8ed856279c9737206bf725bf36935d8666ead7aa69b52be55af369d193402" +checksum = "e6c10a63a0fa32252be49d21e7709d4d4baf8d231c2dbce1eaa8141b9b127d88" dependencies = [ "ppv-lite86", "rand_core", @@ -384,94 +286,44 @@ dependencies = [ [[package]] name = "rand_core" -version = "0.5.1" +version = "0.6.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "90bde5296fc891b0cef12a6d03ddccc162ce7b2aff54160af9338f8d40df6d19" +checksum = "ec0be4795e2f6a28069bec0b5ff3e2ac9bafc99e6a9a7dc3547996c5c816922c" dependencies = [ "getrandom", ] -[[package]] -name = "rand_hc" -version = "0.2.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "ca3129af7b92a17112d59ad498c6f81eaf463253766b90396d39ea7a39d6613c" -dependencies = [ - "rand_core", -] - -[[package]] -name = "rayon" -version = "1.5.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "c06aca804d41dbc8ba42dfd964f0d01334eceb64314b9ecf7c5fad5188a06d90" -dependencies = [ - "autocfg", - "crossbeam-deque", - "either", - "rayon-core", -] - -[[package]] -name = "rayon-core" -version = "1.9.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d78120e2c850279833f1dd3582f730c4ab53ed95aeaaaa862a2a5c71b1656d8e" -dependencies = [ - "crossbeam-channel", - "crossbeam-deque", - "crossbeam-utils", - "lazy_static", - "num_cpus", -] - [[package]] name = "regex" -version = "1.5.6" +version = "1.7.3" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d83f127d94bdbcda4c8cc2e50f6f84f4b611f69c902699ca385a39c3a75f9ff1" +checksum = "8b1f693b24f6ac912f4893ef08244d70b6067480d2f1a46e950c9691e6749d1d" dependencies = [ "regex-syntax", ] -[[package]] -name = "regex-automata" -version = "0.1.10" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6c230d73fb8d8c1b9c0b3135c5142a8acee3a0558fb8db5cf1cb65f8d7862132" - [[package]] name = "regex-syntax" -version = "0.6.26" +version = "0.6.29" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "49b3de9ec5dc0a3417da371aab17d729997c15010e7fd24ff707773a33bddb64" +checksum = "f162c6dd7b008981e4d40210aca20b4bd0f9b60ca9271061b07f78537722f2e1" [[package]] name = "rust-lapper" -version = "1.3.0" +version = "2.0.0-beta.1" dependencies = [ "bincode", - "cpu-time", "criterion", "num-traits", "rand", "serde", ] -[[package]] -name = "rustc_version" -version = "0.4.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "bfa0f585226d2e68097d4f95d113b15b83a82e819ab25717ec0590d9584ef366" -dependencies = [ - "semver", -] - [[package]] name = "ryu" -version = "1.0.5" +version = "1.0.18" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "71d301d4193d031abdd79ff7e3dd721168a9572ef3fe51a1517aba235bd8f86e" +checksum = "f3cb5ba0dc43242ce17de99c180e96db90b235b8a9fdc9543c96d2209116bd9f" [[package]] name = "same-file" @@ -482,42 +334,20 @@ dependencies = [ "winapi-util", ] -[[package]] -name = "scopeguard" -version = "1.1.0" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d29ab0c6d3fc0ee92fe66e2d99f700eab17a8d57d1c1d3b748380fb20baa78cd" - -[[package]] -name = "semver" -version = "1.0.4" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "568a8e6258aa33c13358f81fd834adb854c6f7c9468520910a9b1e8fac068012" - [[package]] name = "serde" -version = "1.0.127" +version = "1.0.203" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f03b9878abf6d14e6779d3f24f07b2cfa90352cfec4acc5aab8f1ac7f146fae8" +checksum = "7253ab4de971e72fb7be983802300c30b5a7f0c2e56fab8abfc6a214307c0094" dependencies = [ "serde_derive", ] -[[package]] -name = "serde_cbor" -version = "0.11.1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1e18acfa2f90e8b735b2836ab8d538de304cbb6729a7360729ea5a895d15a622" -dependencies = [ - "half", - "serde", -] - [[package]] name = "serde_derive" -version = "1.0.127" +version = "1.0.203" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "a024926d3432516606328597e0f224a51355a493b49fdd67e9209187cbe55ecc" +checksum = "500cbc0ebeb6f46627f50f3f5811ccf6bf00643be300b4c3eabc0ef55dc5b5ba" dependencies = [ "proc-macro2", "quote", @@ -526,9 +356,9 @@ dependencies = [ [[package]] name = "serde_json" -version = "1.0.66" +version = "1.0.117" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "336b10da19a12ad094b59d870ebde26a45402e5b470add4b5fd03c5048a32127" +checksum = "455182ea6142b14f93f4bc5320a2b31c1f266b66a4a5c858b013302a5d8cbfc3" dependencies = [ "itoa", "ryu", @@ -537,23 +367,20 @@ dependencies = [ [[package]] name = "syn" -version = "1.0.74" +version = "2.0.48" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "1873d832550d4588c3dbc20f01361ab00bfe741048f71e3fecf145a7cc18b29c" +checksum = "0f3531638e407dfc0814761abb7c00a5b54992b849452a0646b7f65c9f770f3f" dependencies = [ "proc-macro2", "quote", - "unicode-xid", + "unicode-ident", ] [[package]] name = "textwrap" -version = "0.11.0" +version = "0.16.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "d326610f408c7a4eb6f51c37c330e496b08506c9457c9d34287ecc38809fb060" -dependencies = [ - "unicode-width", -] +checksum = "222a222a5bfe1bba4a77b45ec488a741b3cb8872e5e499451fd7d0129c9c7c3d" [[package]] name = "tinytemplate" @@ -566,97 +393,26 @@ dependencies = [ ] [[package]] -name = "unicode-width" -version = "0.1.8" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9337591893a19b88d8d87f2cec1e73fad5cdfd10e5a6f349f498ad6ea2ffb1e3" - -[[package]] -name = "unicode-xid" -version = "0.2.2" +name = "unicode-ident" +version = "1.0.12" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "8ccb82d61f80a663efe1f787a51b16b5a51e3314d6ac365b08639f52387b33f3" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" [[package]] name = "walkdir" -version = "2.3.2" +version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "808cf2735cd4b6866113f648b791c6adc5714537bc222d9347bb203386ffda56" +checksum = "29790946404f91d9c5d06f9874efddea1dc06c5efe94541a7d6863108e3a5e4b" dependencies = [ "same-file", - "winapi", "winapi-util", ] [[package]] name = "wasi" -version = "0.9.0+wasi-snapshot-preview1" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "cccddf32554fecc6acb585f82a32a72e28b48f8c4c1883ddfeeeaa96f7d8e519" - -[[package]] -name = "wasm-bindgen" -version = "0.2.75" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b608ecc8f4198fe8680e2ed18eccab5f0cd4caaf3d83516fa5fb2e927fda2586" -dependencies = [ - "cfg-if", - "wasm-bindgen-macro", -] - -[[package]] -name = "wasm-bindgen-backend" -version = "0.2.75" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "580aa3a91a63d23aac5b6b267e2d13cb4f363e31dce6c352fca4752ae12e479f" -dependencies = [ - "bumpalo", - "lazy_static", - "log", - "proc-macro2", - "quote", - "syn", - "wasm-bindgen-shared", -] - -[[package]] -name = "wasm-bindgen-macro" -version = "0.2.75" +version = "0.11.1+wasi-snapshot-preview1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "171ebf0ed9e1458810dfcb31f2e766ad6b3a89dbda42d8901f2b268277e5f09c" -dependencies = [ - "quote", - "wasm-bindgen-macro-support", -] - -[[package]] -name = "wasm-bindgen-macro-support" -version = "0.2.75" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "6c2657dd393f03aa2a659c25c6ae18a13a4048cebd220e147933ea837efc589f" -dependencies = [ - "proc-macro2", - "quote", - "syn", - "wasm-bindgen-backend", - "wasm-bindgen-shared", -] - -[[package]] -name = "wasm-bindgen-shared" -version = "0.2.75" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2e0c4a743a309662d45f4ede961d7afa4ba4131a59a639f29b0069c3798bbcc2" - -[[package]] -name = "web-sys" -version = "0.3.52" -source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "01c70a82d842c9979078c772d4a1344685045f1a5628f677c2b2eab4dd7d2696" -dependencies = [ - "js-sys", - "wasm-bindgen", -] +checksum = "ccf3ec651a847eb01de73ccad15eb7d99f80485de043efb2f370cd654f4ea44b" [[package]] name = "winapi" @@ -676,9 +432,9 @@ checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" [[package]] name = "winapi-util" -version = "0.1.5" +version = "0.1.6" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "70ec6ce85bb158151cae5e5c87f95a8e97d2c0c4b001223f33a334e3ce5de178" +checksum = "f29e6f9198ba0d26b4c9f07dbe6f9ed633e1f3d5b8b414090084349e46a52596" dependencies = [ "winapi", ] diff --git a/Cargo.toml b/Cargo.toml index fbf5d16..70ad0d1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,13 +1,24 @@ [package] name = "rust-lapper" -version = "1.3.0" +version = "2.0.0-beta.1" authors = ["Seth Stadick "] edition = "2018" +rust-version = "1.59" license = "MIT" repository = "https://github.com/sstadick/rust-lapper" description = "A fast and easy interval overlap library" documentation = "https://docs.rs/rust-lapper" readme = "README.md" +include = [ + "/src/**", + "/tests/**", + "/benches/**", + "/examples/**", + "/Cargo.toml", + "/Cargo.lock", + "/LICENSE", + "/README.md", +] categories = ["algorithms", "science"] keywords = ["tree", "interval", "bioinformatics", "ranges", "genomic"] @@ -25,12 +36,13 @@ sort_unstable = [] [dependencies] num-traits = "0.2.12" -serde = { version = "1.0.127", features = ["derive"], optional = true } +serde = { version = "1.0.127", features = ["derive"], optional = true } [dev-dependencies] -criterion = "0.3" -rand = "0.7" -cpu-time = "0.1.0" +# Criterion 0.4 is the newest release compatible with the Rust 1.59 MSRV. +criterion = { version = "0.4", default-features = false, features = ["cargo_bench_support"] } +rand = "0.8" +# Retain bincode 1.x to test compatibility with rust-lapper 1.x data. bincode = "1.3.3" [[bench]] diff --git a/README.md b/README.md index 92f5480..252c57d 100644 --- a/README.md +++ b/README.md @@ -9,23 +9,64 @@ [Documentation](https://docs.rs/rust-lapper) [Crates.io](https://crates.io/crates/rust-lapper) -This is a rust port of Brent Pendersen's -[nim-lapper](https://github.com/brentp/nim-lapper). It has a few notable -differences, mostly that the find and seek methods both return -iterators, so all adaptor methods may be used normally. - -This crate works well for most interval data that does not include very long -intervals that engulf a majority of other intervals. It is still fairly -comparable to other methods. If you absolutely need time guarantees in the -worst case, see [COItres](https://github.com/dcjones/coitrees) and [IITree](https://docs.rs/bio/0.32.0/bio/data_structures/interval_tree/struct.ArrayBackedIntervalTree.html). - -However, on more typical datasets, this crate is between 4-10x faster -than other interval overlap methods. - -It should also be noted that the `count` method is agnostic to data -type, and should be about as fast as it is possible to be on any -dataset. It is an implementation of the [BITS -algorithm](https://academic.oup.com/bioinformatics/article/29/1/1/273289) +rust-lapper 2 is currently in beta. To try it, use `2.0.0-beta.1` +explicitly; Cargo will not select it from a `version = "2"` requirement. + +This was originally a Rust port of Brent Pedersen's +[nim-lapper](https://github.com/brentp/nim-lapper). `find()` and `seek()` return +lazy borrowed iterators in ascending start order, so normal iterator adaptors +work without collecting results first. + +All stored intervals and query ranges use half-open `[start, stop)` semantics. +`Lapper` keeps its intervals sorted by start and builds a fixed 32-interval +block index that can skip regions proven not to overlap. Mixed blocks use NEON +on AArch64, runtime-detected AVX2 on x86-64, and an exact scalar fallback +elsewhere. The same algorithm handles both ordinary and pathological datasets +with long intervals that engulf many shorter intervals. + +The `count()` method uses the +[BITS algorithm](https://academic.oup.com/bioinformatics/article/29/1/1/273289) +to count overlaps with two binary searches. + +## API and algorithm compatibility + +The block index and SIMD backends are private implementation details: there is +no mode flag, alternate query method, or architecture-specific API. Existing +call patterns for `find()`, `seek()`, `count()`, `cov()`, `set_cov()`, +`merge_overlaps()`, `depth()`, `union_and_intersect()`, `union()`, and +`intersect()` retain their range semantics and return types. `count()` remains +the independent BITS implementation; methods that use `find()` or `seek()` +internally automatically share the exact indexed query path. + +| Target | Mixed-block backend | Selection | +|---|---|---| +| AArch64 | 128-bit NEON | Baseline for the architecture | +| x86-64 with AVX2 | 256-bit AVX2 | Runtime detected once per iterator | +| x86-64 without AVX2 | Scalar | Automatic fallback | +| Other architectures | Scalar | Automatic fallback | + +NEON and AVX2 cover `u8`, `i8`, `u16`, `i16`, `u32`, `i32`, `u64`, `i64`, +`usize`, and `isize`. The same block algorithm uses exact scalar masks for +`u128`, `i128`, custom `PrimInt` types, and partial vector tails. + +## Minimum Supported Rust Version + +rust-lapper 2 supports Rust 1.59 and newer. Rust 1.59 is the first stable +release that provides the AArch64 intrinsics used by the NEON query backend. + +Query coordinates must be `'static` so private dispatch code can use `TypeId` +before reinterpreting primitive integer slices for SIMD. This includes every +primitive integer and ordinary owned custom numeric type; it does not require a +`Lapper` value to live for the entire program. Non-primitive `PrimInt` types use +the scalar mask implementation. The bound is the API-breaking change that makes +this a major release. + +## Mutation + +Use `insert()` and `merge_overlaps()` for coordinate or structural changes so +the private query index is rebuilt. Directly changing `Lapper::intervals` +coordinates or length leaves derived metadata stale; changing payload values is +safe. ## Serde Support @@ -33,63 +74,22 @@ algorithm](https://academic.oup.com/bioinformatics/article/29/1/1/273289) ```toml [dependencies] -rust-lapper = { version = "*", features = ["with_serde"] } +rust-lapper = { version = "2.0.0-beta.1", features = ["with_serde"] } ``` See `examples/serde.rs` for a brief example. ## Benchmarks -Benchmarking interval tree-ish datastructures is hard -Please see the -[interval_bakeoff](https://github.com/sstadick/interval_bakeoff) project -for details on how the benchmarks were run... It's not fully baked yet -though, and is finiky to run. - -Command to run: - -``` -./target/release/interval_bakeoff fake -a -l RustLapper -l -RustBio -l NestedInterval -n50000 -u100000 - -# This equates to the following params: -# num_intervals 50000 -# universe_size 100000 -# min_interval_size 500 -# max_interval_size 80000 -# add_large_span true (universe spanning) -``` - -Set A / b Creation Times - -| crate/method | A time | B time | -| ---------------- | -------- | -------- | -| rust_lapper | 15.625ms | 31.25ms | -| nested_intervals | 15.625ms | 15.625ms | -| bio | 15.625ms | 31.25ms | - -100% hit rate (A vs A) - -| crate/method | mean time | intersection | -| ---------------------------------- | ---------- | ------------ | -| rust_lapper/find | 4.78125s | 1469068763 | -| rust_lapper/count | 15.625ms | 1469068763 | -| nested_intervals/query_overlapping | 157.4375s | 1469068763 | -| bio/find | 33.296875s | 1469068763 | +The retained v2 release measurements, raw samples, compiler flags, and pinned +competitor revisions live in +[lapper_bakeoff](https://github.com/sstadick/lapper_bakeoff/tree/main/results/avx2-2026-07-28). +On an AMD Ryzen 9 3950X with AVX2, the new implementation improved total time +over rust-lapper 1.3.0 by 37.30%, 98.89%, and 34.64% on the three retained +article cases. All implementations returned identical overlap counts. - -Sub 100% hit rate (A vs B) - -| crate/method | mean time | intersection | -| ---------------------------------- | ---------- | ------------ | -| rust_lapper/find | 531.25ms | 176488436 | -| rust_lapper/count | 15.625ms | 176488436 | -| nested_intervals/query_overlapping | 11.109375s | 196090092 | -| bio/find | 4.3125s | 176488436 | - -[nested_intervals](https://docs.rs/nested_intervals/0.2.0/nested_intervals/) -[rust-bio](https://docs.rs/bio/0.28.2/bio/) -*Note that rust-bio has a new interval tree structure which should be faster than what is shown here* +Benchmark results are workload- and hardware-specific; use the linked harness +and raw data when making comparisons. ## Example @@ -104,7 +104,7 @@ fn main() { start: 70, stop: 120, val: 0, - }, // max_len = 50 + }, // a long interval Iv { start: 10, stop: 15, @@ -155,9 +155,8 @@ fn main() { // make lapper structure let mut lapper = Lapper::new(data); - // Iterator based find to extract all intervals that overlap 6..7 - // If your queries are coming in start sorted order, use the seek method to retain a cursor for - // a big speedup. + // Find every interval that overlaps [11, 15). + // For queries in nondecreasing start order, seek() can reuse a caller-owned cursor. assert_eq!( lapper.find(11, 15).collect::>(), vec![ @@ -183,9 +182,10 @@ fn main() { }, // overlap end ] ); + assert_eq!(lapper.count(11, 15), 4); - // Merge overlaping regions within the lapper to simplifiy and speed up quries that only depend - // on 'any + // Merge overlapping regions to simplify queries that only depend on whether + // any interval overlaps. lapper.merge_overlaps(); assert_eq!( lapper.find(11, 15).collect::>(), @@ -196,10 +196,10 @@ fn main() { },] ); - // Get the number of positions covered by the lapper tree: + // Get the number of positions covered by the interval collection. assert_eq!(lapper.cov(), 73); - // Get the union and intersect of two different lapper trees + // Get the union and intersection lengths of two interval collections. let data = vec![ Iv { start: 5, @@ -229,10 +229,13 @@ fn main() { ## Release Notes +- `2.0.0-beta.1`: Begin the rust-lapper 2 beta with a portable SIMD block + index, signed coordinates, Rust 1.59 as the MSRV, and the `I: 'static` + coordinate bound. - `1.3.0`: Add the `sort_unstable` feature flag for allocation-sensitive sorting thanks to @jameslkingsley. - `1.1.0`: Added insert functionality thanks to @zaporter -- `0.4.0`: Addition of the BITS count algorithm. -- `0.4.2`: Bugfix in to update starts/stops vectors when overlaps merged -- `0.4.3`: Remove leftover print statement -- `0.5.0`: Make Interval start/stop generic - `1.0.0`: Add serde support via the `with_serde` feature flag +- `0.5.0`: Make Interval start/stop generic +- `0.4.3`: Remove leftover print statement +- `0.4.2`: Bugfix in to update starts/stops vectors when overlaps merged +- `0.4.0`: Addition of the BITS count algorithm. diff --git a/benches/lapper_benchmark.rs b/benches/lapper_benchmark.rs index a52a872..b5b4e8a 100644 --- a/benches/lapper_benchmark.rs +++ b/benches/lapper_benchmark.rs @@ -14,7 +14,7 @@ type Iv = Interval; fn randomi(imin: u32, imax: u32) -> u32 { let mut rng = rand::thread_rng(); - imin + rng.gen_range(0, imax - imin) + imin + rng.gen_range(0..imax - imin) } fn make_random(n: usize, range_max: u32, size_min: u32, size_max: u32) -> Vec { @@ -52,8 +52,8 @@ fn make_random_seeded( let mut rng = StdRng::seed_from_u64(seed); let mut result = Vec::with_capacity(n); for _ in 0..n { - let s = rng.gen_range(0, range_max); - let e = s + rng.gen_range(size_min, size_max); + let s = rng.gen_range(0..range_max); + let e = s + rng.gen_range(size_min..size_max); result.push(Interval { start: s, stop: e, diff --git a/devlog.md b/devlog.md index 4098081..3419787 100644 --- a/devlog.md +++ b/devlog.md @@ -1,7 +1,13 @@ +# Historical development log + +> These notes preserve the 2019 development history and are not current API or +> release guidance. See `README.md` and the generated crate documentation for +> rust-lapper 2 behavior. + ## 9/20/19 - Maybe add max_ends list to the Lapper. Then binary search for upper bounds like AIList does, and iterate backwards checking against the - max_end to deteremine overlap... It could make the hopper actually + max_end to determine overlap... It could make the hopper actually work as well since binary searches in a local area would be useful ## 9/3/19 @@ -13,7 +19,7 @@ - Corrected the range on the IterFind to properly stop looking once stop has been hit. - Added benchmarks for find and seek -- Sped up the intesection method when we know that overlaps have been +- Sped up the intersection method when we know that overlaps have been merged. - Figure out the possible speedup with miss counting!! - Fix up docs to include features list / make less block of text @@ -21,7 +27,7 @@ - Corrected the intersection_and_union function to work regardless of the comparison direction. It's a little slower now and in need of benchmarking. I think there is a gain to be had if we know that both - lappers have had overlaps merged, I'm just now sure how to make it + lappers have had overlaps merged, I'm just not sure how to make it work yet. - Copy the version of intersect_and_union from chromcomp, it works and is faster @@ -30,8 +36,8 @@ ## 8/31/19 - Moved a bunch of functionality in, like merge overlaps and detect intersects and coverage calculation. -- I would like to add some whole Lapper functions like a union an - intersect. -- did this... but needs fine tuning +- I would like to add some whole Lapper functions like a union and + intersection. -- did this... but needs fine tuning - I would like to add benchmarks and identify the slow stuff - Something is now really slow I think. Like maybe the coverage calculation. Bench it and see. I also think the checked_sub might be diff --git a/examples/ex1.rs b/examples/ex1.rs index 80d0cd0..b68750b 100644 --- a/examples/ex1.rs +++ b/examples/ex1.rs @@ -8,7 +8,7 @@ fn main() { start: 70, stop: 120, val: 0, - }, // max_len = 50 + }, // a long interval Iv { start: 10, stop: 15, @@ -59,9 +59,8 @@ fn main() { // make lapper structure let mut lapper = Lapper::new(data); - // Iterator based find to extract all intervals that overlap 6..7 - // If your queries are coming in start sorted order, use the seek method to retain a cursor for - // a big speedup. + // Find every interval that overlaps [11, 15). + // For queries in nondecreasing start order, seek() can reuse a caller-owned cursor. assert_eq!( lapper.find(11, 15).collect::>(), vec![ @@ -87,9 +86,10 @@ fn main() { }, // overlap end ] ); + assert_eq!(lapper.count(11, 15), 4); - // Merge overlaping regions within the lapper to simplifiy and speed up quries that only depend - // on 'any + // Merge overlapping regions to simplify queries that only depend on whether + // any interval overlaps. lapper.merge_overlaps(); assert_eq!( lapper.find(11, 15).collect::>(), @@ -100,10 +100,10 @@ fn main() { },] ); - // Get the number of positions covered by the lapper tree: + // Get the number of positions covered by the interval collection. assert_eq!(lapper.cov(), 73); - // Get the union and intersect of two different lapper trees + // Get the union and intersection lengths of two interval collections. let data = vec![ Iv { start: 5, diff --git a/examples/serde.rs b/examples/serde.rs index a5b6def..342750b 100644 --- a/examples/serde.rs +++ b/examples/serde.rs @@ -8,7 +8,7 @@ fn main() { start: 70, stop: 120, val: 0, - }, // max_len = 50 + }, // a long interval Iv { start: 10, stop: 15, @@ -59,9 +59,8 @@ fn main() { // make lapper structure let mut lapper = Lapper::new(data); - // Iterator based find to extract all intervals that overlap 6..7 - // If your queries are coming in start sorted order, use the seek method to retain a cursor for - // a big speedup. + // Find every interval that overlaps [11, 15). + // For queries in nondecreasing start order, seek() can reuse a caller-owned cursor. assert_eq!( lapper.find(11, 15).collect::>(), vec![ @@ -87,9 +86,10 @@ fn main() { }, // overlap end ] ); + assert_eq!(lapper.count(11, 15), 4); - // Merge overlaping regions within the lapper to simplifiy and speed up quries that only depend - // on 'any + // Merge overlapping regions to simplify queries that only depend on whether + // any interval overlaps. lapper.merge_overlaps(); assert_eq!( lapper.find(11, 15).collect::>(), @@ -100,10 +100,10 @@ fn main() { },] ); - // Get the number of positions covered by the lapper tree: + // Get the number of positions covered by the interval collection. assert_eq!(lapper.cov(), 73); - // Get the union and intersect of two different lapper trees + // Get the union and intersection lengths of two interval collections. let data = vec![ Iv { start: 5, @@ -130,6 +130,7 @@ fn main() { let encoded = bincode::serialize(&lapper).unwrap(); let decoded: Lapper = bincode::deserialize(&encoded[..]).unwrap(); - dbg!(lapper); - dbg!(decoded); + assert_eq!(decoded.intervals, lapper.intervals); + assert_eq!(decoded.cov(), lapper.cov()); + assert_eq!(decoded.count(11, 15), lapper.count(11, 15)); } diff --git a/plans/simd/AARCH64_U32_HAND_TYPE.md b/plans/simd/AARCH64_U32_HAND_TYPE.md new file mode 100644 index 0000000..9343a08 --- /dev/null +++ b/plans/simd/AARCH64_U32_HAND_TYPE.md @@ -0,0 +1,1528 @@ +# Hand-type the whole forward index: original Lapper to portable SIMD + +> This is a historical, checkpoint-based implementation course. Temporary APIs +> and intermediate branches shown here are teaching steps; [`README.md`](../../README.md) +> and the generated crate documentation describe the final public API. + +This is the single end-to-end exercise for understanding the worked +rust-lapper contender. It starts with the original start-sorted linear iterator, +adds the first per-interval skip index, replaces it with the 32-entry block +index, adds exact block routes and a saved scalar mask, types the AArch64 NEON +path one operation at a time, and finishes with the portable, mutation-safe +implementation. + +Do not begin with the SIMD intrinsics. SIMD is only a faster implementation of +one integer-valued function. The data structure, skip proofs, iterator state, +and output order must already make sense before that function is replaced. + +Every checkpoint has the same four gates: + +1. **Type it:** the smallest useful source change. +2. **State the invariant:** the fact that makes the change correct. +3. **Work an example:** calculate the result without running Rust. +4. **Pass the gate:** compare with a preserved branch and run an exact-order test. + +## Set up a practice worktree + +Keep this worked branch open as the guide and create a second worktree from the +original release state: + +```bash +cd /Users/sethstadick/dev/super_intervals_more/bakeoff/sandboxes/rust-lapper-contenders + +git worktree add \ + -b practice/hand-typed-forward-index \ + ../rust-lapper-hand-typed \ + baseline/v1.3.0 + +cd ../rust-lapper-hand-typed +``` + +The preserved answer checkpoints are: + +| Checkpoint | Branch | What it isolates | +|---|---|---| +| Original | `baseline/v1.3.0` | Start-sorted scalar Lapper | +| First index | `contender/per-interval` | One next-greater-end link per interval | +| Block index | `contender/block-32` | One summary and link per 32 intervals | +| First mask | `contender/simd-mask` | Exact routes, saved mask, first NEON implementation | +| Weighted mask | `contender/simd-mask-addv` | Four truth lanes reduced to four bits | +| Paired loads | `contender/simd-mask-addv-ld1x2` | Eight values loaded per loop | +| Eight-lane pack | `contender/simd-mask-narrow8-ld1x2` | Eight truths reduced once | +| Bounds audit | `contender/simd-mask-narrow8-unchecked-index` | Proof-scoped unchecked sidecars | +| Narrow ARM answer | `tutorial/aarch64-u32-hand-typed` | The ARM-only portion without portable integration | +| Full answer | `worked/portable-simd-index` | Normal API, all primitive types, mutation, serde, and fallbacks | + +Use this as the progress sheet: + +- [ ] 0. Install the forward brute-force oracle. +- [ ] 1. Explain original `max_len` search and its pathological case. +- [ ] 2. Build and prove per-interval next-greater links. +- [ ] 3. Replace them with 32-entry block summaries and links. +- [ ] 4. Binary-search block prefix maxima. +- [ ] 5. Derive the exact skip, dense, and mixed routes. +- [ ] 6. Save and drain a scalar block mask. +- [ ] 7. Reproduce that mask with four-lane NEON comparisons. +- [ ] 8. Pack four truth lanes with weighted horizontal addition. +- [ ] 9. Pair loads, narrow truth, and pack eight lanes once. +- [ ] 10. Find the intended instructions in release assembly. +- [ ] 11. Move the mask state behind normal `find()` and `seek()`. +- [ ] 12. Add exact type and CPU dispatch with scalar fallback. +- [ ] 13. Rebuild all derived state after mutation and deserialize. +- [ ] 14. Repair the scalar algorithms for signed coordinates. +- [ ] 15. Remove only bounds checks covered by private invariants. + +At any point, inspect the historical delta rather than copying the whole file: + +```bash +git diff baseline/v1.3.0..contender/per-interval -- src/lib.rs +git diff contender/per-interval..contender/block-32 -- src/lib.rs +git diff contender/block-32..contender/simd-mask -- src/lib.rs tests/block_index.rs +``` + +The historical index and SIMD branches through checkpoint 12 were immutable +query experiments: their older mutation methods do not rebuild every new +sidecar. Use `new()` plus queries while learning those stages. Checkpoint 13 is +where the exercise becomes mutation-safe. + +## The contract that never changes + +Intervals and queries are half-open: + +```text +interval = [interval.start, interval.stop) +query = [query.start, query.stop) +``` + +They overlap only when both strict comparisons pass: + +```rust +interval.start < query_stop && interval.stop > query_start +``` + +Equality means no overlap. An interval ending at `query_start` is to the left; +an interval starting at `query_stop` is to the right. + +The observable identity of this design is also fixed: + +- `intervals` is the one canonical vector; +- it is sorted by `(start, stop)` during construction; +- queries yield borrowed entries from that vector; +- results appear in increasing vector position, hence start order; and +- no index may omit a true overlap or invent a false one. + +The index may skip, classify, or compare many entries together. It may not +reorder the canonical vector or return results backward. + +## Checkpoint 0: make brute force the executable specification + +Before adding an index, add an integration test that calculates the answer in the +most obvious way and compares values in order. This is more important than a +count-only test: two iterators can return the same count while returning the +wrong entries or order. + +Create `tests/block_index.rs`: + +```rust +use rust_lapper::{Interval, Lapper}; + +#[test] +fn block_index_matches_forward_brute_force() { + let mut state = 0x1234_5678_u64; + let mut next = || { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1); + (state >> 32) as u32 + }; + + let mut intervals = vec![Interval { + start: 0, + stop: 1_000_000, + val: 0usize, + }]; + for value in 1..5000 { + let start = next() % 1_000_000; + let len = 1 + next() % 1000; + intervals.push(Interval { + start, + stop: start + len, + val: value, + }); + } + let lapper = Lapper::new(intervals); + + for _ in 0..20_000 { + let start = next() % 1_000_000; + let stop = start + 1 + next() % 2000; + let got: Vec<_> = lapper.find(start, stop).map(|iv| iv.val).collect(); + let expected: Vec<_> = lapper + .intervals + .iter() + .filter(|iv| iv.start < stop && iv.stop > start) + .map(|iv| iv.val) + .collect(); + assert_eq!(got, expected, "query {start}..{stop}"); + } +} +``` + +### Invariant + +The expected vector is produced by filtering the canonical start-sorted vector +without skipping or reordering. It is the oracle for every later iterator. + +### Work it by hand + +For intervals `[0, 3)`, `[2, 8)`, `[5, 7)`, and `[9, 12)`, query `[3, 6)`: + +```text +[0, 3): 0 < 6, but 3 > 3 is false +[2, 8): 2 < 6 and 8 > 3 -> hit +[5, 7): 5 < 6 and 7 > 3 -> hit +[9, 12): 9 < 6 is false; all later starts also fail + +answer: the second interval, then the third interval +``` + +### Pass the gate + +```bash +cargo test --test block_index +``` + +Do not proceed until this passes on `baseline/v1.3.0`. + +## Checkpoint 1: understand the original search + +Original Lapper stores: + +```text +intervals canonical intervals sorted by (start, stop) +starts starts sorted globally +stops stops sorted globally, no longer paired with starts +max_len length of the longest interval +``` + +`find(query_start, query_stop)` computes a conservative earliest possible start: + +```rust +let earliest_start = query_start + .checked_sub(&self.max_len) + .unwrap_or_else(zero::); +let off = Self::lower_bound(earliest_start, &self.intervals); +``` + +Why it is safe: no interval starting before `query_start - max_len` can be long +enough to reach `query_start`. From `off`, the iterator checks intervals one by +one and stops forever when `interval.start >= query_stop`. + +### Invariant + +Sorted starts prove the right-hand stopping condition. `max_len` proves the +left-hand starting condition. + +### Work it by hand + +If `query_start = 1_000` and `max_len = 40`, no interval starting before `960` +can reach the query. But one extremely long interval makes `max_len` enormous, +pushes `earliest_start` back toward zero, and forces nearly every query to scan a +large prefix. That is the pathological case the new index must repair. + +### Pass the gate + +Be able to answer these before continuing: + +```text +Why can the iterator stop on interval.start >= query_stop? +Why can it not stop merely because one interval.stop <= query_start? +Why does one global long interval make max_len a weak lower bound? +``` + +## Checkpoint 2: add one next-greater link per interval + +The first index asks: if interval `i` ends too early, which later interval is the +first one whose end is strictly greater? + +Add one sidecar: + +```rust +jump_index: Vec, +``` + +Build it in linear time with a monotonic stack, scanning right to left: + +```rust +let interval_count = intervals.len(); +let mut jump_index = vec![interval_count; interval_count]; +let mut stack = Vec::::new(); + +for i in (0..interval_count).rev() { + while stack + .last() + .is_some_and(|j| intervals[*j].stop <= intervals[i].stop) + { + stack.pop(); + } + jump_index[i] = stack.last().copied().unwrap_or(interval_count); + stack.push(i); +} +``` + +Why the stack gives the *nearest* greater end: + +- before processing `i`, the stack contains undominated positions to its right; +- viewed from the top outward, their ends are strictly increasing; +- anything popped has an end no greater than `end[i]` and is farther right, so + `i` is both closer and at least as capable of reaching a future query; +- after those dominated entries are removed, the top is the closest surviving + position with an end strictly greater than `end[i]`. + +That domination argument is why construction is linear: every position is +pushed once and popped at most once. + +Then replace the scalar miss step with the link: + +```rust +while self.off < self.inner.intervals.len() { + let interval = &self.inner.intervals[self.off]; + + if interval.start >= self.stop { + break; + } + if interval.stop > self.start { + self.off += 1; + return Some(interval); + } + + self.off = self.inner.jump_index[self.off]; +} +``` + +### Invariant + +Suppose `jump_index[i] = j`. Every index strictly between `i` and `j` has +`stop <= intervals[i].stop`. The link is used only after proving +`intervals[i].stop <= query_start`, so every skipped interval also ends at or +before the query and cannot overlap. + +The link points right, so forward order is preserved. + +### Work it by hand + +For start-ordered end values: + +```text +index: 0 1 2 3 4 +end: 7 3 5 2 9 +link: 4 2 4 4 sentinel +``` + +With `query_start = 6`, a miss at index 1 follows `1 -> 2 -> 4`: + +```text +end[1] = 3 <= 6 +end[2] = 5 <= 6 +end[4] = 9 > 6 +``` + +Index 3 is skipped by `2 -> 4` because its end `2` is no greater than end `5`. + +### Pass the gate + +```bash +cargo test --test block_index +git diff contender/per-interval -- src/lib.rs +``` + +This historical branch is an immutable-query experiment. Its `insert()` and +`merge_overlaps()` paths do not rebuild the new link array. Do not treat that as +production-correct; mutation is repaired in checkpoint 13. + +## Checkpoint 3: amortize the index over fixed blocks + +The per-interval index costs one `usize` per interval and performs irregular link +work for individual misses. Replace it with one summary for each 32 consecutive +intervals: + +```rust +const INDEX_BLOCK_SIZE: usize = 32; + +block_index: Vec, +block_max_ends: Vec, +block_prefix_max_ends: Vec, +``` + +First calculate each block's maximum end: + +```rust +let mut block_max_ends = Vec::new(); +for block in intervals.chunks(INDEX_BLOCK_SIZE) { + let mut max_end = block[0].stop; + for interval in &block[1..] { + max_end = std::cmp::max(max_end, interval.stop); + } + block_max_ends.push(max_end); +} +``` + +Run the same next-greater construction over block maxima: + +```rust +let block_count = block_max_ends.len(); +let mut block_index = vec![block_count; block_count]; +let mut stack = Vec::::new(); + +for block in (0..block_count).rev() { + while stack + .last() + .is_some_and(|next| block_max_ends[*next] <= block_max_ends[block]) + { + stack.pop(); + } + block_index[block] = stack.last().copied().unwrap_or(block_count); + stack.push(block); +} +``` + +At a block boundary, skip only when the entire block is an exact miss: + +```rust +let block = block_start / INDEX_BLOCK_SIZE; +if self.inner.block_max_ends[block] <= self.start { + self.next_block_start = + self.inner.block_index[block] * INDEX_BLOCK_SIZE; + continue; +} +``` + +### Invariant + +`block_max_end <= query_start` proves every end in the block fails the strict +`end > query_start` overlap test. If the link jumps from block `b` to block `j`, +all intermediate block maxima are no greater than `max_end[b]`, so they also +fail. + +### Work it by hand + +Use a paper-only block size of four. The real implementation remains 32. + +| Block | Start-ordered intervals | Minimum end | Maximum end | +|---|---|---:|---:| +| 0 | `[0,4) [2,30) [5,7) [8,9)` | 4 | 30 | +| 1 | `[10,12) [13,14) [15,18) [19,21)` | 12 | 21 | +| 2 | `[22,60) [24,25) [27,35) [40,44)` | 25 | 60 | + +The block maxima are `[30, 21, 60]`, so the links are `[2, 2, sentinel]`. +For query `[23, 28)`, block 1 has `max_end = 21 <= 23` and can jump directly to +block 2. + +### Pass the gate + +```bash +cargo test --test block_index +git diff contender/block-32 -- src/lib.rs +``` + +Explain why a block maximum can prove an all-miss block but cannot prove which +individual intervals hit inside a candidate block. + +## Checkpoint 4: add the prefix maximum entry search + +`find()` should not enter block zero for every query. Build a monotonically +nondecreasing prefix over block maxima: + +```rust +let mut block_prefix_max_ends = block_max_ends.clone(); +for block in 1..block_prefix_max_ends.len() { + block_prefix_max_ends[block] = std::cmp::max( + block_prefix_max_ends[block - 1], + block_prefix_max_ends[block], + ); +} +``` + +Binary-search it at query creation: + +```rust +let first_block = self + .block_prefix_max_ends + .partition_point(|max_end| *max_end <= query_start); +let next_block_start = first_block * INDEX_BLOCK_SIZE; +``` + +### Invariant + +If `prefix_max[b] <= query_start`, every interval in every block through `b` +ends too early. The first prefix value greater than `query_start` is therefore +the first block that could contain an overlap. + +### Work it by hand + +For block maxima `[30, 21, 60]`, the prefix maxima are `[30, 30, 60]`. + +```text +query_start = 23 -> partition point 0 -> block 0 may contain [2,30) +query_start = 31 -> partition point 2 -> blocks 0 and 1 are impossible +query_start = 61 -> partition point 3 -> no block can overlap +``` + +### Pass the gate + +Add assertions for those three partition points or write them on paper before +running: + +```bash +cargo test --test block_index +``` + +Notice that final `find()` no longer depends on global `max_len`; `seek()` still +uses it to advance a cursor for sorted query streams. + +## Checkpoint 5: add exact candidate-block routes + +Maximum ends prove an all-miss route. Add the complementary minimum end: + +```rust +block_min_ends: Vec, +``` + +Calculate minimum and maximum together. Also preserve ends in start order: + +```rust +let (starts, mut stops): (Vec<_>, Vec<_>) = + intervals.iter().map(|iv| (iv.start, iv.stop)).unzip(); +let stops_by_start = stops.clone(); +stops.sort(); +``` + +These two end arrays are intentionally different: + +```text +stops globally sorted; used by count() +stops_by_start lane i belongs to starts[i]; used by block masks +``` + +Every candidate block now has three exact routes: + +```text +max_end <= query_start + -> all miss; follow the forward block link + +min_end > query_start + -> every end passes; return only starts < query_stop + +otherwise + -> end results are mixed; calculate every lane exactly +``` + +For the all-ends-pass route, sorted starts identify one dense prefix: + +```rust +let starts = &self.inner.starts[block_start..block_end]; +let active_len = starts.partition_point(|start| *start < self.stop); +self.dense_next = block_start; +self.dense_end = block_start + active_len; +``` + +### Invariant + +`min_end > query_start` proves the end half of overlap for every lane. It does +not prove the start half. The active prefix is still needed because starts at or +after `query_stop` do not overlap. + +This answers the active-prefix question precisely: + +- dense route: yes, find the active start prefix; +- mixed route: no separate prefix is needed once the mask checks both predicates; +- whole query: stop forever when the first start of a block is at or beyond + `query_stop`. + +### Work it by hand + +Return to query `[23, 28)` in the three toy blocks: + +```text +block 0: min 4 <= 23 < max 30 -> mixed +block 1: max 21 <= 23 -> jump +block 2: min 25 > 23 -> every end passes +``` + +Block 2 starts are `[22, 24, 27, 40]`. Its active prefix is the first three +lanes because `22`, `24`, and `27` are less than `28`; `40` is not. + +### Pass the gate + +Before adding SIMD, write down why both inequalities are strict: + +```text +max_end == query_start -> all miss +min_end == query_start -> not all ends pass +start == query_stop -> lane misses +``` + +## Checkpoint 6: save one scalar block mask + +Add a temporary `find_block_mask()` iterator while retaining normal `find()` as +the control. This is scaffolding for the experiment, not a permanent API or a +mode switch. + +Give the iterator memory: + +```rust +next_block_start: usize, +mask_block_start: usize, +mask: u32, +dense_next: usize, +dense_end: usize, +start: u32, +stop: u32, +``` + +Start with a scalar mask that is obviously equivalent to brute force: + +```rust +fn overlap_mask( + starts: &[u32], + stops: &[u32], + query_start: u32, + query_stop: u32, +) -> u32 { + let mut mask = 0_u32; + for lane in 0..starts.len() { + if stops[lane] > query_start && starts[lane] < query_stop { + mask |= 1 << lane; + } + } + mask +} +``` + +Drain the lowest saved bit across later `next()` calls: + +```rust +if self.mask != 0 { + let lane = self.mask.trailing_zeros() as usize; + self.mask &= self.mask - 1; + return Some(&self.inner.intervals[self.mask_block_start + lane]); +} +``` + +The full safe state machine is: + +```rust +loop { + if self.dense_next < self.dense_end { + let index = self.dense_next; + self.dense_next += 1; + return Some(&self.inner.intervals[index]); + } + + if self.mask != 0 { + let lane = self.mask.trailing_zeros() as usize; + self.mask &= self.mask - 1; + return Some(&self.inner.intervals[self.mask_block_start + lane]); + } + + let block_start = self.next_block_start; + if block_start >= self.inner.starts.len() + || self.inner.starts[block_start] >= self.stop + { + return None; + } + + let block = block_start / INDEX_BLOCK_SIZE; + if self.inner.block_max_ends[block] <= self.start { + self.next_block_start = + self.inner.block_index[block] * INDEX_BLOCK_SIZE; + continue; + } + + let block_end = + (block_start + INDEX_BLOCK_SIZE).min(self.inner.starts.len()); + + if self.inner.block_min_ends[block] > self.start { + let starts = &self.inner.starts[block_start..block_end]; + let active_len = + starts.partition_point(|lane_start| *lane_start < self.stop); + self.dense_next = block_start; + self.dense_end = block_start + active_len; + self.next_block_start = if active_len == starts.len() { + block_end + } else { + self.inner.starts.len() + }; + continue; + } + + self.mask_block_start = block_start; + self.next_block_start = block_end; + self.mask = overlap_mask( + &self.inner.starts[block_start..block_end], + &self.inner.stops_by_start[block_start..block_end], + self.start, + self.stop, + ); +} +``` + +### Invariant + +Bit `k` describes lane `k` in one block. `trailing_zeros()` selects the lowest +remaining lane, and `mask & (mask - 1)` erases exactly that bit. Blocks only move +right. Therefore output remains in canonical vector order. + +### Work it by hand + +For toy block 0 and query `[23, 28)`: + +```text +lane 0 [0,4) -> 0 +lane 1 [2,30) -> 1 +lane 2 [5,7) -> 0 +lane 3 [8,9) -> 0 + +mask = 0b0010 +trailing_zeros(mask) = 1 +next() returns canonical index block_start + 1 +``` + +Across all blocks the result indices are `[1, 8, 9, 10]`, still in start order. + +### Pass the gate + +Temporarily extend the oracle test: + +```rust +let got_mask: Vec<_> = lapper + .find_block_mask(start, stop) + .map(|iv| iv.val) + .collect(); +assert_eq!(got_mask, expected); +``` + +Then run: + +```bash +cargo test --test block_index +``` + +Do not type an intrinsic until the scalar mask passes. + +## Checkpoint 7: compute four `u32` lanes with NEON + +On AArch64, one 128-bit NEON register holds four `u32` values. Inside an +`unsafe` block, broadcast the query bounds once, then load and compare four +paired intervals: + +```rust +use std::arch::aarch64::*; + +let query_start_v = vdupq_n_u32(query_start); +let query_stop_v = vdupq_n_u32(query_stop); + +let lane_starts = vld1q_u32(starts.as_ptr().add(lane)); +let lane_stops = vld1q_u32(stops.as_ptr().add(lane)); + +let overlaps = vandq_u32( + vcgtq_u32(lane_stops, query_start_v), + vcgtq_u32(query_stop_v, lane_starts), +); +``` + +Read that literally: + +```text +lane_stops > query_start +query_stop > lane_starts +``` + +NEON comparison results are not scalar booleans. A true `u32` lane contains +`0xffff_ffff`; a false lane contains zero. + +For the first version, shift truth down to zero/one and extract each lane: + +```rust +let bits = vshrq_n_u32::<31>(overlaps); +mask |= vgetq_lane_u32::<0>(bits) << lane; +mask |= vgetq_lane_u32::<1>(bits) << (lane + 1); +mask |= vgetq_lane_u32::<2>(bits) << (lane + 2); +mask |= vgetq_lane_u32::<3>(bits) << (lane + 3); +``` + +Keep a scalar tail for a final block whose length is not divisible by four. + +### Invariant + +The SIMD function must return the exact same `u32` as checkpoint 6. It does not +get a new definition of overlap, block membership, or order. + +### Work it by hand + +Truth lanes `[true, false, true, true]` become: + +```text +[0xffff_ffff, 0, 0xffff_ffff, 0xffff_ffff] +shift by 31 -> [1, 0, 1, 1] +packed mask -> 0b1101 +``` + +### Pass the gate + +```bash +cargo test --test block_index +git diff contender/simd-mask -- src/lib.rs tests/block_index.rs +``` + +On non-AArch64, retain the scalar mask under `cfg` so the branch still compiles. + +## Checkpoint 8: replace lane extraction with weighted bits + +Each true lane is already all ones. Give each lane the numeric value of its final +bit: + +```rust +let weights = vld1q_u32([1_u32, 2, 4, 8].as_ptr()); +let bits = vaddvq_u32(vandq_u32(overlaps, weights)); +mask |= bits << lane; +``` + +This is the complete explanation of the apparent magic: + +```text +truths: [all ones, zero, all ones, all ones] +weights: [1, 2, 4, 8 ] +after AND: [1, 0, 4, 8 ] +horizontal +: 13 +binary 13: 0b1101 +``` + +The sum is a mask, not a count. Powers of two cannot carry into each other +because each weight appears at most once. + +### Invariant + +Lane `k` contributes exactly `1 << k` if and only if it overlaps. + +### Work it by hand + +Calculate these without Rust: + +```text +[false, true, false, true] -> 2 + 8 = 10 -> 0b1010 +[true, true, true, true] -> 1 + 2 + 4 + 8 = 15 -> 0b1111 +[false, false, false, false] -> 0 +``` + +### Pass the gate + +```bash +cargo test --test block_index +git diff contender/simd-mask-addv -- src/lib.rs +``` + +## Checkpoint 9: load and pack eight lanes + +First load two adjacent NEON registers at once: + +```rust +let lane_starts = vld1q_u32_x2(starts.as_ptr().add(lane)); +let lane_stops = vld1q_u32_x2(stops.as_ptr().add(lane)); +``` + +Compare `.0` and `.1` independently. The intermediate paired-load checkpoint +uses two four-lane reductions: + +```rust +let bits0 = vaddvq_u32(vandq_u32(overlaps0, weights4)); +let bits1 = vaddvq_u32(vandq_u32(overlaps1, weights4)); +mask |= (bits0 | (bits1 << 4)) << lane; +``` + +Then narrow the two all-ones/all-zero vectors into eight `u16` truth lanes and +reduce only once: + +```rust +let eight_truths = vcombine_u16( + vmovn_u32(overlaps0), + vmovn_u32(overlaps1), +); +let weights8 = vld1q_u16( + [1_u16, 2, 4, 8, 16, 32, 64, 128].as_ptr(), +); +let bits = vaddvq_u16(vandq_u16(eight_truths, weights8)); +mask |= u32::from(bits) << lane; +``` + +Narrowing preserves truth: the low 16 bits of all ones are still all ones, and +zero stays zero. + +Retain this loop structure: + +```text +while at least 8 lanes remain -> paired eight-lane path +while at least 4 lanes remain -> single four-lane path +while any lanes remain -> scalar tail +``` + +### Invariant + +The eight-lane result owns eight consecutive mask bits beginning at `lane`. +Neither paired loads nor narrowing changes which interval belongs to a lane. + +### Work it by hand + +For truth lanes: + +```text +[true, false, true, true, false, false, true, false] +``` + +the kept weights are `[1, 0, 4, 8, 0, 0, 64, 0]`, whose sum is `77`, or +`0b0100_1101`. + +### Pass the gate + +```bash +cargo test --test block_index +git diff contender/simd-mask-addv-ld1x2 -- src/lib.rs +git diff contender/simd-mask-narrow8-ld1x2 -- src/lib.rs +``` + +## Checkpoint 10: inspect the generated AArch64 instructions + +Do not add handwritten assembly because the intrinsic source looks verbose. +Inspect what LLVM actually emitted: + +```bash +RUSTFLAGS='-C target-cpu=native' \ + cargo rustc --release --lib -- --emit=asm + +rg -n 'ldp|uzp1|addv' target/release/deps/rust_lapper-*.s +``` + +On the measured Apple M3 path: + +```text +ldp q..., q... paired vector loads +uzp1.8h narrow/join the relevant halfwords +addv.8h horizontal weighted reduction +``` + +### Invariant + +Intrinsics define semantics; assembly inspection tells you the final cost. Add +inline assembly only when a specific unwanted instruction remains and a measured +replacement wins. That condition was not met here. + +### Work it by hand + +For one eight-lane group, account for the conceptual work before reading the +assembly: + +```text +two vector registers of starts +two vector registers of stops +two end comparisons +two start comparisons +two ANDs +one narrow/join +one weighted AND +one horizontal reduction +``` + +Then identify which instructions LLVM combined or folded. In particular, verify +that you do not see eight scalar loads or eight scalar branches. + +### Pass the gate + +Find the paired loads and one eight-halfword reduction in your release assembly. +The randomized correctness test is still mandatory; assembly shape cannot prove +the lane-to-bit mapping. + +## Checkpoint 11: make the mask iterator the normal API + +The temporary `find_block_mask()` method was useful for A/B testing. Remove it +and put its state into the existing generic `IterFind`. `find()` initializes the +first block from prefix maxima: + +```rust +let off = self + .block_prefix_max_ends + .partition_point(|max_end| *max_end <= start) + * INDEX_BLOCK_SIZE; + +IterFind { + inner: self, + next_block_start: off, + mask_block_start: 0, + mask: 0, + dense_next: 0, + dense_end: 0, + backend: detect_backend(), + start, + stop, +} +``` + +`seek()` still uses its sorted-query cursor and `max_len`, but it must round down +to the containing block: + +```rust +next_block_start: (*cursor / INDEX_BLOCK_SIZE) * INDEX_BLOCK_SIZE, +``` + +Rounding up would miss earlier lanes in the cursor's block that can still overlap +the query. The mask will reject lanes that are actually too early. + +### Invariant + +The API has one query algorithm. Backend selection changes how a mixed mask is +calculated, not whether the block index is used. + +### Work it by hand + +If `cursor = 45` and block size is 32, the containing block begins at 32, not 64. +A long interval at lane 35 may overlap even though the cursor has advanced to 45. + +### Pass the gate + +Remove the temporary `got_mask` test path. The ordinary call must now pass the +same oracle: + +```bash +cargo test --test block_index +rg -n 'find_block_mask' src tests +``` + +The final `rg` should return nothing. + +## Checkpoint 12: separate the mask backend and support integer widths + +Move mask generation to private `src/simd.rs`. Keep one semantic entry point: + +```rust +pub(crate) fn overlap_mask( + backend: MaskBackend, + starts: &[I], + stops: &[I], + query_start: I, + query_stop: I, +) -> u32 +where + I: PrimInt + 'static, +``` + +Select only CPU capability: + +```rust +pub(crate) enum MaskBackend { + Scalar, + #[cfg(target_arch = "aarch64")] + Neon, + #[cfg(target_arch = "x86_64")] + Avx2, +} + +pub(crate) fn detect_backend() -> MaskBackend { + #[cfg(target_arch = "aarch64")] + { + return MaskBackend::Neon; + } + + #[cfg(target_arch = "x86_64")] + { + if std::is_x86_feature_detected!("avx2") { + return MaskBackend::Avx2; + } + } + + #[allow(unreachable_code)] + MaskBackend::Scalar +} +``` + +Use exact `TypeId` checks before reinterpreting a generic slice as a primitive +slice. This is why the generic query implementation acquires `I: 'static`: + +```rust +if TypeId::of::() == TypeId::of::() { + let starts = unsafe { + std::slice::from_raw_parts( + starts.as_ptr().cast::(), + starts.len(), + ) + }; + // Convert stops and query bounds under the same exact-type proof. + return unsafe { neon::mask_u32(starts, stops, query_start, query_stop) }; +} +``` + +Do not dispatch from `size_of::()` alone. A custom integer wrapper can have +the same size without having a primitive integer's layout. Unknown `PrimInt` +implementations must use the scalar mask. + +The AArch64 packing map is: + +| Type | Values in paired 128-bit loads | Packing idea | +|---|---:|---| +| `u8` / `i8` | 32 | weighted byte halves | +| `u16` / `i16` | 16 | narrow truth to bytes | +| `u32` / `i32` | 8 | narrow truth to halfwords | +| `u64` / `i64` | 4 | narrow truth to words | + +On x86-64 AVX2: + +| Type | Values per 256-bit vector | Bit extraction | +|---|---:|---| +| `u8` / `i8` | 32 | `_mm256_movemask_epi8` | +| `u16` / `i16` | 16 | byte movemask, then compact duplicate bits | +| `u32` / `i32` | 8 | `_mm256_movemask_ps` | +| `u64` / `i64` | 4 | `_mm256_movemask_pd` | + +AVX2 integer greater-than operations are signed. For unsigned types, XOR both +operands with the sign bit before the signed comparison: + +```text +unsigned order: 0 ............ MAX +XOR sign bit: MIN_SIGNED ... MAX_SIGNED +``` + +The transform preserves order while moving it into the signed domain. + +### Invariant + +Every backend returns the same `u32` lane mask. CPU dispatch is not workload +selection. `usize`/`isize` are mapped only to their exact 64-bit representation +on these 64-bit targets. + +### Work it by hand + +Explain why an AVX2 `u16` comparison produces two identical sign bits in the +byte movemask for each true lane, and why those bits must be compacted to one bit +per interval. + +### Pass the gate + +```bash +cargo test --all-features +cargo test --target x86_64-apple-darwin --all-features --lib +cargo clippy --target wasm32-wasip1 --lib --all-features -- -D warnings +``` + +The Rosetta x86 test may take the scalar backend if AVX2 is not advertised. A +forced-Haswell assembly build is a separate instruction check, not native x86 +performance evidence. + +## Checkpoint 13: make every derived array rebuildable + +At this point the query is fast but historical mutation paths leave new metadata +stale. Treat `intervals` as the source of truth and centralize all sidecar +construction: + +```rust +fn rebuild_derived(&mut self) { + let (starts, stops_by_start): (Vec<_>, Vec<_>) = self + .intervals + .iter() + .map(|interval| (interval.start, interval.stop)) + .unzip(); + + self.starts = starts; + self.stops = stops_by_start.clone(); + self.stops_by_start = stops_by_start; + + self.max_len = self + .intervals + .iter() + .map(|interval| { + interval + .stop + .checked_sub(&interval.start) + .unwrap_or_else(zero::) + }) + .max() + .unwrap_or_else(zero::); + + self.stops.sort(); + + self.block_max_ends.clear(); + self.block_min_ends.clear(); + let block_count = self.intervals.len().div_ceil(INDEX_BLOCK_SIZE); + + for block in self.intervals.chunks(INDEX_BLOCK_SIZE) { + let mut max_end = block[0].stop; + let mut min_end = max_end; + for interval in &block[1..] { + max_end = std::cmp::max(max_end, interval.stop); + min_end = std::cmp::min(min_end, interval.stop); + } + self.block_max_ends.push(max_end); + self.block_min_ends.push(min_end); + } + + self.block_index.clear(); + self.block_index.resize(block_count, block_count); + let mut stack = Vec::::new(); + for block in (0..block_count).rev() { + while stack.last().is_some_and(|next| { + self.block_max_ends[*next] <= self.block_max_ends[block] + }) { + stack.pop(); + } + self.block_index[block] = + stack.last().copied().unwrap_or(block_count); + stack.push(block); + } + + self.block_prefix_max_ends + .clone_from(&self.block_max_ends); + for block in 1..block_count { + self.block_prefix_max_ends[block] = std::cmp::max( + self.block_prefix_max_ends[block - 1], + self.block_prefix_max_ends[block], + ); + } +} +``` + +Preserve the existing `sort_unstable` feature branches around the two sorts in +the actual source. + +Call `rebuild_derived()` from: + +```text +new() +insert(), after inserting into canonical order +merge_overlaps(), after replacing intervals +deserialize(), by constructing through new() +``` + +Do not serialize CPU- or implementation-specific sidecars. Custom serde writes +the original six fields and rebuilds the new private fields from `intervals` on +read. + +The serialization side deliberately names only the original fields: + +```rust +let mut state = serializer.serialize_struct("Lapper", 6)?; +state.serialize_field("intervals", &self.intervals)?; +state.serialize_field("starts", &self.starts)?; +state.serialize_field("stops", &self.stops)?; +state.serialize_field("max_len", &self.max_len)?; +state.serialize_field("cov", &self.cov)?; +state.serialize_field("overlaps_merged", &self.overlaps_merged)?; +state.end() +``` + +Deserialize those six fields into a helper, but trust only the canonical +intervals as structural input: + +```rust +let serialized = SerializedLapper::::deserialize(deserializer)?; +let mut lapper = Self::new(serialized.intervals); +lapper.cov = serialized.cov; +lapper.overlaps_merged = serialized.overlaps_merged; +Ok(lapper) +``` + +The serialized `starts`, `stops`, and `max_len` fields remain in the helper for +wire compatibility, but `new()` recalculates them and every new sidecar. This +also prevents stale or inconsistent serialized metadata from entering the unsafe +query proof. + +### Invariant + +After every supported mutation: + +```text +starts.len() == stops_by_start.len() +metadata.len() == ceil(starts.len() / 32) +starts[i] and stops_by_start[i] describe canonical interval i +every block link is a valid later block or the block-count sentinel +``` + +### Work it by hand + +Insert a new interval at canonical index 17. Every later lane position changes, +and potentially every block minimum, maximum, link, and prefix maximum changes. +Updating only `starts` and `stops` cannot be correct. A full rebuild is simple and +appropriate because `insert()` is already documented as inefficient. + +### Pass the gate + +Add tests that cross several block boundaries: + +```bash +cargo test --all-features insert_rebuilds_every_query_index +cargo test --all-features merge_rebuilds_every_query_index +cargo test --all-features serde_keeps_the_v1_six_field_representation +``` + +## Checkpoint 14: remove unsigned assumptions + +Supporting signed primitive coordinates requires more than writing signed SIMD +comparisons. + +Remove the `Unsigned` bound, then repair the scalar assumptions below. + +First, a valid signed interval can have a mathematical length larger than +`I::max_value()`. Record that case while rebuilding instead of treating the +failed subtraction as a zero-length interval: + +```rust +self.max_len = zero::(); +self.max_len_overflowed = false; +for interval in &self.intervals { + match interval.stop.checked_sub(&interval.start) { + Some(length) => self.max_len = std::cmp::max(self.max_len, length), + None if interval.stop >= interval.start => self.max_len_overflowed = true, + None => {} + } +} +``` + +For ordinary `seek()` calls, saturation belongs at the type minimum, not zero. +When `max_len_overflowed` is set, use the exact prefix index because no value of +`I` can provide a conservative length bound: + +```rust +if self.max_len_overflowed { + *cursor = self + .block_prefix_max_ends + .partition_point(|max_end| *max_end <= start) + * INDEX_BLOCK_SIZE; +} else { + let earliest_start = start + .checked_sub(&self.max_len) + .unwrap_or_else(I::min_value); + // Retain the normal monotonic-cursor search. +} +``` + +For `count()`, avoid `start + 1`, which can overflow at the coordinate maximum: + +```rust +let ends_before_or_at_start = self + .stops + .partition_point(|interval_stop| *interval_stop <= start); +let starts_before_stop = Self::bsearch_seq(stop, &self.starts); + +let count = starts_before_stop - ends_before_or_at_start; +``` + +The source expresses the same arithmetic through total length and excluded +suffix count. The identity is: + +```text +overlaps = intervals with start < query_stop + - intervals with stop <= query_start +``` + +For `depth()`, zero cannot be an uninitialized sentinel. Add an explicit +`initialized: bool` so a coverage interval crossing zero starts at its real +negative coordinate. After advancing to the merged interval's stop, break +before constructing another one-unit query; that stop may be `I::max_value()`. + +### Invariant + +Zero is an ordinary coordinate. Overflow behavior and iterator initialization +must not give it special semantic meaning. + +### Work it by hand + +For ends `[-10, -2, 4, 9]` and query start `-2`, exactly the first two ends are +`<= -2` and cannot overlap. No `+1` conversion is required. + +### Pass the gate + +```bash +cargo test --all-features every_primitive_integer_type_and_block_tail_matches_forward_brute_force +cargo test --all-features signed_seek_saturates_at_the_coordinate_minimum +cargo test --all-features signed_seek_keeps_intervals_whose_length_exceeds_the_coordinate_type +cargo test --all-features signed_depth_crosses_zero_once +cargo test --all-features signed_depth_stops_at_the_coordinate_maximum +``` + +## Checkpoint 15: remove only proven bounds checks + +Finish the safe implementation first and inspect its optimized assembly. The +SIMD helpers already use pointer loads; they do not retain a Rust slice check per +vector lane. The remaining useful reduction was private sidecar indexing in the +block loop. + +The proof boundary is: + +```rust +let block_start = self.next_block_start; +if block_start >= self.inner.starts.len() { + return None; +} +``` + +After construction or supported mutation, private invariants prove: + +```text +stops_by_start.len() == starts.len() +each metadata vector has ceil(starts.len() / 32) entries +block = block_start / 32 names the metadata for this boundary +block_end is clamped to starts.len() +``` + +Only then replace repeated private accesses: + +```rust +let first_start = unsafe { + *self.inner.starts.get_unchecked(block_start) +}; +let max_end = unsafe { + *self.inner.block_max_ends.get_unchecked(block) +}; +let starts = unsafe { + self.inner.starts.get_unchecked(block_start..block_end) +}; +``` + +Keep returned interval indexing safe: + +```rust +return Some(&self.inner.intervals[index]); +``` + +The public `intervals` vector could be edited directly, which has always made +derived metadata stale. Bounding unsafe traversal by the private `starts.len()` +prevents direct growth from extending unchecked sidecars; safe result indexing +turns direct shrinkage into a possible panic rather than unchecked result access. + +### Invariant + +Every unsafe access must cite a private length relation established by +`rebuild_derived()`. "The benchmark did not crash" is not a safety proof. + +### Work it by hand + +If `starts.len() = 65`, metadata length is `ceil(65 / 32) = 3`. Valid block +starts are 0, 32, and 64; their block numbers are 0, 1, and 2. The final +`block_end` is `min(96, 65) = 65`, so the final SIMD slice has one lane. + +### Pass the gate + +```bash +cargo test --all-features +cargo clippy --all-targets --all-features -- -D warnings +git diff contender/simd-mask-narrow8-unchecked-index -- src/lib.rs +``` + +The measured specialized AArch64 function fell from 442 to 382 static +instructions and from nine bounds-panic edges to two. Making result yields +unchecked reached 363/zero but improved only about 0.3-0.6%, so that broader +unsafe surface was rejected. + +## Final trace: run the complete search in your head + +Use the toy blocks from checkpoint 3 and query `[23, 28)`. + +1. Prefix maxima `[30, 30, 60]` return block 0 as the first possible block. +2. Block 0 has `min=4`, `max=30`, so it is mixed. +3. Its mask is `0b0010`; the iterator saves it. +4. `next()` drains lane 1 and returns canonical interval `[2,30)`. +5. The mask is empty, so the iterator enters block 1. +6. Block 1 has `max=21 <= 23`; its next-greater link jumps to block 2. +7. Block 2 has `min=25 > 23`; every end passes. +8. The start prefix below query stop 28 contains lanes 0, 1, and 2. +9. Three later `next()` calls return `[22,60)`, `[24,25)`, and `[27,35)`. +10. Lane 3 starts at 40, so the globally sorted starts prove the query is done. + +The returned canonical indices are `[1, 8, 9, 10]`. That single trace exercises +the mixed mask, saved iterator state, forward block skip, dense prefix, and global +right-hand stopping proof. + +## The query algorithms side by side + +The sidecars support three distinct public query shapes: + +| API | Entry strategy | Result strategy | +|---|---|---| +| `find(start, stop)` | Binary-search `block_prefix_max_ends` | Lazily enumerate borrowed intervals through block routes | +| `seek(start, stop, cursor)` | Reuse a monotonic-query cursor bounded by `max_len`, or use the exact block-prefix search if a signed length overflowed, then round down to a block | Use the same lazy block iterator as `find()` | +| `count(start, stop)` | Binary-search global `starts` and global `stops` | Subtract the two excluded endpoint populations without enumeration | + +`find()` answers arbitrary queries independently. `seek()` is the same overlap +problem with an extra promise from the caller: query starts arrive in sorted +order, so a cursor can avoid repeating the entry binary search. It still returns +the same references in the same order. + +`count()` uses a different identity because it does not need interval values: + +```text +starts_before_query_stop = number of interval starts < query_stop +ends_at_or_before_query_start = number of interval stops <= query_start + +overlap_count = starts_before_query_stop - ends_at_or_before_query_start +``` + +For a valid half-open query, every interval in the second population is also in +the first, so the subtraction is nonnegative. `union_and_intersect()` and +`depth()` compose the enumerating APIs; they do not introduce another overlap +index. + +## The final data structure at a glance + +| State | Order | Purpose | +|---|---|---| +| `intervals` | `(start, stop)` | Canonical storage and returned references | +| `starts` | Same as `intervals` | Prefix searches, stopping, SIMD starts | +| `stops` | Globally sorted | `count()` binary search | +| `stops_by_start` | Same as `intervals` | SIMD ends paired with lanes | +| `block_max_ends` | One per 32 | Prove an all-miss block | +| `block_min_ends` | One per 32 | Prove every end passes | +| `block_index` | One per 32 | First later block with greater max end | +| `block_prefix_max_ends` | Running max | Binary-search first possible block | +| `max_len` | Scalar | Conservative cursor movement in `seek()` | +| `max_len_overflowed` | Scalar flag | Select exact prefix entry when no `I`-sized length bound is conservative | + +Complexity is: + +```text +build: O(n log n) sorting + O(n) sidecar construction +space: O(n) coordinate sidecars + O(n / 32) block metadata +count: O(log n) +find: O(log blocks) entry search + candidate block work + output count +seek: cursor reuse for sorted queries + the same block iterator +``` + +## Why this remains Lapper-shaped + +The final index does not decompose intervals into containment-derived lists like +AIList, build an NCList containment hierarchy, follow SuperIntervals' +previous-greater links backward, or reorder nodes into a tree. It retains one +canonical start-sorted vector, fixed positional block summaries, forward links, +and a query-local mask drained low bit first. + +Maximum-end summaries, monotonic stacks, binary search, and SIMD comparison are +general techniques. The structure's identity comes from how they are composed +and from the preserved forward iterator contract. + +## Final comprehension checkpoint + +Do not look at the answers above while explaining these aloud: + +1. Why is `end == query_start` a miss? +2. What does a next-greater link prove about every skipped entry? +3. Why did per-interval links become 32-entry block links? +4. Why is a prefix maximum needed in addition to each block maximum? +5. Why does the dense route still need an active start prefix? +6. Why does the mixed route not need a separate active prefix? +7. Why must `stops` and `stops_by_start` both exist? +8. Why does clearing the lowest set bit preserve forward order? +9. Why is the weighted horizontal sum a mask rather than a count? +10. Why is narrowing all-ones truth from `u32` to `u16` lossless? +11. Why must `seek()` round its cursor down to a block boundary? +12. What makes CPU dispatch different from a workload mode switch? +13. Which invariant makes each `get_unchecked` access valid? +14. Why are returned intervals still safely indexed? +15. Which arrays must change after one insertion, and why is rebuilding simpler? + +Once you can answer all fifteen comprehension questions and reproduce the final +toy trace, read +[`PORTABLE_SIMD_INDEX.md`](PORTABLE_SIMD_INDEX.md) for the engineering summary +and compare your complete practice worktree with +`worked/portable-simd-index`. diff --git a/plans/simd/CONTENDER.md b/plans/simd/CONTENDER.md new file mode 100644 index 0000000..5185cd5 --- /dev/null +++ b/plans/simd/CONTENDER.md @@ -0,0 +1,67 @@ +# Worked portable SIMD block index + +This branch integrates the best forward-order contender into rust-lapper's +normal API. `find()` and `seek()` retain their existing signatures and lazy, +ascending-start output. There is no workload score, data classifier, mode flag, +or alternate public query method. + +The implementation keeps one canonical start-sorted interval vector and adds +fixed 32-entry block facts: + +- a maximum end and a next-greater-block link for exact forward skips; +- a minimum end for an exact all-ends-pass route; +- a prefix maximum for finding the first possible block; +- ends in start order for mixed-block masks. + +Mixed blocks use NEON on AArch64, AVX2 when available on x86-64, and an exact +scalar fallback elsewhere. The SIMD dispatch covers all primitive integer +coordinate types from 8 through 64 bits, including `usize` and `isize`. +Custom `PrimInt` implementations use the scalar mask. + +`insert()`, `merge_overlaps()`, construction, and serde deserialization all run +the same derived-index rebuild. Signed coordinates, minimum-value queries, +`count()` at a type maximum, and depth spans crossing zero have dedicated tests. +Serde keeps the original six-field representation and reconstructs the new +private sidecars when reading it. + +The hot loop uses `get_unchecked` only after a length invariant has proved the +private sidecars cover the candidate block. Returned intervals remain safely +indexed. The assembly audit found no reason to add handwritten assembly: the +AArch64 compiler output already contains paired loads, narrowing, and horizontal +reduction, while forced-Haswell output contains the expected AVX2 compares and +movemask instructions. + +Apple M3 AArch64 five-library total medians: + +| Case | SuperIntervals | Worked Lapper | COITrees | rust-bio IITree | rust-bio AVL | +|---|---:|---:|---:|---:|---:| +| `1-2` | 6.063 ms | **5.729 ms** | 10.415 ms | 13.311 ms | 42.062 ms | +| `7-3` | 68.881 ms | **66.385 ms** | 89.782 ms | 149.975 ms | 344.913 ms | +| `8-7` | **550.372 ms** | 588.766 ms | 834.219 ms | 1259.076 ms | 2140.363 ms | + +Alternating direct trials confirmed the shape: Lapper won total time on 15/15 +`1-2` pairs, was 1.76% faster on `7-3`, and trailed by about 9.5% total on +dense `8-7`. + +AMD Ryzen 9 3950X AVX2 five-library total medians: + +| Case | SuperIntervals | Worked Lapper | COITrees | rust-bio IITree | rust-bio AVL | +|---|---:|---:|---:|---:|---:| +| `1-2` | 11.968 ms | **8.600 ms** | 13.886 ms | 18.150 ms | 68.153 ms | +| `7-3` | **90.807 ms** | 97.084 ms | 136.789 ms | 221.536 ms | 539.133 ms | +| `8-7` | **685.273 ms** | 779.564 ms | 1101.420 ms | 1783.182 ms | 3302.566 ms | + +Runtime dispatch selected AVX2 and every implementation returned the expected +overlap counts. Alternating SuperIntervals comparisons put worked totals 19.70% +ahead on `1-2`, 5.38% behind on `7-3`, and 12.67% behind on `8-7`. Against +rust-lapper 1.3.0, worked totals improved by 37.30%, 98.89%, and 34.64%; +pathological `7-3` was 90.2 times faster. Raw samples and complete host/build +metadata are retained in the +[native AVX2 bakeoff record](https://github.com/sstadick/lapper_bakeoff/tree/main/results/avx2-2026-07-28). + +See [`PORTABLE_SIMD_INDEX.md`](PORTABLE_SIMD_INDEX.md) for the full design, +safety argument, compatibility notes, tests, measurements, and identity audit. +Use [`AARCH64_U32_HAND_TYPE.md`](AARCH64_U32_HAND_TYPE.md) on this branch as the +single hand-typing course from original scalar Lapper through the full portable +implementation. At its ARM checkpoints, the intentionally narrow answer code is +the diff from `contender/block-32` to `tutorial/aarch64-u32-hand-typed`. diff --git a/plans/simd/LAPPER_BENCHMARK_HANDOFF.md b/plans/simd/LAPPER_BENCHMARK_HANDOFF.md new file mode 100644 index 0000000..933c937 --- /dev/null +++ b/plans/simd/LAPPER_BENCHMARK_HANDOFF.md @@ -0,0 +1,138 @@ +# rust-lapper 2 benchmark handoff + +> This is a reproducibility record for the exact pre-release benchmark commit +> named below, not current installation or API guidance. The release-facing +> summary lives in [`README.md`](../../README.md). + +## Objective + +Reproduce the rust-lapper 2 release performance check without relying on the +existing `lapper_bakeoff` checkout or its dirty working tree. The required +implementation revisions, workload semantics, data source, expected counts, +and reference results are below. + +The release question is: does the portable block/SIMD implementation avoid a +regression against rust-lapper 1.3.0 on the three retained article workloads, +especially the pathological `7-3` case? + +## Exact revisions and build + +- Candidate query implementation: `sstadick/rust-lapper`, branch + `worked/portable-simd-index`, commit + `0eeeafdd0773e10fc1ef7f007e69f471917562d7`. +- Baseline: rust-lapper 1.3.0, commit + `c545d40` (also tag `v1.3.0`). +- Compiler used for the retained x86 record: Rust 1.95.0. +- Release flags: `RUSTFLAGS="-C target-cpu=native"`. +- Release profile: thin LTO and one code-generation unit. + +The later local 2.0.0-beta.1 release-preparation edits change metadata, +documentation, tests, and development dependencies, not the measured query +implementation, so `0eeeafd` is the correct reproducible performance revision. + +Prefer a physical x86-64 host with AVX2. Record `lscpu`, OS/kernel, rustc +version, and the full build flags. Before timing the candidate, verify dispatch: + +```sh +cargo +1.95.0 test --all-features --locked \ + simd::dispatch_tests::selected_backend_matches_the_host -- --exact --nocapture +cargo +1.95.0 test --all-features --locked \ + simd::avx2::tests::primitive_masks_match_scalar_when_avx2_is_available \ + -- --exact --nocapture +``` + +The first test must report/select AVX2 on an AVX2 benchmark host, and the direct +primitive-mask test must pass. + +## Data + +The source archive is approximately 401 MB: + +```text +https://drive.usercontent.google.com/download?id=1lctmude31mSAh9fWjI60K1bDrbeDPGfm&export=download&confirm=t +``` + +It contains Parquet directories under `databio/`. Convert these three columns +to headerless, tab-separated BED using Polars 1.32.3: + +```text +contig, pos_start, pos_end +``` + +Required directories and integrity counts: + +| Article ID | Source directory | BED name | Intervals | +|---|---|---|---:| +| 1 | `fBrain-DS14718` | `fBrain.bed` | 198,621 | +| 2 | `exons` | `exons.bed` | 438,694 | +| 3 | `chainOrnAna1` | `chainOrnAna1.bed` | 1,956,864 | +| 7 | `ex-anno` | `ex-anno.bed` | 1,194,285 | +| 8 | `ex-rna` | `ex-rna.bed` | 9,944,559 | + +All coordinates are zero-based, half-open `[start, stop)` ranges. + +## Harness contract + +A small temporary Rust binary is sufficient; the five-library bakeoff is not +required for the v1-versus-v2 release gate. + +1. Read both BED files before starting any timer. +2. Group ranges by chromosome in a `BTreeMap`. +3. For each reported repeat, build a fresh + `BTreeMap>`, one `Lapper` per chromosome, and time + construction separately. +4. Time queries by iterating every query range and evaluating + `index.find(start, stop).count()`. Sum and black-box the result. +5. Do **not** substitute `Lapper::count`; the workload measures result + enumeration through `find`. +6. Emit raw TSV fields: + `revision, repeat, database_intervals, query_intervals, overlaps, build_seconds, query_seconds`. +7. Run at least three reported repeats per revision and case. Retain every + sample. Report build and query medians independently, then add them for the + total. + +Use each case in this physical database/query orientation: + +| Case | Indexed database | Query input | Expected overlaps | +|---|---|---|---:| +| `1-2` | `exons.bed` | `fBrain.bed` | 54,246 | +| `7-3` | `chainOrnAna1.bed` | `ex-anno.bed` | 4,408,383 | +| `8-7` | `ex-anno.bed` | `ex-rna.bed` | 307,184,634 | + +The orientation is intentional. Case `7-3` reproduces rust-lapper 1.3.0's +global-maximum-interval-length worst case. Abort if either revision returns a +different overlap count. + +For clean dependency switching, create two temporary worktrees at the exact +commits and build the same harness once against each path. Do not benchmark a +debug build. Alternate revision execution order where practical, and keep the +machine otherwise idle. + +## Reference result + +The retained native AVX2 run used an AMD Ryzen 9 3950X, Ubuntu 24.04.3, Linux +6.8.0-100-generic, rustc 1.95.0, thin LTO, one codegen unit, and +`-C target-cpu=native`. Independent three-sample medians were: + +| Case | v1.3.0 total | Candidate total | Candidate change | +|---|---:|---:|---:| +| `1-2` | 13.715 ms | 8.600 ms | -37.30% | +| `7-3` | 8,756.427 ms | 97.084 ms | -98.89% | +| `8-7` | 1,192.626 ms | 779.564 ms | -34.64% | + +The `7-3` total improved by 90.2 times. Absolute timings will vary by host; the +release gate is identical counts and no unexplained regression versus v1.3.0. + +## Deliverables + +Retain: + +- raw TSV samples and stderr/load logs; +- CPU, OS/kernel, compiler, commit IDs, and exact flags; +- interval and overlap-count checks; +- the untrimmed median calculation; and +- a short Markdown summary stating whether any case regressed. + +The repository's existing `cargo bench --bench lapper_benchmark` is useful as a +synthetic smoke benchmark, but it does not replace these three release +workloads. diff --git a/plans/simd/PORTABLE_SIMD_INDEX.md b/plans/simd/PORTABLE_SIMD_INDEX.md new file mode 100644 index 0000000..83f4893 --- /dev/null +++ b/plans/simd/PORTABLE_SIMD_INDEX.md @@ -0,0 +1,450 @@ +# Portable forward-order SIMD index + +> This document is the retained design and validation record. See +> [`README.md`](../../README.md) and the generated crate documentation for the +> concise user-facing contract. + +## Goal + +This branch takes the fastest retained rust-lapper contender and makes it a +complete implementation rather than a benchmark-only `u32` experiment. + +The constraints are deliberate: + +- keep `Lapper::find()` and `Lapper::seek()` as the query API; +- yield the same interval references in ascending start order; +- use the same exact algorithm for every workload; +- rebuild every derived array after supported mutation; +- support signed and unsigned primitive coordinates from 8 through 64 bits; +- accelerate AArch64 and x86-64 without making other targets incorrect; and +- use `unsafe` only where a documented invariant removes measured hot-loop work. + +CPU feature dispatch is not a query mode. A query always follows the same block +algorithm. Only the instruction sequence used to calculate an exact mixed-block +mask changes with the CPU and coordinate type. + +## Public behavior + +The production path is the ordinary one: + +```rust +for interval in lapper.find(query_start, query_stop) { + // Same borrowed Interval, same ascending-start order. +} + +for interval in lapper.seek(query_start, query_stop, &mut cursor) { + // Same cursor-facing API. +} +``` + +The experimental `find_block_mask()` method is gone. No mode setting or score is +exposed. `Lapper`, `Interval`, `IterFind`, mutation methods, iterators, and serde +field names remain available as before. + +The generic implementation now requires `I: 'static` so it can use `TypeId` to +select only exact built-in integer representations before any SIMD pointer cast. +This is invisible for all primitive coordinate types. It is a theoretical source +compatibility constraint for an exotic lifetime-carrying custom `PrimInt`. + +## Derived layout + +The canonical `intervals` vector remains sorted by `(start, stop)`. The index does +not reorder that vector. Construction derives: + +| Array | Length | Purpose | +|---|---:|---| +| `starts` | `n` | Start binary searches and SIMD input | +| `stops` | `n` | Globally sorted ends for `count()` | +| `stops_by_start` | `n` | End paired with each start for a mask lane | +| `block_max_ends` | `ceil(n / 32)` | Prove a whole block misses | +| `block_min_ends` | `ceil(n / 32)` | Prove every active lane's end passes | +| `block_index` | `ceil(n / 32)` | Next block to the right with a greater maximum end | +| `block_prefix_max_ends` | `ceil(n / 32)` | Binary-search the first possible block | + +The block size remains 32. It maps one block result to one `u32` mask, amortizes +metadata, and was the best compromise across the measured sparse and dense cases. +Making SIMD lanes wider changes how quickly a mixed block is classified; it does +not remove the metadata, output, and partial-block tradeoffs that made 64-entry +blocks regress in an earlier trial. + +`block_index` is built with a reverse monotonic stack. If block `b` has +`max_end <= query.start`, its link points to the first later block with a strictly +greater maximum. Every skipped intermediate block has a maximum no greater than +block `b`, so all of them are also exact misses. + +## One always-on traversal + +`find()` binary-searches `block_prefix_max_ends` to skip the leading prefix whose +intervals all end at or before `query.start`. `seek()` retains its cursor behavior, +then rounds the candidate position down to its containing 32-entry block. The +normal block loop has three exact routes: + +```text +block.max_end <= query.start + -> all 32 miss; follow the next-greater block link + +block.min_end > query.start + -> every end passes; return starts < query.stop as a dense prefix + +otherwise + -> starts and ends are mixed; calculate and save an exact overlap mask +``` + +Before any route, the first block start is compared with `query.stop`. Starts are +globally sorted, so failure ends the iterator. In the dense route a partition +point finds the active start prefix. In the mixed route the mask checks both +half-open overlap predicates: + +```text +interval.stop > query.start +interval.start < query.stop +``` + +This removes the old active-prefix special case from the mixed path. SIMD checks +the candidate slice that exists, including start failures; those lanes simply +become zero. + +## Why order is unchanged + +The iterator stores one mixed-block mask. Each `next()` call removes its least +significant set bit: + +```rust +let lane = self.mask.trailing_zeros() as usize; +self.mask &= self.mask - 1; +return Some(&self.inner.intervals[self.mask_block_start + lane]); +``` + +The lowest bit is the lowest position in the start-sorted block. Blocks are only +visited to the right, dense prefixes are walked from their first item, and skip +links only point right. Therefore results remain references into the canonical +vector in ascending start order. SIMD changes how a block answer is calculated, +not the observable traversal identity. + +## SIMD backends and integer widths + +`src/simd.rs` owns one private operation: + +```rust +overlap_mask(starts, stops, query_start, query_stop) -> u32 +``` + +It always returns the same lane-to-bit mapping. Exact `TypeId` checks allow SIMD +only for these built-ins: + +```text +u8 i8 u16 i16 u32 i32 u64 i64 usize isize +``` + +Any other `PrimInt` takes the scalar implementation. Partial final vectors also +use a small scalar tail. + +### AArch64 NEON + +NEON is part of baseline AArch64. Its registers are 128 bits, so one register has +16 byte lanes, 8 halfword lanes, 4 word lanes, or 2 doubleword lanes. Paired +loads process two adjacent registers when possible. + +| Coordinate | Values per paired load | Truth-to-bit packing | +|---|---:|---| +| `u8` / `i8` | 32 | two 16-lane weighted byte reductions | +| `u16` / `i16` | 16 | narrow to bytes, then weighted reduction | +| `u32` / `i32` | 8 | narrow to halfwords, then one `addv.8h` | +| `u64` / `i64` | 4 | narrow to words, then one weighted reduction | + +The central `u32` trick is: + +```rust +let eight_truths = vcombine_u16( + vmovn_u32(overlaps_low), + vmovn_u32(overlaps_high), +); +let weights = vld1q_u16([1, 2, 4, 8, 16, 32, 64, 128].as_ptr()); +let bits = vaddvq_u16(vandq_u16(eight_truths, weights)); +``` + +A comparison produces all ones for true and zero for false. Narrowing preserves +that truth. AND keeps a lane's power-of-two weight only when it is true; the +horizontal sum is therefore the bit mask, not a match count. + +### x86-64 AVX2 + +Every iterator performs the standard runtime AVX2 feature check once. AVX2 uses +256-bit registers and native movemask instructions: + +| Coordinate | Values per vector | Mask extraction | +|---|---:|---| +| `u8` / `i8` | 32 | `vpmovmskb` | +| `u16` / `i16` | 16 | `vpmovmskb`, then compact duplicate byte bits | +| `u32` / `i32` | 8 | `vmovmskps` on comparison bits | +| `u64` / `i64` | 4 | `vmovmskpd` on comparison bits | + +AVX2 integer greater-than comparisons are signed. Unsigned inputs are XORed with +their type's sign bit before comparing. That order-preserving transform maps the +unsigned domain onto the signed domain without changing the overlap predicate. + +### Other CPUs and custom integers + +The scalar backend calculates the same `u32` mask lane by lane. This keeps the +algorithm and output identical on targets such as `wasm32-wasip1` and on x86-64 +machines without AVX2. It is a correctness fallback, not a workload-dependent +choice. + +SVE/SVE2 and AVX-512 are possible later backends, but are not present here. Their +predicate/mask facilities could classify a 32-entry block in fewer vector groups. +They need native hardware measurements before becoming a recommendation. + +## Mutation, construction, and serde + +All sidecars are derived state. `rebuild_derived()` is the single construction +path and is called by: + +- `Lapper::new()`; +- `insert()` after its sorted insertion; +- `merge_overlaps()` after it replaces the interval vector; and +- serde deserialization through `Lapper::new()`. + +The rebuild uses one bulk `unzip()` for starts and start-order stops, clones the +stops once for global sorting, then builds block summaries. Restoring this bulk +path recovered build time that an earlier incremental prototype had lost. + +With `with_serde`, serialization deliberately writes the original six fields: +`intervals`, `starts`, `stops`, `max_len`, `cov`, and `overlaps_merged`. New +private sidecars and the length-overflow flag are not serialized. Deserialization +treats `intervals` as the source of truth and rebuilds every derived value, then +restores the cached coverage and merge flag. Existing serialized shape therefore +does not acquire architecture- or implementation-specific metadata. + +## Signed-coordinate corrections + +Generalizing beyond unsigned coordinates exposed several assumptions in the old +code: + +- derived state records when a valid interval length cannot fit in `I`; in that + rare case, `seek()` uses the exact block-prefix search instead of an + underestimated `max_len`; +- otherwise, `seek()` uses `checked_sub(max_len)` and falls back to + `I::min_value()`; +- `count()` uses `partition_point(stop <= start)` instead of forming `start + 1`, + which can overflow at the coordinate maximum; and +- `depth()` has an explicit initialization flag instead of using zero as a + sentinel, and stops at a merged endpoint before forming a one-unit query past + `I::max_value()`. + +These are scalar correctness repairs needed by the wider type support; they are +not query heuristics. + +## Bounds-check and unsafe audit + +Release assembly showed that the SIMD helpers already use pointer vector loads; +there were no hidden per-lane Rust slice checks in the vector core. The remaining +hot-loop checks were on private sidecar indexing. + +The accepted branch uses `get_unchecked` for those sidecars after one public +invariant boundary: + +1. `block_start < starts.len()` is checked. +2. Constructors, supported mutation, and deserialization rebuild `starts` and + `stops_by_start` to equal lengths. +3. They also build one entry in every block metadata vector for each 32-entry + start block. +4. `block_end` is clamped to `starts.len()`. +5. A link is either a valid later block or the block-count sentinel. + +The block link is guarded by a debug assertion. Returned `intervals` remain +safely indexed, so unsupported direct edits to the public vector can at worst +produce stale behavior or a panic, not turn a result yield into unchecked memory +access. Directly appending to `intervals` is ignored by the private derived bound; +supported mutation must continue to use `insert()` or `merge_overlaps()`. + +For the specialized AArch64 query, this change reduced the inlined function from +442 to 382 static instructions and bounds-panic edges from nine to two. Making +the final result access unchecked reduced it further to 363/zero but improved +queries by only about 0.3-0.6%, so that broader unsafe surface was rejected. + +Handwritten assembly was also rejected. LLVM already emits the intended AArch64 +`ldp`, `uzp1.8h`, and `addv.8h` sequence. A forced Haswell build contains +`vpmovmskb`, `vmovmskps`, `vmovmskpd`, `vpcmpgt*`, `vpxor`, and `vpand`. Inline +assembly would add review and register-allocation risk without removing a known +instruction in these cores. + +The detailed audit is in the bakeoff's +[`ASSEMBLY_BOUNDS_AUDIT.md`](https://github.com/sstadick/lapper_bakeoff/blob/main/ASSEMBLY_BOUNDS_AUDIT.md). + +## Correctness and portability checks + +On 2026-07-28, the designated native x86 benchmark host, an AMD Ryzen 9 3950X, +passed: + +```text +cargo test --all-features --locked +cargo clippy --all-targets --all-features --locked -- -D warnings +``` + +Runtime dispatch selected AVX2, and direct calls to every signed and unsigned +primitive AVX2 mask matched the scalar result. The final GitHub matrix also +passed native AVX2 execution on an AMD EPYC 7763, native AArch64 NEON execution +on macOS, and the complete x86-64 scalar suite under a QEMU Nehalem CPU model. +It covers Rust 1.59 and current stable, every feature configuration, Windows, +and scalar-only compilation for i686, PowerPC64LE, and Wasm. + +Coverage includes: + +- randomized `find()`, `seek()`, and `count()` against brute force; +- exact reference order for all primitive integer coordinate types, including + the scalar `u128` and `i128` fallbacks; +- negative intervals, a query at the signed minimum, an interval whose length + exceeds the signed coordinate's positive range, and depth ending at the + signed maximum; +- insert and merge rebuilding across more than four blocks; +- serde round trips with derived metadata reconstruction; +- signed depth across zero; +- safety when the public interval vector is directly extended; +- exhaustive AVX2 `u16` movemask compaction; and +- forced calls to each AVX2 primitive function when AVX2 is available. + +Rosetta provided the initial x86-64 scalar-dispatch check because it did not +advertise AVX2. Native AMD CI and the Ryzen host now supersede that limitation +for AVX2 correctness. The Ryzen host also completed the three retained +performance cases; those results follow the original AArch64 record below. + +## Performance record + +All times below are medians in milliseconds. Totals combine independently +reported build and query medians. + +### Apple M3 AArch64 + +The five-library worked run used the normal `find()` API and native Apple M3 +compilation. + +| Case | SuperIntervals | Worked Lapper | COITrees | rust-bio IITree | rust-bio AVL | +|---|---:|---:|---:|---:|---:| +| `1-2` | 6.063 | **5.729** | 10.415 | 13.311 | 42.062 | +| `7-3` | 68.881 | **66.385** | 89.782 | 149.975 | 344.913 | +| `8-7` | **550.372** | 588.766 | 834.219 | 1259.076 | 2140.363 | + +Direct alternating Lapper/SuperIntervals process pairs give the more reliable +close-call interpretation: + +| Case | Worked Lapper versus SuperIntervals total | Pair wins | +|---|---:|---:| +| `1-2` | -7.71% paired median | 15 / 15 | +| `7-3` | -1.76% paired median | 12 / 15 | +| `8-7` | +9.50% paired median | 0 / 10 | + +On `7-3`, Lapper's query remained about 13% slower but its cheaper build erased +that difference for one build plus one query batch. On dense `8-7`, output work +dominates and SuperIntervals retains a clear advantage. + +The worked generic integration stayed within about 1.3% of the specialized +unchecked `u32` branch in paired query trials: + +| Case | Worked query change versus specialized | Worked wins | +|---|---:|---:| +| `1-2` | -2.14% | 11 / 15 | +| `7-3` | +1.53% | 1 / 15 | +| `8-7` | +1.31% | 2 / 10 | + +For orientation, the original v1.3.0 query medians were 5.648, 5049.867, and +884.546 ms on `1-2`, `7-3`, and `8-7`. Those original-versus-worked values came +from separate controlled runs and should not be read as paired microbenchmark +percentages. Raw files and methods live in the bakeoff's +[`results/forward-simd-2026-07-18`](https://github.com/sstadick/lapper_bakeoff/tree/main/results/forward-simd-2026-07-18) +directory. + +### AMD Ryzen 9 3950X AVX2 + +The native x86-64 run used rustc 1.95.0 and +`RUSTFLAGS="-C target-cpu=native"`. Runtime dispatch selected AVX2, and the +direct primitive mask tests matched the scalar implementation. The worked +binary and all four pinned competitors returned identical overlap counts. + +| Case | SuperIntervals | Worked Lapper | COITrees | rust-bio IITree | rust-bio AVL | +|---|---:|---:|---:|---:|---:| +| `1-2` | 11.968 | **8.600** | 13.886 | 18.150 | 68.153 | +| `7-3` | **90.807** | 97.084 | 136.789 | 221.536 | 539.133 | +| `8-7` | **685.273** | 779.564 | 1101.420 | 1783.182 | 3302.566 | + +Worked Lapper ranked first on `1-2` and second on the other two cases. Relative +to SuperIntervals, its three-repeat totals were 28.14% faster, 6.91% slower, and +13.76% slower. Alternating paired totals give the more reliable direct changes: +19.70% faster, 5.38% slower, and 12.67% slower. + +The same harness also measured rust-lapper 1.3.0 at `c545d40`: + +| Case | Original v1.3.0 | Worked Lapper | Worked total change | +|---|---:|---:|---:| +| `1-2` | 13.715 | 8.600 | **-37.30%** | +| `7-3` | 8756.427 | 97.084 | **-98.89%** | +| `8-7` | 1192.626 | 779.564 | **-34.64%** | + +The pathological `7-3` total improved by 90.2 times. No article case regressed +against the release baseline. The complete method, compiler correction, +alternating comparisons, build/query medians, logs, and raw TSV files are +retained in the bakeoff's +[`results/avx2-2026-07-28`](https://github.com/sstadick/lapper_bakeoff/tree/main/results/avx2-2026-07-28) +directory. + +## Why this is not AIList or another renamed index + +This conclusion is structural, not based on names such as "maximum end." + +| Implementation | Canonical layout | Index shape | Query direction / state | +|---|---|---|---| +| This Lapper | One vector sorted by start | Fixed positional block summaries and forward next-greater links | Forward blocks; query-local low-bit-first mask | +| AIList | Intervals decomposed into multiple containment-derived lists | Per-list running maximum ends | Searches/scans each component backward | +| NCList | Intervals organized by containment | Nested containment hierarchy | Hierarchical sublist traversal | +| SuperIntervals | Separate starts, ends, payloads, and branch arrays | Previous-greater/equal end links | Begins near the query boundary and follows links backward | +| COITrees | Augmented tree nodes in cache-oriented order | Centered interval tree | Tree traversal over augmented nodes | +| rust-bio IITree / AVL | Tree-owned interval nodes | Interval-tree augmentation | Tree traversal | + +This branch does not perform AIList's containment decomposition, create component +lists, or scan those lists backward. It does not build NCList's containment +hierarchy. It does not use SuperIntervals' previous-greater link direction or +separate canonical storage. It does not reorder the intervals into a tree. + +What is shared is general interval-index vocabulary: sorting endpoints, storing +end extrema, skipping regions proved unable to overlap, and using SIMD to evaluate +independent predicates. Those shared techniques do not make the data structures +the same. This comparison is not a patent-clearance or exhaustive prior-art +opinion. + +Primary references used for the structural check: + +- [AIList paper](https://pmc.ncbi.nlm.nih.gov/articles/PMC6901075/) +- [AIList C construction and query](https://github.com/databio/AIList/blob/366fe0a6f5c77cdd165c1452f49adf0ecbbef4e6/src/AIList.c#L130-L317) +- [gtars Rust AIList](https://github.com/databio/gtars/blob/cd75994966f039aa52d9359d64a3a201eee25de7/gtars-overlaprs/src/ailist.rs#L92-L375) +- [SuperIntervals 0.3.6 source](https://github.com/biodatageeks/sequila-native/blob/745d40f77da7ced5d540f9285eb5123ba12682ff/sequila/sequila-core/superintervals/src/superintervals.rs) +- [COITrees 0.4.0 non-SIMD source](https://github.com/dcjones/coitrees/blob/4afeeebe6e5d3229c2c191d5676ed71310d953e0/src/nosimd.rs) + +## Branch map + +| Branch | Commit | Purpose | +|---|---|---| +| `baseline/v1.3.0` | `c545d40` | Original rust-lapper release state | +| `contender/per-interval` | `0af2d23` | Reconstructed next-greater index per interval | +| `contender/block-32` | `559c183` | Preserved always-on 32-entry block index | +| `contender/simd-mask-narrow8-ld1x2` | `84872e7` | Safe specialized AArch64 mask | +| `contender/simd-mask-narrow8-unchecked-index` | `b6ba527` | Specialized best with audited unchecked sidecars | +| `tutorial/aarch64-u32-hand-typed` | `dfc5f01` | Narrow answer key for typing the ARM path once | +| `reference/aarch64-u32-hand-typed` | `b6ba527` | Exact full-reference comparison point | +| `worked/portable-simd-index` | this branch | Mutation-safe, multi-type, multi-architecture integration | + +The unified [`AARCH64_U32_HAND_TYPE.md`](AARCH64_U32_HAND_TYPE.md) course starts +at `baseline/v1.3.0` and walks through every index and SIMD checkpoint. The +`tutorial/aarch64-u32-hand-typed` branch remains a deliberately narrow code +answer for the ARM portion: it excludes generic dispatch, AVX2, mutation +rebuilding, and unchecked indexing so those operations can be learned once +without production scaffolding. + +## Post-release opportunities + +- Physical non-AVX2 x86 and Intel-branded measurements may add hardware-diverse + performance data, but are not v2 release requirements. QEMU already executes + the scalar x86-64 path and native AMD hosts execute the AVX2 path. +- Add native big-endian and non-Apple AArch64 CI if those are support targets. +- Measure memory/cache behavior across coordinate and payload widths. +- Treat SVE2 and AVX-512 as separate CPU backends only after native evidence. +- Preserve safe result indexing unless a new measurement justifies expanding the + unsafe proof. diff --git a/plans/simd/productionization.md b/plans/simd/productionization.md new file mode 100644 index 0000000..d0520cb --- /dev/null +++ b/plans/simd/productionization.md @@ -0,0 +1,134 @@ +# Portable SIMD index production plan + +The five implementation checks below are complete. A short beta is the +remaining step before 2.0.0. + +The implementation on `worked/portable-simd-index` is ready for +2.0.0-beta.1. It preserves forward borrowed iteration order, exact overlap +semantics, and the always-on query algorithm. The accepted `I: 'static` bound +is the intentional major-version API change. + +## 1. Validate AVX2 on native x86-64 hardware + +Status: complete. On 2026-07-28, GitHub's native x86-64 runner reported an AMD +EPYC 7763 with AVX2. Runtime dispatch selected AVX2 and every signed and unsigned +AVX2 mask matched the scalar result. The designated native benchmark host is an +AMD Ryzen 9 3950X. Its locked all-feature suite passes, runtime dispatch selects +AVX2, and the direct primitive mask suite matches the scalar implementation. +Rosetta separately exercised x86-64 scalar selection, and the final CI matrix +ran the complete test suite under a QEMU Nehalem CPU model without AVX2. + +The Ryzen host then ran the three retained article cases with native CPU features +against rust-lapper 1.3.0 and the four pinned Rust competitors. Worked Lapper's +total medians were 8.600, 97.084, and 779.564 ms. That is 37.30%, 98.89%, and +34.64% faster than rust-lapper 1.3.0; the pathological `7-3` total improved by +90.2 times. It ranked first on `1-2` and second on `7-3` and `8-7`. Alternating +comparisons with SuperIntervals put the worked total 19.70% ahead, 5.38% behind, +and 12.67% behind. All implementations returned identical overlap counts. CPU, +compiler, flags, method, raw samples, and medians are retained in the +[native AVX2 bakeoff record](https://github.com/sstadick/lapper_bakeoff/tree/main/results/avx2-2026-07-28). + +Completed: + +- Run the complete test suite on native x86-64 with AVX2 available. +- Verify that the AVX2 backend is selected and executed, rather than only + inspecting forced-target assembly. +- Exercise the complete x86-64 scalar suite under a CPU model without AVX2. +- Run the three retained datasets on the Ryzen 9 3950X against rust-lapper + 1.3.0, the worked Lapper, and the pinned Rust competitors. +- Record CPU, compiler, build flags, raw samples, medians, and overlap counts. + +Release decision: physical non-AVX2 and Intel-branded hosts are not additional +gates. QEMU executes the complete scalar x86-64 suite with AVX2 hidden, while +native AMD hosts execute and benchmark the vendor-neutral AVX2 instruction +path. Physical Intel measurements would add vendor-diverse performance data, +not exercise a different implementation. + +Pass condition met: native AVX2 and modeled non-AVX2 paths are correct, and the +native AVX2 record has no regression against rust-lapper 1.3.0. + +## 2. Decide and document the effective MSRV + +Status: complete. The declared MSRV is Rust 1.59.0. The locked all-feature +library builds on 1.59.0, while 1.58.1 fails because its AArch64 `std::arch` +intrinsics are still unstable. Rust-lapper 1.3.0 had no declared MSRV and its +locked all-feature library still builds on Rust 1.56.1. + +- Determine the oldest compiler supported by rust-lapper 1.3.0. +- Identify the oldest compiler accepted by the worked implementation and its + dependencies. +- Either retain the existing effective MSRV or choose a deliberate new one. +- Add the decision to package metadata, CI, and release documentation. + +Pass condition: the declared MSRV builds and tests the supported feature set, +and any increase is intentional and documented. + +## 3. Resolve the `I: 'static` API-bound addition + +Status: complete and accepted. Stable Rust has no specialization mechanism that +can select primitive SIMD kernels while retaining a blanket scalar path for +every custom `PrimInt`. Safe pointer reinterpretation therefore uses exact +`TypeId` checks, which require `I: 'static`. This excludes only custom coordinate +types carrying non-static borrows; it does not require a `Lapper` value to live +for the program lifetime. Primitive types use SIMD, while `u128`, `i128`, and +other unmatched owned `PrimInt` types use the tested scalar fallback. + +- Measure the public API difference from rust-lapper 1.3.0. +- Determine whether safe primitive-type SIMD dispatch can avoid `TypeId` and + the corresponding `I: 'static` bound without restricting custom `PrimInt` + implementations. +- If it cannot be removed without a larger compatibility or performance cost, + document and test the bound as an intentional compatibility decision. + +Pass condition: the bound is either eliminated or explicitly accepted with its +actual user impact documented. + +## 4. Measure cached backend and type selection + +Status: complete and rejected. A construction-time dispatch experiment stored a +typed mask function pointer in each `Lapper`, removing per-iterator backend +selection and per-mixed-block `TypeId` selection. On an Apple M3, paired +alternating measurements showed query regressions of 4.21% on `1-2`, 2.04% on +`7-3`, and 1.49% on `8-7`. The cached variant won only 1/15, 0/15, and 0/10 +query pairs respectively. The indirect call costs more than the current +compiler-folded checks, so the existing dispatch is retained. + +- Compare the current per-iterator backend detection and per-mask type + selection with selection cached at `Lapper` construction. +- Measure construction, query, and total time on all three retained datasets. +- Reject function-pointer or cached-dispatch designs that are slower, more + fragile, or materially more complex without a demonstrated benefit. + +Pass condition: retain the fastest defensible dispatch arrangement based on +paired measurements, without adding a workload heuristic or user-visible mode. + +## 5. Run the final CI target matrix + +Status: complete. On 2026-07-28, the final GitHub matrix passed native AArch64 +NEON on macOS, native x86-64 AVX2 on Linux, and the complete x86-64 scalar test +suite under a QEMU Nehalem CPU model. It passed Rust 1.59 and current stable, +default, `with_serde`, `sort_unstable`, and all-feature configurations, plus +scalar-only i686, PowerPC64LE, and Wasm compilation. + +- Test AArch64 NEON, x86-64 AVX2, and x86-64 scalar execution. +- Compile and test the scalar fallback on other supported targets. +- Cover default features, `with_serde`, and `sort_unstable`. +- Include the declared MSRV and current stable Rust. + +Pass condition: every supported target and feature combination is green, with +native execution for the SIMD backends claimed by the release. + +## 6. Run the beta + +Status: pending. + +- Publish `2.0.0-beta.1` and confirm the crate installs from crates.io. +- Confirm that docs.rs builds the beta documentation. +- Ask downstream users to compile and test against the beta explicitly. +- Leave the beta open for one or two weeks for correctness and compatibility + reports. +- Promote the tested code to `2.0.0` if no blockers appear. Publish + `2.0.0-beta.2` instead if the beta requires code changes. + +Pass condition: the published crate and documentation work, downstream builds +pass, and there are no unresolved correctness or compatibility issues. diff --git a/src/lib.rs b/src/lib.rs index 0fcc7b2..2fc12d4 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,50 +1,59 @@ -//! This module provides a simple data-structure for fast interval searches. +//! This module provides a simple data structure for fast interval searches. //! ## Features -//! - Extremely fast on most genomic datasets. (3-4x faster than other methods) -//! - Extremely fast on in order queries. (10x faster than other methods) -//! - Extremely fast intersections count method based on the +//! - Extremely fast overlap queries on both ordinary genomic datasets and +//! datasets with long intervals that engulf many shorter intervals. +//! - Extremely fast in order queries through the cursor-based `seek()` method. +//! - Extremely fast intersection counts based on the //! [BITS](https://arxiv.org/pdf/1208.3407.pdf) algorithm -//! - Parallel friendly. Queries are on an immutable structure, even for seek -//! - Consumer / Adapter paradigm, Iterators are returned and serve as the main API for interacting -//! with the lapper +//! - NEON acceleration on AArch64, runtime-detected AVX2 on x86-64, and an exact +//! scalar fallback everywhere else. +//! - Parallel friendly. Queries are on an immutable structure, even for `seek()`. +//! - Consumer / Adapter paradigm. Iterators are returned and serve as the main +//! API for interacting with the Lapper. //! //! ## Details: //! //! ```text -//! 0 1 2 3 4 5 6 7 8 9 10 11 -//! (0,10]X X X X X X X X X X -//! (2,5] X X X -//! (3,8] X X X X X -//! (3,8] X X X X X -//! (3,8] X X X X X -//! (3,8] X X X X X -//! (5,9] X X X X -//! (8,11] X X X +//! 0 1 2 3 4 5 6 7 8 9 10 11 +//! [0, 10) X X X X X X X X X X +//! [2, 5) X X X +//! [3, 8) X X X X X +//! [3, 8) X X X X X +//! [3, 8) X X X X X +//! [3, 8) X X X X X +//! [5, 9) X X X X +//! [8, 11) X X X //! -//! Query: (8, 11] -//! Answer: ((0,10], (5,9], (8,11]) +//! Query: [8, 11) +//! Answer: [0, 10), [5, 9), [8, 11) //! ``` //! -//! Most interaction with this crate will be through the [`Lapper`](struct.Lapper.html) struct -//! The main methods are [`find`](struct.Lapper.html#method.find), -//! [`seek`](struct.Lapper.html#method.seek), and [`count`](struct.Lapper.html#method.count) -//! where both `seek` and `count` are special cases allowing for very fast queries in certain scenarios. +//! Most interaction with this crate will be through the [`Lapper`] struct. The +//! main methods are [`Lapper::find`], [`Lapper::seek`], and [`Lapper::count`]. +//! `find()` handles independent queries, `seek()` reuses a caller-owned cursor +//! when query starts arrive in order, and `count()` is used when only the number +//! of overlaps is needed. //! -//! The overlap function for this assumes a zero based genomic coordinate system. So [start, stop) -//! is not inclusive of the stop position for neither the queries, nor the Intervals. +//! Ranges are half-open: `[start, stop)`. Two ranges overlap when +//! `interval.start < query.stop` and `interval.stop > query.start`, so adjacent +//! ranges such as `[0, 10)` and `[10, 20)` do not overlap. This matches the +//! usual zero-based genomic coordinate system. Signed and unsigned primitive +//! coordinates are supported. //! -//! Lapper does not use an interval tree, instead, it operates on the assumtion that most intervals are -//! of similar length; or, more exactly, that the longest interval in the set is not long compred to -//! the average distance between intervals. +//! Lapper does not use an interval tree. It keeps intervals sorted by start and +//! builds a small index over fixed blocks of 32 intervals. A prefix maximum +//! finds the first block that could overlap; each block's minimum and maximum +//! end positions then prove whether the block is a miss or a dense prefix, and +//! a next-greater link skips runs of blocks that cannot overlap. Mixed blocks +//! produce an exact 32-bit overlap mask with NEON, AVX2, or the scalar fallback. +//! Mask bits are drained from low to high, so results remain borrowed and in +//! ascending start order. //! -//! For cases where this holds true (as it often does with genomic data), we can sort by start and -//! use binary search on the starts, accounting for the length of the longest interval. The advantage -//! of this approach is simplicity of implementation and speed. In realistic tests queries returning -//! the overlapping intervals are 1000 times faster than brute force and queries that merely check -//! for the overlaps are > 5000 times faster. -//! -//! When this is not the case, if possible in your scenario, use merge_overlaps first, and then use -//! `find` or `seek`. The `count` method will be fast in all scenarios. +//! The same block algorithm handles ordinary data and the old worst case where +//! one long interval engulfs many shorter intervals. There is no workload mode +//! to configure. `merge_overlaps()` remains useful when callers want merged +//! coverage, while `count()` remains the independent BITS implementation and is +//! fast regardless of interval shape. //! //! # Examples //! @@ -77,13 +86,17 @@ //! ``` use num_traits::{ identities::{one, zero}, - PrimInt, Unsigned, + PrimInt, }; use std::cmp::Ordering::{self}; use std::collections::VecDeque; +mod simd; + +use simd::{detect_backend, overlap_mask, MaskBackend, BLOCK_SIZE as INDEX_BLOCK_SIZE}; + #[cfg(feature = "with_serde")] -use serde::{Deserialize, Serialize}; +use serde::{ser::SerializeStruct, Deserialize, Deserializer, Serialize, Serializer}; /// Represent a range from [start, stop) /// Inclusive start, exclusive of stop @@ -91,7 +104,7 @@ use serde::{Deserialize, Serialize}; #[derive(Eq, Debug, Clone)] pub struct Interval where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { pub start: I, @@ -99,32 +112,107 @@ where pub val: T, } -/// Primary object of the library. The public intervals holds all the intervals and can be used for -/// iterating / pulling values out of the tree. -#[cfg_attr(feature = "with_serde", derive(Serialize, Deserialize))] +/// Primary interval collection and query index. +/// +/// The public interval vector is the canonical storage and can be read or used +/// to mutate payload values. Coordinate or structural changes must use +/// [`Lapper::insert`] or [`Lapper::merge_overlaps`] so the private query index +/// is rebuilt. #[derive(Debug, Clone)] pub struct Lapper where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { - /// List of intervals + /// Intervals in ascending start order. + /// + /// Directly changing coordinates or vector length leaves the private query + /// index stale. Payload-only changes are safe. pub intervals: Vec>, /// Sorted list of start positions, starts: Vec, /// Sorted list of end positions, stops: Vec, + /// End positions in the same order as `intervals`, for block-mask queries. + stops_by_start: Vec, + /// Index of the next block with a strictly larger maximum end, or the + /// number of blocks when no such block exists. + block_index: Vec, + /// Maximum end position in each fixed-size block. + block_max_ends: Vec, + /// Minimum end position in each fixed-size block. + block_min_ends: Vec, + /// Inclusive prefix maximum of `block_max_ends`, used to find the first + /// candidate block for a query. + block_prefix_max_ends: Vec, /// The length of the longest interval max_len: I, + /// Whether a valid interval length exceeds the positive range of `I`. + max_len_overflowed: bool, /// The calculated number of positions covered by the intervals cov: Option, /// Whether or not overlaps have been merged pub overlaps_merged: bool, } +#[cfg(feature = "with_serde")] +impl Serialize for Lapper +where + I: PrimInt + Ord + Clone + Send + Sync + Serialize, + T: Eq + Clone + Send + Sync + Serialize, +{ + fn serialize(&self, serializer: S) -> Result + where + S: Serializer, + { + let mut state = serializer.serialize_struct("Lapper", 6)?; + state.serialize_field("intervals", &self.intervals)?; + state.serialize_field("starts", &self.starts)?; + state.serialize_field("stops", &self.stops)?; + state.serialize_field("max_len", &self.max_len)?; + state.serialize_field("cov", &self.cov)?; + state.serialize_field("overlaps_merged", &self.overlaps_merged)?; + state.end() + } +} + +#[cfg(feature = "with_serde")] +#[derive(Deserialize)] +#[allow(dead_code)] +struct SerializedLapper +where + I: PrimInt + Ord + Clone + Send + Sync, + T: Eq + Clone + Send + Sync, +{ + intervals: Vec>, + starts: Vec, + stops: Vec, + max_len: I, + cov: Option, + overlaps_merged: bool, +} + +#[cfg(feature = "with_serde")] +impl<'de, I, T> Deserialize<'de> for Lapper +where + I: PrimInt + Ord + Clone + Send + Sync + Deserialize<'de> + 'static, + T: Eq + Clone + Send + Sync + Deserialize<'de>, +{ + fn deserialize(deserializer: D) -> Result + where + D: Deserializer<'de>, + { + let serialized = SerializedLapper::::deserialize(deserializer)?; + let mut lapper = Self::new(serialized.intervals); + lapper.cov = serialized.cov; + lapper.overlaps_merged = serialized.overlaps_merged; + Ok(lapper) + } +} + impl Interval where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { /// Compute the intsect between two intervals @@ -144,7 +232,7 @@ where impl Ord for Interval where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { #[inline] @@ -159,7 +247,7 @@ where impl PartialOrd for Interval where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { #[inline] @@ -170,7 +258,7 @@ where impl PartialEq for Interval where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, T: Eq + Clone + Send + Sync, { #[inline] @@ -181,7 +269,7 @@ where impl Lapper where - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync + 'static, T: Eq + Clone + Send + Sync, { /// Create a new instance of Lapper by passing in a vector of Intervals. This vector will @@ -200,38 +288,87 @@ where #[cfg(not(feature = "sort_unstable"))] intervals.sort(); - let (mut starts, mut stops): (Vec<_>, Vec<_>) = - intervals.iter().map(|x| (x.start, x.stop)).unzip(); + let mut lapper = Lapper { + intervals, + starts: Vec::new(), + stops: Vec::new(), + stops_by_start: Vec::new(), + block_index: Vec::new(), + block_max_ends: Vec::new(), + block_min_ends: Vec::new(), + block_prefix_max_ends: Vec::new(), + max_len: zero::(), + max_len_overflowed: false, + cov: None, + overlaps_merged: false, + }; + lapper.rebuild_derived(); + lapper + } - #[cfg(feature = "sort_unstable")] - { - starts.sort_unstable(); - stops.sort_unstable(); + fn rebuild_derived(&mut self) { + let (starts, stops_by_start): (Vec<_>, Vec<_>) = self + .intervals + .iter() + .map(|interval| (interval.start, interval.stop)) + .unzip(); + self.starts = starts; + self.stops = stops_by_start.clone(); + self.stops_by_start = stops_by_start; + + self.max_len = zero::(); + self.max_len_overflowed = false; + for interval in &self.intervals { + match interval.stop.checked_sub(&interval.start) { + Some(length) => self.max_len = std::cmp::max(self.max_len, length), + None if interval.stop >= interval.start => self.max_len_overflowed = true, + None => {} + } } + #[cfg(feature = "sort_unstable")] + self.stops.sort_unstable(); + #[cfg(not(feature = "sort_unstable"))] - { - starts.sort(); - stops.sort(); + self.stops.sort(); + + self.block_max_ends.clear(); + self.block_min_ends.clear(); + let interval_count = self.intervals.len(); + let block_count = + interval_count / INDEX_BLOCK_SIZE + usize::from(interval_count % INDEX_BLOCK_SIZE != 0); + self.block_max_ends.reserve(block_count); + self.block_min_ends.reserve(block_count); + for block in self.intervals.chunks(INDEX_BLOCK_SIZE) { + let mut max_end = block[0].stop; + let mut min_end = max_end; + for interval in &block[1..] { + max_end = std::cmp::max(max_end, interval.stop); + min_end = std::cmp::min(min_end, interval.stop); + } + self.block_max_ends.push(max_end); + self.block_min_ends.push(min_end); } - let mut max_len = zero::(); - for interval in intervals.iter() { - let i_len = interval - .stop - .checked_sub(&interval.start) - .unwrap_or_else(zero::); - if i_len > max_len { - max_len = i_len; + self.block_index.clear(); + self.block_index.resize(block_count, block_count); + let mut stack = Vec::::new(); + for block in (0..block_count).rev() { + while stack.last().map_or(false, |&next| { + self.block_max_ends[next] <= self.block_max_ends[block] + }) { + stack.pop(); } + self.block_index[block] = stack.last().copied().unwrap_or(block_count); + stack.push(block); } - Lapper { - intervals, - starts, - stops, - max_len, - cov: None, - overlaps_merged: false, + + self.block_prefix_max_ends.clone_from(&self.block_max_ends); + for block in 1..block_count { + self.block_prefix_max_ends[block] = std::cmp::max( + self.block_prefix_max_ends[block - 1], + self.block_prefix_max_ends[block], + ); } } @@ -258,16 +395,9 @@ where /// /// ``` pub fn insert(&mut self, elem: Interval) { - let starts_insert_index = Self::bsearch_seq(elem.start, &self.starts); - let stops_insert_index = Self::bsearch_seq(elem.stop, &self.stops); let intervals_insert_index = Self::bsearch_seq_ref(&elem, &self.intervals); - let i_len = elem.stop.checked_sub(&elem.start).unwrap_or_else(zero::); - if i_len > self.max_len { - self.max_len = i_len; - } - self.starts.insert(starts_insert_index, elem.start); - self.stops.insert(stops_insert_index, elem.stop); self.intervals.insert(intervals_insert_index, elem); + self.rebuild_derived(); self.cov = None; self.overlaps_merged = false; } @@ -315,7 +445,7 @@ where } } - /// Get the number fo positions covered by the intervals in Lapper and store it. If you are + /// Get the number of positions covered by the intervals in Lapper and store it. If you are /// going to be using the coverage, you should set it to avoid calculating it over and over. pub fn set_cov(&mut self) -> I { let cov = self.calculate_coverage(); @@ -358,7 +488,7 @@ where } } - /// Merge any intervals that overlap with eachother within the Lapper. This is an easy way to + /// Merge any intervals that overlap with each other within the Lapper. This is an easy way to /// speed up queries. pub fn merge_overlaps(&mut self) { let mut stack: VecDeque<&mut Interval> = VecDeque::new(); @@ -389,30 +519,7 @@ where }) .collect(); } - // Fix the starts and stops used by counts - let (mut starts, mut stops): (Vec<_>, Vec<_>) = - self.intervals.iter().map(|x| (x.start, x.stop)).unzip(); - - #[cfg(feature = "sort_unstable")] - { - starts.sort_unstable(); - stops.sort_unstable(); - } - - #[cfg(not(feature = "sort_unstable"))] - { - starts.sort(); - stops.sort(); - } - - self.starts = starts; - self.stops = stops; - self.max_len = self - .intervals - .iter() - .map(|x| x.stop.checked_sub(&x.start).unwrap_or_else(zero::)) - .max() - .unwrap_or_else(zero::); + self.rebuild_derived(); } /// Determine the first index that we should start checking for overlaps for via a binary @@ -465,15 +572,15 @@ where cursor } - /// Find the union and the intersect of two lapper objects. - /// Union: The set of positions found in both lappers - /// Intersect: The number of positions where both lappers intersect. Note that a position only - /// counts one time, multiple Intervals covering the same position don't add up. + /// Return the number of positions in the union and intersection of two Lappers. + /// + /// The union counts each position covered by either Lapper once. The intersection counts each + /// position covered by both Lappers once, regardless of how many intervals cover it. /// ``` rust /// use rust_lapper::{Lapper, Interval}; /// type Iv = Interval; /// let data1: Vec = vec![ - /// Iv{start: 70, stop: 120, val: 0}, // max_len = 50 + /// Iv{start: 70, stop: 120, val: 0}, // a long interval /// Iv{start: 10, stop: 15, val: 0}, // exact overlap /// Iv{start: 12, stop: 15, val: 0}, // inner overlap /// Iv{start: 14, stop: 16, val: 0}, // overlap end @@ -591,13 +698,15 @@ where inner: self, merged: merged_lapper, curr_merged_pos: zero::(), + initialized: false, curr_pos: 0, cursor: 0, end: merged_len, } } - /// Count all intervals that overlap start .. stop. This performs two binary search in order to + /// Count all intervals that overlap the half-open query `[start, stop)`. + /// This performs two binary searches in order to /// find all the excluded elements, and then deduces the intersection from there. See /// [BITS](https://arxiv.org/pdf/1208.3407.pdf) for more details. /// ``` @@ -610,14 +719,15 @@ where #[inline] pub fn count(&self, start: I, stop: I) -> usize { let len = self.intervals.len(); - // Plus one to account for half-openness of lapper intervals compared to BITS paper - let first = Self::bsearch_seq(start + one::(), &self.stops); + let first = self + .stops + .partition_point(|interval_stop| *interval_stop <= start); let last = Self::bsearch_seq(stop, &self.starts); let num_cant_after = len - last; len - first - num_cant_after } - /// Find all intervals that overlap start .. stop + /// Find all intervals that overlap the half-open query `[start, stop)`. /// ``` /// use rust_lapper::{Lapper, Interval}; /// let lapper = Lapper::new((0..100).step_by(5) @@ -627,22 +737,29 @@ where /// ``` #[inline] pub fn find(&self, start: I, stop: I) -> IterFind<'_, I, T> { + let off = self + .block_prefix_max_ends + .partition_point(|max_end| *max_end <= start) + * INDEX_BLOCK_SIZE; IterFind { inner: self, - off: Self::lower_bound( - start.checked_sub(&self.max_len).unwrap_or_else(zero::), - &self.intervals, - ), + next_block_start: off, + mask_block_start: 0, + mask: 0, + dense_next: 0, + dense_end: 0, + backend: detect_backend(), start, stop, } } - /// Find all intevals that overlap start .. stop. This method will work when queries - /// to this lapper are in sorted (start) order. It uses a linear search from the last query - /// instead of a binary search. A reference to a cursor must be passed in. This reference will - /// be modified and should be reused in the next query. This allows seek to not need to make - /// the lapper object mutable, and thus use the same lapper accross threads. + /// Find all intervals that overlap the half-open query `[start, stop)`. + /// + /// Use this method when query starts arrive in nondecreasing order. A caller-owned cursor + /// narrows the first candidate block, after which `seek()` uses the same block traversal as + /// [`Lapper::find`]. Keeping the cursor outside `Lapper` allows immutable queries and preserves + /// `Sync` when `T` and `I` are `Sync`. /// ``` /// use rust_lapper::{Lapper, Interval}; /// let lapper = Lapper::new((0..100).step_by(5) @@ -655,24 +772,36 @@ where /// ``` #[inline] pub fn seek<'a>(&'a self, start: I, stop: I, cursor: &mut usize) -> IterFind<'a, I, T> { - if *cursor == 0 || (*cursor < self.intervals.len() && self.intervals[*cursor].start > start) - { - *cursor = Self::lower_bound( - start.checked_sub(&self.max_len).unwrap_or_else(zero::), - &self.intervals, - ); - } + if self.max_len_overflowed { + *cursor = self + .block_prefix_max_ends + .partition_point(|max_end| *max_end <= start) + * INDEX_BLOCK_SIZE; + } else { + let earliest_start = start + .checked_sub(&self.max_len) + .unwrap_or_else(I::min_value); + if *cursor == 0 + || (*cursor < self.intervals.len() && self.intervals[*cursor].start > start) + { + *cursor = Self::lower_bound(earliest_start, &self.intervals); + } - while *cursor + 1 < self.intervals.len() - && self.intervals[*cursor + 1].start - < start.checked_sub(&self.max_len).unwrap_or_else(zero::) - { - *cursor += 1; + while *cursor + 1 < self.intervals.len() + && self.intervals[*cursor + 1].start < earliest_start + { + *cursor += 1; + } } IterFind { inner: self, - off: *cursor, + next_block_start: (*cursor / INDEX_BLOCK_SIZE) * INDEX_BLOCK_SIZE, + mask_block_start: 0, + mask: 0, + dense_next: 0, + dense_end: 0, + backend: detect_backend(), start, stop, } @@ -684,36 +813,119 @@ where pub struct IterFind<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { inner: &'a Lapper, - off: usize, + next_block_start: usize, + mask_block_start: usize, + mask: u32, + dense_next: usize, + dense_end: usize, + backend: MaskBackend, start: I, stop: I, } +impl<'a, I, T> IterFind<'a, I, T> +where + T: Eq + Clone + Send + Sync + 'a, + I: PrimInt + Ord + Clone + Send + Sync + 'static, +{ + #[inline(always)] + fn next_blockwise(&mut self) -> Option<&'a Interval> { + loop { + // Return the next pending match from the dense branch. + if self.dense_next < self.dense_end { + let index = self.dense_next; + self.dense_next += 1; + return Some(&self.inner.intervals[index]); + } + + // Return the next pending match from a mixed block's mask. + if self.mask != 0 { + let lane = self.mask.trailing_zeros() as usize; + self.mask &= self.mask - 1; + return Some(&self.inner.intervals[self.mask_block_start + lane]); + } + + let block_start = self.next_block_start; + // No blocks remain, so there are no more matches. + if block_start >= self.inner.starts.len() { + return None; + } + + // Private arrays are rebuilt together by constructors, mutations, + // and deserialization. Public interval edits cannot extend this bound. + // Starts are sorted, so no interval in this or a later block can match. + if unsafe { *self.inner.starts.get_unchecked(block_start) } >= self.stop { + return None; + } + + let block = block_start / INDEX_BLOCK_SIZE; + // No ends in this block reach past the query start. Skip to the next + // block whose maximum end might match. + if unsafe { *self.inner.block_max_ends.get_unchecked(block) } <= self.start { + let next_block = unsafe { *self.inner.block_index.get_unchecked(block) }; + debug_assert!( + next_block <= self.inner.block_index.len(), + "Lapper block index is corrupt" + ); + self.next_block_start = next_block * INDEX_BLOCK_SIZE; + continue; + } + + let block_end = (block_start + INDEX_BLOCK_SIZE).min(self.inner.starts.len()); + // All ends match. Take the dense branch for the prefix whose starts + // are before the query stop. + if unsafe { *self.inner.block_min_ends.get_unchecked(block) } > self.start { + let starts = unsafe { self.inner.starts.get_unchecked(block_start..block_end) }; + // If the last start matches, the whole block is dense. Otherwise, + // find the matching prefix. + let active_len = if unsafe { *starts.get_unchecked(starts.len() - 1) } < self.stop { + starts.len() + } else { + starts.partition_point(|lane_start| *lane_start < self.stop) + }; + self.dense_next = block_start; + self.dense_end = block_start + active_len; + // Continue after a full block. A partial prefix means every later + // start is outside the query. + self.next_block_start = if active_len == starts.len() { + block_end + } else { + self.inner.starts.len() + }; + continue; + } + + // Some ends match. Build the exact overlap mask for this mixed block. + self.mask_block_start = block_start; + self.next_block_start = block_end; + self.mask = overlap_mask( + self.backend, + unsafe { self.inner.starts.get_unchecked(block_start..block_end) }, + unsafe { + self.inner + .stops_by_start + .get_unchecked(block_start..block_end) + }, + self.start, + self.stop, + ); + } + } +} + impl<'a, I, T> Iterator for IterFind<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync + 'static, { type Item = &'a Interval; #[inline] - // interval.start < stop && interval.stop > start fn next(&mut self) -> Option { - while self.off < self.inner.intervals.len() { - //let mut generator = self.inner.intervals[self.off..].iter(); - //while let Some(interval) = generator.next() { - let interval = &self.inner.intervals[self.off]; - self.off += 1; - if interval.overlap(self.start, self.stop) { - return Some(interval); - } else if interval.start >= self.stop { - break; - } - } - None + self.next_blockwise() } } @@ -722,28 +934,30 @@ where pub struct IterDepth<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { inner: &'a Lapper, merged: Lapper, // A lapper that is the merged_lapper of inner curr_merged_pos: I, // Current start position in current interval - curr_pos: usize, // In merged list of non-overlapping intervals - cursor: usize, // cursor for seek over inner lapper - end: usize, // len of merged + initialized: bool, + curr_pos: usize, // In merged list of non-overlapping intervals + cursor: usize, // cursor for seek over inner lapper + end: usize, // len of merged } impl<'a, I, T> Iterator for IterDepth<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync + 'static, { type Item = Interval; #[inline] fn next(&mut self) -> Option { let mut interval: &Interval = &self.merged.intervals[self.curr_pos]; - if self.curr_merged_pos == zero::() { + if !self.initialized { self.curr_merged_pos = interval.start; + self.initialized = true; } if interval.stop == self.curr_merged_pos { if self.curr_pos + 1 != self.end { @@ -766,6 +980,9 @@ where let mut new_depth_at_point = depth_at_point; while new_depth_at_point == depth_at_point && self.curr_merged_pos < interval.stop { self.curr_merged_pos = self.curr_merged_pos + one::(); + if self.curr_merged_pos == interval.stop { + break; + } new_depth_at_point = self .inner .seek( @@ -786,7 +1003,7 @@ where pub struct IterLapper<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { inner: &'a Lapper, pos: usize, @@ -795,7 +1012,7 @@ where impl<'a, I, T> Iterator for IterLapper<'a, I, T> where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { type Item = &'a Interval; @@ -812,7 +1029,7 @@ where impl IntoIterator for Lapper where T: Eq + Clone + Send + Sync, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { type Item = Interval; type IntoIter = ::std::vec::IntoIter; @@ -825,7 +1042,7 @@ where impl<'a, I, T> IntoIterator for &'a Lapper where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { type Item = &'a Interval; type IntoIter = std::slice::Iter<'a, Interval>; @@ -838,7 +1055,7 @@ where impl<'a, I, T> IntoIterator for &'a mut Lapper where T: Eq + Clone + Send + Sync + 'a, - I: PrimInt + Unsigned + Ord + Clone + Send + Sync, + I: PrimInt + Ord + Clone + Send + Sync, { type Item = &'a mut Interval; type IntoIter = std::slice::IterMut<'a, Interval>; @@ -1363,7 +1580,7 @@ mod tests { } // When there is a very long interval that spans many little intervals, test that the little - // intevals still get returne properly + // Intervals still get returned properly. #[test] fn test_bad_skips() { let data = vec![ @@ -1384,9 +1601,20 @@ mod tests { assert_eq!(lapper.count(28974798, 33141355), 1); } + #[cfg(feature = "with_serde")] + #[derive(Debug, PartialEq, serde::Serialize, serde::Deserialize)] + struct LegacyLapper { + intervals: Vec, + starts: Vec, + stops: Vec, + max_len: usize, + cov: Option, + overlaps_merged: bool, + } + #[cfg(feature = "with_serde")] #[test] - fn serde_test() { + fn serde_keeps_the_v1_six_field_representation() { let data = vec![ Iv{start:25264912, stop: 25264986, val: 0}, Iv{start:27273024, stop: 27273065 , val: 0}, @@ -1398,14 +1626,26 @@ mod tests { ]; let lapper = Lapper::new(data); - let serialized = bincode::serialize(&lapper).unwrap(); - let deserialzed: Lapper = bincode::deserialize(&serialized).unwrap(); - - let found = deserialzed.find(28974798, 33141355).collect::>(); + let legacy = LegacyLapper { + intervals: lapper.intervals.clone(), + starts: lapper.starts.clone(), + stops: lapper.stops.clone(), + max_len: lapper.max_len, + cov: lapper.cov, + overlaps_merged: lapper.overlaps_merged, + }; + let legacy_bytes = bincode::serialize(&legacy).unwrap(); + let deserialized: Lapper = bincode::deserialize(&legacy_bytes).unwrap(); + let current_bytes = bincode::serialize(&deserialized).unwrap(); + assert_eq!(current_bytes, legacy_bytes); + let legacy_again: LegacyLapper = bincode::deserialize(¤t_bytes).unwrap(); + assert_eq!(legacy_again, legacy); + + let found = deserialized.find(28974798, 33141355).collect::>(); assert_eq!(found, vec![ &Iv{start:28866309, stop: 33141404 , val: 0}, ]); - assert_eq!(deserialzed.count(28974798, 33141355), 1); + assert_eq!(deserialized.count(28974798, 33141355), 1); } } diff --git a/src/simd.rs b/src/simd.rs new file mode 100644 index 0000000..376a358 --- /dev/null +++ b/src/simd.rs @@ -0,0 +1,789 @@ +use num_traits::PrimInt; +#[cfg(any(target_arch = "aarch64", target_arch = "x86_64"))] +use std::any::TypeId; + +pub(crate) const BLOCK_SIZE: usize = 32; + +#[derive(Clone, Copy, Debug)] +pub(crate) enum MaskBackend { + Scalar, + #[cfg(target_arch = "aarch64")] + Neon, + #[cfg(target_arch = "x86_64")] + Avx2, +} + +#[inline] +pub(crate) fn detect_backend() -> MaskBackend { + #[cfg(target_arch = "aarch64")] + { + return MaskBackend::Neon; + } + + #[cfg(target_arch = "x86_64")] + { + if std::is_x86_feature_detected!("avx2") { + return MaskBackend::Avx2; + } + } + + #[allow(unreachable_code)] + MaskBackend::Scalar +} + +#[inline(always)] +pub(crate) fn overlap_mask( + backend: MaskBackend, + starts: &[I], + stops: &[I], + query_start: I, + query_stop: I, +) -> u32 +where + I: PrimInt + 'static, +{ + debug_assert_eq!(starts.len(), stops.len()); + debug_assert!(starts.len() <= BLOCK_SIZE); + let _ = backend; + + #[cfg(target_arch = "aarch64")] + if matches!(backend, MaskBackend::Neon) { + macro_rules! neon_as { + ($source:ty, $repr:ty, $function:path) => { + if TypeId::of::() == TypeId::of::<$source>() { + // TypeId proves the source type. On AArch64, usize/isize have the + // same representation and alignment as u64/i64 respectively. + let starts = unsafe { + std::slice::from_raw_parts(starts.as_ptr().cast::<$repr>(), starts.len()) + }; + let stops = unsafe { + std::slice::from_raw_parts(stops.as_ptr().cast::<$repr>(), stops.len()) + }; + let query_start = unsafe { *(&query_start as *const I).cast::<$repr>() }; + let query_stop = unsafe { *(&query_stop as *const I).cast::<$repr>() }; + return unsafe { $function(starts, stops, query_start, query_stop) }; + } + }; + } + + neon_as!(u8, u8, neon::mask_u8); + neon_as!(i8, i8, neon::mask_i8); + neon_as!(u16, u16, neon::mask_u16); + neon_as!(i16, i16, neon::mask_i16); + neon_as!(u32, u32, neon::mask_u32); + neon_as!(i32, i32, neon::mask_i32); + neon_as!(u64, u64, neon::mask_u64); + neon_as!(i64, i64, neon::mask_i64); + neon_as!(usize, u64, neon::mask_u64); + neon_as!(isize, i64, neon::mask_i64); + } + + #[cfg(target_arch = "x86_64")] + if matches!(backend, MaskBackend::Avx2) { + macro_rules! avx2_as { + ($source:ty, $repr:ty, $function:path) => { + if TypeId::of::() == TypeId::of::<$source>() { + // TypeId proves the source type. On x86-64, usize/isize have the + // same representation and alignment as u64/i64 respectively. + let starts = unsafe { + std::slice::from_raw_parts(starts.as_ptr().cast::<$repr>(), starts.len()) + }; + let stops = unsafe { + std::slice::from_raw_parts(stops.as_ptr().cast::<$repr>(), stops.len()) + }; + let query_start = unsafe { *(&query_start as *const I).cast::<$repr>() }; + let query_stop = unsafe { *(&query_stop as *const I).cast::<$repr>() }; + return unsafe { $function(starts, stops, query_start, query_stop) }; + } + }; + } + + avx2_as!(u8, u8, avx2::mask_u8); + avx2_as!(i8, i8, avx2::mask_i8); + avx2_as!(u16, u16, avx2::mask_u16); + avx2_as!(i16, i16, avx2::mask_i16); + avx2_as!(u32, u32, avx2::mask_u32); + avx2_as!(i32, i32, avx2::mask_i32); + avx2_as!(u64, u64, avx2::mask_u64); + avx2_as!(i64, i64, avx2::mask_i64); + avx2_as!(usize, u64, avx2::mask_u64); + avx2_as!(isize, i64, avx2::mask_i64); + } + + scalar_mask(starts, stops, query_start, query_stop) +} + +#[inline(always)] +fn scalar_mask(starts: &[I], stops: &[I], query_start: I, query_stop: I) -> u32 { + let mut mask = 0_u32; + for lane in 0..starts.len() { + if stops[lane] > query_start && starts[lane] < query_stop { + mask |= 1 << lane; + } + } + mask +} + +#[cfg(test)] +mod dispatch_tests { + use super::*; + + #[test] + fn selected_backend_matches_the_host() { + let backend = detect_backend(); + + #[cfg(target_arch = "aarch64")] + { + assert!(matches!(backend, MaskBackend::Neon)); + eprintln!("runtime backend: NEON"); + } + + #[cfg(target_arch = "x86_64")] + if std::is_x86_feature_detected!("avx2") { + assert!(matches!(backend, MaskBackend::Avx2)); + eprintln!("runtime backend: AVX2"); + } else { + assert!(matches!(backend, MaskBackend::Scalar)); + eprintln!("runtime backend: scalar"); + } + + #[cfg(not(any(target_arch = "aarch64", target_arch = "x86_64")))] + assert!(matches!(backend, MaskBackend::Scalar)); + } +} + +#[cfg(target_arch = "aarch64")] +mod neon { + use std::arch::aarch64::*; + + #[inline(always)] + unsafe fn bits_u8x16(overlap: uint8x16_t) -> u16 { + let weights = vld1_u8([1_u8, 2, 4, 8, 16, 32, 64, 128].as_ptr()); + let low = vaddv_u8(vand_u8(vget_low_u8(overlap), weights)); + let high = vaddv_u8(vand_u8(vget_high_u8(overlap), weights)); + u16::from(low) | (u16::from(high) << 8) + } + + #[inline(always)] + unsafe fn bits_u16x8(overlap: uint16x8_t) -> u8 { + let weights = vld1_u8([1_u8, 2, 4, 8, 16, 32, 64, 128].as_ptr()); + vaddv_u8(vand_u8(vmovn_u16(overlap), weights)) + } + + #[inline(always)] + unsafe fn bits_u16x16(low: uint16x8_t, high: uint16x8_t) -> u16 { + bits_u8x16(vcombine_u8(vmovn_u16(low), vmovn_u16(high))) + } + + #[inline(always)] + unsafe fn bits_u32x4(overlap: uint32x4_t) -> u32 { + let weights = vld1q_u32([1_u32, 2, 4, 8].as_ptr()); + vaddvq_u32(vandq_u32(overlap, weights)) + } + + #[inline(always)] + unsafe fn bits_u32x8(low: uint32x4_t, high: uint32x4_t) -> u16 { + let weights = vld1q_u16([1_u16, 2, 4, 8, 16, 32, 64, 128].as_ptr()); + vaddvq_u16(vandq_u16( + vcombine_u16(vmovn_u32(low), vmovn_u32(high)), + weights, + )) + } + + #[inline(always)] + unsafe fn bits_u64x4(low: uint64x2_t, high: uint64x2_t) -> u32 { + let weights = vld1q_u32([1_u32, 2, 4, 8].as_ptr()); + vaddvq_u32(vandq_u32( + vcombine_u32(vmovn_u64(low), vmovn_u64(high)), + weights, + )) + } + + #[inline(always)] + unsafe fn scalar_tail( + starts: &[I], + stops: &[I], + query_start: I, + query_stop: I, + mut lane: usize, + mut mask: u32, + ) -> u32 { + while lane < starts.len() { + if *stops.get_unchecked(lane) > query_start && *starts.get_unchecked(lane) < query_stop + { + mask |= 1 << lane; + } + lane += 1; + } + mask + } + + pub(super) unsafe fn mask_u8( + starts: &[u8], + stops: &[u8], + query_start: u8, + query_stop: u8, + ) -> u32 { + let query_start_v = vdupq_n_u8(query_start); + let query_stop_v = vdupq_n_u8(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 32 <= starts.len() { + let lane_starts = vld1q_u8_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u8_x2(stops.as_ptr().add(lane)); + let low = vandq_u8( + vcgtq_u8(lane_stops.0, query_start_v), + vcgtq_u8(query_stop_v, lane_starts.0), + ); + let high = vandq_u8( + vcgtq_u8(lane_stops.1, query_start_v), + vcgtq_u8(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u8x16(low)) << lane; + mask |= u32::from(bits_u8x16(high)) << (lane + 16); + lane += 32; + } + while lane + 16 <= starts.len() { + let lane_starts = vld1q_u8(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u8(stops.as_ptr().add(lane)); + let overlap = vandq_u8( + vcgtq_u8(lane_stops, query_start_v), + vcgtq_u8(query_stop_v, lane_starts), + ); + mask |= u32::from(bits_u8x16(overlap)) << lane; + lane += 16; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_i8( + starts: &[i8], + stops: &[i8], + query_start: i8, + query_stop: i8, + ) -> u32 { + let query_start_v = vdupq_n_s8(query_start); + let query_stop_v = vdupq_n_s8(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 32 <= starts.len() { + let lane_starts = vld1q_s8_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s8_x2(stops.as_ptr().add(lane)); + let low = vandq_u8( + vcgtq_s8(lane_stops.0, query_start_v), + vcgtq_s8(query_stop_v, lane_starts.0), + ); + let high = vandq_u8( + vcgtq_s8(lane_stops.1, query_start_v), + vcgtq_s8(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u8x16(low)) << lane; + mask |= u32::from(bits_u8x16(high)) << (lane + 16); + lane += 32; + } + while lane + 16 <= starts.len() { + let lane_starts = vld1q_s8(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s8(stops.as_ptr().add(lane)); + let overlap = vandq_u8( + vcgtq_s8(lane_stops, query_start_v), + vcgtq_s8(query_stop_v, lane_starts), + ); + mask |= u32::from(bits_u8x16(overlap)) << lane; + lane += 16; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_u16( + starts: &[u16], + stops: &[u16], + query_start: u16, + query_stop: u16, + ) -> u32 { + let query_start_v = vdupq_n_u16(query_start); + let query_stop_v = vdupq_n_u16(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 16 <= starts.len() { + let lane_starts = vld1q_u16_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u16_x2(stops.as_ptr().add(lane)); + let low = vandq_u16( + vcgtq_u16(lane_stops.0, query_start_v), + vcgtq_u16(query_stop_v, lane_starts.0), + ); + let high = vandq_u16( + vcgtq_u16(lane_stops.1, query_start_v), + vcgtq_u16(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u16x16(low, high)) << lane; + lane += 16; + } + while lane + 8 <= starts.len() { + let lane_starts = vld1q_u16(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u16(stops.as_ptr().add(lane)); + let overlap = vandq_u16( + vcgtq_u16(lane_stops, query_start_v), + vcgtq_u16(query_stop_v, lane_starts), + ); + mask |= u32::from(bits_u16x8(overlap)) << lane; + lane += 8; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_i16( + starts: &[i16], + stops: &[i16], + query_start: i16, + query_stop: i16, + ) -> u32 { + let query_start_v = vdupq_n_s16(query_start); + let query_stop_v = vdupq_n_s16(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 16 <= starts.len() { + let lane_starts = vld1q_s16_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s16_x2(stops.as_ptr().add(lane)); + let low = vandq_u16( + vcgtq_s16(lane_stops.0, query_start_v), + vcgtq_s16(query_stop_v, lane_starts.0), + ); + let high = vandq_u16( + vcgtq_s16(lane_stops.1, query_start_v), + vcgtq_s16(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u16x16(low, high)) << lane; + lane += 16; + } + while lane + 8 <= starts.len() { + let lane_starts = vld1q_s16(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s16(stops.as_ptr().add(lane)); + let overlap = vandq_u16( + vcgtq_s16(lane_stops, query_start_v), + vcgtq_s16(query_stop_v, lane_starts), + ); + mask |= u32::from(bits_u16x8(overlap)) << lane; + lane += 8; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_u32( + starts: &[u32], + stops: &[u32], + query_start: u32, + query_stop: u32, + ) -> u32 { + let query_start_v = vdupq_n_u32(query_start); + let query_stop_v = vdupq_n_u32(query_stop); + let simd_len = starts.len() / 4 * 4; + let mut lane = 0; + let mut mask = 0; + while lane + 8 <= simd_len { + let lane_starts = vld1q_u32_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u32_x2(stops.as_ptr().add(lane)); + let low = vandq_u32( + vcgtq_u32(lane_stops.0, query_start_v), + vcgtq_u32(query_stop_v, lane_starts.0), + ); + let high = vandq_u32( + vcgtq_u32(lane_stops.1, query_start_v), + vcgtq_u32(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u32x8(low, high)) << lane; + lane += 8; + } + while lane < simd_len { + let lane_starts = vld1q_u32(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u32(stops.as_ptr().add(lane)); + let overlap = vandq_u32( + vcgtq_u32(lane_stops, query_start_v), + vcgtq_u32(query_stop_v, lane_starts), + ); + mask |= bits_u32x4(overlap) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_i32( + starts: &[i32], + stops: &[i32], + query_start: i32, + query_stop: i32, + ) -> u32 { + let query_start_v = vdupq_n_s32(query_start); + let query_stop_v = vdupq_n_s32(query_stop); + let simd_len = starts.len() / 4 * 4; + let mut lane = 0; + let mut mask = 0; + while lane + 8 <= simd_len { + let lane_starts = vld1q_s32_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s32_x2(stops.as_ptr().add(lane)); + let low = vandq_u32( + vcgtq_s32(lane_stops.0, query_start_v), + vcgtq_s32(query_stop_v, lane_starts.0), + ); + let high = vandq_u32( + vcgtq_s32(lane_stops.1, query_start_v), + vcgtq_s32(query_stop_v, lane_starts.1), + ); + mask |= u32::from(bits_u32x8(low, high)) << lane; + lane += 8; + } + while lane < simd_len { + let lane_starts = vld1q_s32(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s32(stops.as_ptr().add(lane)); + let overlap = vandq_u32( + vcgtq_s32(lane_stops, query_start_v), + vcgtq_s32(query_stop_v, lane_starts), + ); + mask |= bits_u32x4(overlap) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_u64( + starts: &[u64], + stops: &[u64], + query_start: u64, + query_stop: u64, + ) -> u32 { + let query_start_v = vdupq_n_u64(query_start); + let query_stop_v = vdupq_n_u64(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 4 <= starts.len() { + let lane_starts = vld1q_u64_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_u64_x2(stops.as_ptr().add(lane)); + let low = vandq_u64( + vcgtq_u64(lane_stops.0, query_start_v), + vcgtq_u64(query_stop_v, lane_starts.0), + ); + let high = vandq_u64( + vcgtq_u64(lane_stops.1, query_start_v), + vcgtq_u64(query_stop_v, lane_starts.1), + ); + mask |= bits_u64x4(low, high) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + pub(super) unsafe fn mask_i64( + starts: &[i64], + stops: &[i64], + query_start: i64, + query_stop: i64, + ) -> u32 { + let query_start_v = vdupq_n_s64(query_start); + let query_stop_v = vdupq_n_s64(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 4 <= starts.len() { + let lane_starts = vld1q_s64_x2(starts.as_ptr().add(lane)); + let lane_stops = vld1q_s64_x2(stops.as_ptr().add(lane)); + let low = vandq_u64( + vcgtq_s64(lane_stops.0, query_start_v), + vcgtq_s64(query_stop_v, lane_starts.0), + ); + let high = vandq_u64( + vcgtq_s64(lane_stops.1, query_start_v), + vcgtq_s64(query_stop_v, lane_starts.1), + ); + mask |= bits_u64x4(low, high) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } +} + +#[cfg(target_arch = "x86_64")] +mod avx2 { + use std::arch::x86_64::*; + + #[inline(always)] + fn compact_u16_movemask(mut bits: u32) -> u32 { + bits &= 0x5555_5555; + bits = (bits | (bits >> 1)) & 0x3333_3333; + bits = (bits | (bits >> 2)) & 0x0f0f_0f0f; + bits = (bits | (bits >> 4)) & 0x00ff_00ff; + (bits | (bits >> 8)) & 0x0000_ffff + } + + #[inline(always)] + unsafe fn scalar_tail( + starts: &[I], + stops: &[I], + query_start: I, + query_stop: I, + mut lane: usize, + mut mask: u32, + ) -> u32 { + while lane < starts.len() { + if *stops.get_unchecked(lane) > query_start && *starts.get_unchecked(lane) < query_stop + { + mask |= 1 << lane; + } + lane += 1; + } + mask + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_u8( + starts: &[u8], + stops: &[u8], + query_start: u8, + query_stop: u8, + ) -> u32 { + let bias = _mm256_set1_epi8(i8::MIN); + let query_start_v = _mm256_xor_si256(_mm256_set1_epi8(query_start as i8), bias); + let query_stop_v = _mm256_xor_si256(_mm256_set1_epi8(query_stop as i8), bias); + let mut lane = 0; + let mut mask = 0; + while lane + 32 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi8(_mm256_xor_si256(lane_stops, bias), query_start_v), + _mm256_cmpgt_epi8(query_stop_v, _mm256_xor_si256(lane_starts, bias)), + ); + mask |= (_mm256_movemask_epi8(overlap) as u32) << lane; + lane += 32; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_i8( + starts: &[i8], + stops: &[i8], + query_start: i8, + query_stop: i8, + ) -> u32 { + let query_start_v = _mm256_set1_epi8(query_start); + let query_stop_v = _mm256_set1_epi8(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 32 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi8(lane_stops, query_start_v), + _mm256_cmpgt_epi8(query_stop_v, lane_starts), + ); + mask |= (_mm256_movemask_epi8(overlap) as u32) << lane; + lane += 32; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_u16( + starts: &[u16], + stops: &[u16], + query_start: u16, + query_stop: u16, + ) -> u32 { + let bias = _mm256_set1_epi16(i16::MIN); + let query_start_v = _mm256_xor_si256(_mm256_set1_epi16(query_start as i16), bias); + let query_stop_v = _mm256_xor_si256(_mm256_set1_epi16(query_stop as i16), bias); + let mut lane = 0; + let mut mask = 0; + while lane + 16 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi16(_mm256_xor_si256(lane_stops, bias), query_start_v), + _mm256_cmpgt_epi16(query_stop_v, _mm256_xor_si256(lane_starts, bias)), + ); + mask |= compact_u16_movemask(_mm256_movemask_epi8(overlap) as u32) << lane; + lane += 16; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_i16( + starts: &[i16], + stops: &[i16], + query_start: i16, + query_stop: i16, + ) -> u32 { + let query_start_v = _mm256_set1_epi16(query_start); + let query_stop_v = _mm256_set1_epi16(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 16 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi16(lane_stops, query_start_v), + _mm256_cmpgt_epi16(query_stop_v, lane_starts), + ); + mask |= compact_u16_movemask(_mm256_movemask_epi8(overlap) as u32) << lane; + lane += 16; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_u32( + starts: &[u32], + stops: &[u32], + query_start: u32, + query_stop: u32, + ) -> u32 { + let bias = _mm256_set1_epi32(i32::MIN); + let query_start_v = _mm256_xor_si256(_mm256_set1_epi32(query_start as i32), bias); + let query_stop_v = _mm256_xor_si256(_mm256_set1_epi32(query_stop as i32), bias); + let mut lane = 0; + let mut mask = 0; + while lane + 8 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi32(_mm256_xor_si256(lane_stops, bias), query_start_v), + _mm256_cmpgt_epi32(query_stop_v, _mm256_xor_si256(lane_starts, bias)), + ); + mask |= (_mm256_movemask_ps(_mm256_castsi256_ps(overlap)) as u32) << lane; + lane += 8; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_i32( + starts: &[i32], + stops: &[i32], + query_start: i32, + query_stop: i32, + ) -> u32 { + let query_start_v = _mm256_set1_epi32(query_start); + let query_stop_v = _mm256_set1_epi32(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 8 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi32(lane_stops, query_start_v), + _mm256_cmpgt_epi32(query_stop_v, lane_starts), + ); + mask |= (_mm256_movemask_ps(_mm256_castsi256_ps(overlap)) as u32) << lane; + lane += 8; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_u64( + starts: &[u64], + stops: &[u64], + query_start: u64, + query_stop: u64, + ) -> u32 { + let bias = _mm256_set1_epi64x(i64::MIN); + let query_start_v = _mm256_xor_si256(_mm256_set1_epi64x(query_start as i64), bias); + let query_stop_v = _mm256_xor_si256(_mm256_set1_epi64x(query_stop as i64), bias); + let mut lane = 0; + let mut mask = 0; + while lane + 4 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi64(_mm256_xor_si256(lane_stops, bias), query_start_v), + _mm256_cmpgt_epi64(query_stop_v, _mm256_xor_si256(lane_starts, bias)), + ); + mask |= (_mm256_movemask_pd(_mm256_castsi256_pd(overlap)) as u32) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[target_feature(enable = "avx2")] + pub(super) unsafe fn mask_i64( + starts: &[i64], + stops: &[i64], + query_start: i64, + query_stop: i64, + ) -> u32 { + let query_start_v = _mm256_set1_epi64x(query_start); + let query_stop_v = _mm256_set1_epi64x(query_stop); + let mut lane = 0; + let mut mask = 0; + while lane + 4 <= starts.len() { + let lane_starts = _mm256_loadu_si256(starts.as_ptr().add(lane).cast()); + let lane_stops = _mm256_loadu_si256(stops.as_ptr().add(lane).cast()); + let overlap = _mm256_and_si256( + _mm256_cmpgt_epi64(lane_stops, query_start_v), + _mm256_cmpgt_epi64(query_stop_v, lane_starts), + ); + mask |= (_mm256_movemask_pd(_mm256_castsi256_pd(overlap)) as u32) << lane; + lane += 4; + } + scalar_tail(starts, stops, query_start, query_stop, lane, mask) + } + + #[cfg(test)] + mod tests { + use super::*; + + #[test] + fn u16_movemask_compaction_preserves_every_lane_bit() { + for expected in 0_u32..=u16::MAX.into() { + let mut expanded = 0_u32; + for lane in 0..16 { + if expected & (1 << lane) != 0 { + expanded |= 0b11 << (lane * 2); + } + } + assert_eq!(compact_u16_movemask(expanded), expected); + } + } + + #[test] + fn primitive_masks_match_scalar_when_avx2_is_available() { + if !std::is_x86_feature_detected!("avx2") { + return; + } + + macro_rules! check_unsigned { + ($type:ty, $function:path) => {{ + let starts: Vec<$type> = (0..32).map(|lane| (lane * 4) as $type).collect(); + let stops: Vec<$type> = starts.iter().map(|start| *start + 20).collect(); + let expected = + super::super::scalar_mask(&starts, &stops, 37 as $type, 91 as $type); + assert_eq!( + unsafe { $function(&starts, &stops, 37 as $type, 91 as $type) }, + expected + ); + }}; + } + + macro_rules! check_signed { + ($type:ty, $function:path) => {{ + let starts: Vec<$type> = (0..32).map(|lane| (lane * 4 - 64) as $type).collect(); + let stops: Vec<$type> = starts.iter().map(|start| *start + 20).collect(); + let expected = + super::super::scalar_mask(&starts, &stops, -11 as $type, 43 as $type); + assert_eq!( + unsafe { $function(&starts, &stops, -11 as $type, 43 as $type) }, + expected + ); + }}; + } + + check_unsigned!(u8, mask_u8); + check_unsigned!(u16, mask_u16); + check_unsigned!(u32, mask_u32); + check_unsigned!(u64, mask_u64); + check_signed!(i8, mask_i8); + check_signed!(i16, mask_i16); + check_signed!(i32, mask_i32); + check_signed!(i64, mask_i64); + } + } +} diff --git a/tests/block_index.rs b/tests/block_index.rs new file mode 100644 index 0000000..071c537 --- /dev/null +++ b/tests/block_index.rs @@ -0,0 +1,40 @@ +use rust_lapper::{Interval, Lapper}; + +#[test] +fn block_index_matches_forward_brute_force() { + let mut state = 0x1234_5678_u64; + let mut next = || { + state = state.wrapping_mul(6364136223846793005).wrapping_add(1); + (state >> 32) as u32 + }; + + let mut intervals = Vec::new(); + intervals.push(Interval { + start: 0, + stop: 1_000_000, + val: 0usize, + }); + for value in 1..5000 { + let start = next() % 1_000_000; + let len = 1 + next() % 1000; + intervals.push(Interval { + start, + stop: start + len, + val: value, + }); + } + let lapper = Lapper::new(intervals); + + for _ in 0..20_000 { + let start = next() % 1_000_000; + let stop = start + 1 + next() % 2000; + let got: Vec<_> = lapper.find(start, stop).map(|iv| iv.val).collect(); + let expected: Vec<_> = lapper + .intervals + .iter() + .filter(|iv| iv.start < stop && iv.stop > start) + .map(|iv| iv.val) + .collect(); + assert_eq!(got, expected, "query {start}..{stop}"); + } +} diff --git a/tests/portable_index.rs b/tests/portable_index.rs new file mode 100644 index 0000000..e03e510 --- /dev/null +++ b/tests/portable_index.rs @@ -0,0 +1,229 @@ +use num_traits::PrimInt; +use rust_lapper::{Interval, Lapper}; + +fn assert_queries_match(lapper: &Lapper, query_min: i64, query_max: i64) +where + I: PrimInt + Ord + Clone + Send + Sync + 'static, +{ + let mut cursor = 0; + for raw_start in query_min..query_max { + let start = I::from(raw_start).unwrap(); + let stop = I::from(raw_start + 7).unwrap(); + let expected: Vec<_> = lapper + .intervals + .iter() + .filter(|interval| interval.start < stop && interval.stop > start) + .map(|interval| interval.val) + .collect(); + let found: Vec<_> = lapper + .find(start, stop) + .map(|interval| interval.val) + .collect(); + let sought: Vec<_> = lapper + .seek(start, stop, &mut cursor) + .map(|interval| interval.val) + .collect(); + assert_eq!(found, expected, "find query {raw_start}..{}", raw_start + 7); + assert_eq!( + sought, + expected, + "seek query {raw_start}..{}", + raw_start + 7 + ); + assert_eq!( + lapper.count(start, stop), + expected.len(), + "count query {raw_start}..{}", + raw_start + 7 + ); + } +} + +fn exercise_unsigned() +where + I: PrimInt + Ord + Clone + Send + Sync + 'static, +{ + for interval_count in 0..=65 { + let intervals = (0..interval_count) + .map(|value| { + let raw_start = (value * 2) % 181; + let raw_stop = raw_start + 1 + (value * 7) % 19; + Interval { + start: I::from(raw_start).unwrap(), + stop: I::from(raw_stop).unwrap(), + val: value, + } + }) + .collect(); + assert_queries_match(&Lapper::new(intervals), 0, 190); + } +} + +fn exercise_signed() +where + I: PrimInt + Ord + Clone + Send + Sync + 'static, +{ + for interval_count in 0..=65 { + let intervals = (0..interval_count) + .map(|value| { + let raw_start = (value * 2) as i64 % 121 - 60; + let raw_stop = raw_start + 1 + (value * 7) as i64 % 19; + Interval { + start: I::from(raw_start).unwrap(), + stop: I::from(raw_stop).unwrap(), + val: value, + } + }) + .collect(); + assert_queries_match(&Lapper::new(intervals), -70, 70); + } +} + +#[test] +fn every_primitive_integer_type_and_block_tail_matches_forward_brute_force() { + exercise_unsigned::(); + exercise_unsigned::(); + exercise_unsigned::(); + exercise_unsigned::(); + exercise_unsigned::(); + exercise_unsigned::(); + exercise_signed::(); + exercise_signed::(); + exercise_signed::(); + exercise_signed::(); + exercise_signed::(); + exercise_signed::(); +} + +#[test] +fn insert_rebuilds_every_query_index() { + let mut lapper = Lapper::::new(Vec::new()); + for value in (0..129).rev() { + let start = (value * 13 % 257) as u32; + lapper.insert(Interval { + start, + stop: start + 1 + (value * 11 % 31) as u32, + val: value, + }); + } + assert_queries_match(&lapper, 0, 280); +} + +#[test] +fn merge_rebuilds_every_query_index() { + let intervals = (0..160) + .map(|value| { + let start = (value * 3) as u32; + Interval { + start, + stop: start + 5 + (value % 9) as u32, + val: value, + } + }) + .collect(); + let mut lapper = Lapper::new(intervals); + lapper.merge_overlaps(); + assert_queries_match(&lapper, 0, 500); +} + +#[test] +fn signed_seek_saturates_at_the_coordinate_minimum() { + let lapper = Lapper::new(vec![ + Interval { + start: i8::MIN, + stop: -100, + val: 0_usize, + }, + Interval { + start: -110, + stop: -90, + val: 1, + }, + ]); + let mut cursor = 0; + let found: Vec<_> = lapper + .seek(-125, -105, &mut cursor) + .map(|interval| interval.val) + .collect(); + assert_eq!(found, vec![0, 1]); +} + +#[test] +fn signed_seek_keeps_intervals_whose_length_exceeds_the_coordinate_type() { + let mut intervals = vec![Interval { + start: i8::MIN, + stop: i8::MAX, + val: 0_usize, + }]; + intervals.extend((-127_i16..=99).enumerate().map(|(value, raw_start)| { + let start = raw_start as i8; + Interval { + start, + stop: start + 1, + val: value + 1, + } + })); + + let lapper = Lapper::new(intervals); + let expected: Vec<_> = lapper.find(100, 101).map(|interval| interval.val).collect(); + let mut cursor = 0; + let sought: Vec<_> = lapper + .seek(100, 101, &mut cursor) + .map(|interval| interval.val) + .collect(); + + assert_eq!(expected, vec![0]); + assert_eq!(sought, expected); +} + +#[test] +fn signed_depth_crosses_zero_once() { + let lapper = Lapper::new(vec![Interval { + start: -2_i16, + stop: 3, + val: (), + }]); + let depth: Vec<_> = lapper.depth().collect(); + assert_eq!( + depth, + vec![Interval { + start: -2, + stop: 3, + val: 1 + }] + ); +} + +#[test] +fn signed_depth_stops_at_the_coordinate_maximum() { + let lapper = Lapper::new(vec![Interval { + start: 126_i8, + stop: i8::MAX, + val: (), + }]); + + assert_eq!( + lapper.depth().collect::>(), + vec![Interval { + start: 126, + stop: i8::MAX, + val: 1, + }] + ); +} + +#[test] +fn direct_structural_growth_cannot_extend_unchecked_indexing() { + let mut lapper = Lapper::new(vec![Interval { + start: 1_u32, + stop: 3, + val: 0_usize, + }]); + lapper.intervals.push(Interval { + start: 4, + stop: 6, + val: 1, + }); + let found: Vec<_> = lapper.find(0, 10).map(|interval| interval.val).collect(); + assert_eq!(found, vec![0]); +}