diff --git a/CHANGELOG.md b/CHANGELOG.md index 4838685..8b5a006 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,39 @@ Release notes for the [`caps-sa`](https://crates.io/crates/caps-sa) crate. +## Unreleased + +### Added + +- Opt-in `PackedPrefixSeedPolicy` for seeding external-memory phase-1 sorts + with segment-aware fixed-depth `u64` prefix keys. The default is + `Disabled`; `DenseAlphabetOnly` never allocates a second text-sized buffer, + while `remap(max_extra_bytes)` explicitly bounds an order-preserving ranked + copy for gapped byte alphabets. +- `LimitProvider::boundary_rank()` and `BoundaryRank` let a provider declare + whether segment ends sort below or above real symbols. This semantic + capability is separate from the `ExtMemOpts` activation policy and defaults + to `None` for custom providers. + +### Changed + +- Eligible packed-prefix builds resolve most phase-1 comparisons from one + segment-bounded key, use exact key-derived LCPs between runs, and invoke the + full comparator only inside equal-key groups. On the complete 6.56-billion- + symbol ruSTAR GRCh38 + GENCODE v50 fixture, the seed reduced phase 1 from + 49.038 to 12.204 seconds and the memoized build from 171.205 to 134.618 + seconds (21.4%) at 32 physical cores, with identical output. +- Packed-prefix eligibility is checked before alphabet scanning. Non-`u8` + symbols, finite contexts, providers without a representable boundary order, + and over-budget remaps fall back to the existing comparison sort. + +### Memory + +- The packed seed holds one `(u64, I)` record per selected suffix in each + active phase-1 task. On the annotated GRCh38 `u64` run this added 366-376 + MiB (about 4.1%) peak RSS. Callers must opt in so this bounded worker scratch + and any explicitly budgeted ranked-text copy are never imposed silently. + ## [v0.7.0](https://github.com/COMBINE-lab/caps-sa/releases/tag/v0.7.0) — 2026-08-13 ### Added diff --git a/README.md b/README.md index 5d9a382..45c614d 100644 --- a/README.md +++ b/README.md @@ -73,16 +73,46 @@ build_ext_mem(&text, &opts, |sa_pos| { })?; ``` +Byte-valued workloads can optionally seed each external-memory phase-1 sort +with a fixed-depth packed prefix key. Pre-encoded dense alphabets use no +text-sized copy: + +```rust +use caps_sa::{ExtMemOpts, PackedPrefixSeedPolicy}; + +let opts = ExtMemOpts::default() + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); +``` + +The seed requires unbounded comparisons and a `LimitProvider` whose +`boundary_rank()` describes how segment ends sort. `PlainText` and +`SegmentedText` provide the standard shorter-first declaration; custom +comparators such as STAR's longer-first convention declare their own. Other +symbol types, finite contexts, unsupported boundary semantics, and ineligible +alphabets fall back to the comparison sort. + +The policy is disabled by default because each active task additionally holds +one `(u64, I)` key record per selected suffix in its subarray. A gapped byte +alphabet can use `PackedPrefixSeedPolicy::remap(max_extra_bytes)`, which makes +the possible text-sized ranked copy explicit and declines it when the budget +is insufficient. + Inputs with many repeated long contexts can opt into bounded geometric LCP memoization during the final partition merges: ```rust -use caps_sa::LcpMemoizationPolicy; +use caps_sa::{LcpMemoizationPolicy, PackedPrefixSeedPolicy}; let opts = ExtMemOpts::default() + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly) .lcp_memoization(LcpMemoizationPolicy::geometric()); ``` +The two policies act on different phases and compose on the ruSTAR workload: +the current 32-core GRCh38 + GENCODE v50 A/B measured 171.205 seconds with +memoization alone and 134.618 seconds with both enabled (21.4% faster), with +identical output. + The direct path remains the default: memoization pays only when the workload contains enough repeated long contexts. See the [user guide](https://combine-lab.github.io/caps-sa/concepts/geometric-memoization/) diff --git a/bench/README.md b/bench/README.md index 6ee901f..9a5a381 100644 --- a/bench/README.md +++ b/bench/README.md @@ -73,7 +73,8 @@ cargo build --release --example rustar_segmented_bench CAPS_SA_PROFILE=1 taskset -c 0-31 \ target/release/examples/rustar_segmented_bench /path/to/fixture \ - --threads 32 --repeat 3 --work-dir /path/to/fast/temp + --threads 32 --repeat 3 --work-dir /path/to/fast/temp \ + --packed-prefix-seed ``` Omit `taskset` on platforms where it is unavailable. Use `--in-mem` only for @@ -82,6 +83,13 @@ checksums are useful regression signals, but hashes are not proofs of equality; release validation should compare emitted position streams exactly or use a direct suffix comparator on a smaller fixture. +For a packed-prefix A/B, use `--packed-prefix-seed` for the candidate and +`--no-packed-prefix-seed` for the comparison; these override the corresponding +`ExtMemOpts::from_env()` policy. Combined with `--packed-prefix-seed`, +`--no-boundary-rank` is a separate negative control that keeps the policy +enabled but makes the provider semantically ineligible, verifying the +automatic fallback path. + ## Results Machine: 64-core x86_64 Linux node, 1 socket, AVX2 enabled. diff --git a/docs/src/content/docs/concepts/algorithm.md b/docs/src/content/docs/concepts/algorithm.md index 8875e0f..37f13d0 100644 --- a/docs/src/content/docs/concepts/algorithm.md +++ b/docs/src/content/docs/concepts/algorithm.md @@ -18,7 +18,7 @@ When a fallback scan *is* needed, it runs through the **SIMD LCP fast path** (se For inputs too large for a single merge-sort pass, caps-sa wraps the kernel in a **sample sort**. The in-memory path uses the conventional four stages; the external path fuses the first three to avoid a complete intermediate spill. With `p` subproblems: 1. **Presample pivots.** Sort a small deterministic position sample and pick `p − 1` evenly-spaced pivots. These define `p` partition ranges that together cover the whole SA. -2. **Sort + distribute.** Split positions into `p` subarrays. Sort each in an outer Rayon task, binary-search its pivot splits, and write each sorted slice directly to its final partition bucket. +2. **Sort + distribute.** Split positions into `p` subarrays. Sort each in an outer Rayon task, binary-search its pivot splits, and write each sorted slice directly to its final partition bucket. Eligible byte texts may opt into a segment-bounded packed-prefix seed: keys decide short prefixes, while equal-key groups retain the complete LCP merge-sort. 3. **Per-partition merge.** Load each partition's bucket, cascade 2-way LCP-enhanced merges over its sub-slices, and emit the resulting sorted positions through the caller's closure. Because the partitions are globally ordered, emitting them in turn yields the full SA in lexicographic order, and **peak RAM stays bounded at `~O(text + n/p)` per worker** regardless of input size. @@ -27,7 +27,7 @@ Because the partitions are globally ordered, emitting them in turn yields the fu The external-memory path (`build_ext_mem`) is the default for production-scale genomes. Final partition buckets are **disk-spilling**. Positions are read back partition-by-partition only when that partition is merged, then streamed straight out—the suffix array is **never fully materialised in memory**. -The bucket pool collapses the `p` logical partition buckets onto a small set of physical temp files (one per worker by default) to keep kernel-level write contention bounded. Tuning knobs—subproblem count, working directory, physical file count, and LCP memoization policy—are on [`ExtMemOpts`](/caps-sa/reference/api/#extmemopts). +The bucket pool collapses the `p` logical partition buckets onto a small set of physical temp files (one per worker by default) to keep kernel-level write contention bounded. Tuning knobs—subproblem count, working directory, physical file count, packed-prefix seed, and LCP memoization policy—are on [`ExtMemOpts`](/caps-sa/reference/api/#extmemopts). :::note[Positioned I/O] The pooled bucket path uses positioned reads/writes (`pread`/`pwrite` on Unix, `seek_read`/`seek_write` on Windows) so many workers can share one file handle without a shared cursor. It is portable across Unix and Windows as of v0.6.1. diff --git a/docs/src/content/docs/getting-started/quick-start.md b/docs/src/content/docs/getting-started/quick-start.md index dcfda6f..2dc0414 100644 --- a/docs/src/content/docs/getting-started/quick-start.md +++ b/docs/src/content/docs/getting-started/quick-start.md @@ -36,15 +36,35 @@ build_ext_mem(&text, &opts, |sa_pos| { })?; ``` +### Optional: seed phase 1 from packed prefixes + +Dense byte alphabets can opt into segment-aware fixed-depth keys for the +external-memory phase-1 sorts: + +```rust +use caps_sa::{ExtMemOpts, PackedPrefixSeedPolicy}; + +let opts = ExtMemOpts::default() + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); +``` + +The mode is disabled by default because it adds one key record per selected +suffix in each active phase-1 task. It requires unbounded comparisons and a +`LimitProvider` with a representable `boundary_rank()`; otherwise caps-sa +falls back automatically. `DenseAlphabetOnly` never creates a text-sized copy. +See the [library API](/caps-sa/reference/api/#packed-prefix-phase-1-seed) for +gapped alphabets and custom boundary conventions. + ### Optional: reuse repeated long contexts For inputs with many long repeated contexts, opt into geometric LCP memoization for the final partition merges: ```rust -use caps_sa::{ExtMemOpts, LcpMemoizationPolicy}; +use caps_sa::{ExtMemOpts, LcpMemoizationPolicy, PackedPrefixSeedPolicy}; let opts = ExtMemOpts::default() + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly) .lcp_memoization(LcpMemoizationPolicy::geometric()); ``` @@ -56,6 +76,11 @@ repetitive inputs can be neutral or slightly slower. See [Geometric LCP memoization](/caps-sa/concepts/geometric-memoization/) before enabling it broadly. +The policies compose: packed prefixes reduce phase 1, while geometric +memoization reduces phase 4. On complete ruSTAR-shaped GRCh38 + GENCODE v50, +the memoized build improved from 171.205 to 134.618 seconds when the packed +seed was also enabled, with identical output. + ## Library: sort only a subset When many positions should be excluded from the sort (e.g. `N`s or inter-sequence spacers in a genome), hand only the positions you want sorted to a `*_for_positions` entry point. The others never enter the sort: diff --git a/docs/src/content/docs/reference/api.md b/docs/src/content/docs/reference/api.md index 8717ffb..e0a6de3 100644 --- a/docs/src/content/docs/reference/api.md +++ b/docs/src/content/docs/reference/api.md @@ -80,6 +80,7 @@ Sample-sort / external-memory tuning. `Default::default()` is tuned for genome-s | `work_dir` | `std::env::temp_dir()` | Directory for the temporary bucket files. | | `physical_file_count` | `0` → auto | Physical temp files in the bucket pool. `0` picks one per worker; the `p` logical partition buckets collapse onto this pool. | | `ordered_phase4_emit` | `false` | Opt into the bounded ordered emitter (lower transient residency on skewed inputs, slightly slower on balanced ones). | +| `packed_prefix_seed` | `Disabled` | Optional fixed-depth packed-key seed for external-memory phase 1. `DenseAlphabetOnly` adds no text-sized copy. | | `lcp_memoization` | `Disabled` | Optional exact long-LCP reuse. `LcpMemoizationPolicy::geometric()` selects the GRCh38-tuned defaults. | `ExtMemOpts` is non-exhaustive. Construct it with `default()` or `from_env()` @@ -87,6 +88,56 @@ and use its builder methods rather than an external struct literal. The `CAPS_SA_N_PHYS` environment variable overrides `physical_file_count` for one-off runs. +#### Packed-prefix phase-1 seed + +Packed-prefix seeding is an external-memory-only, opt-in acceleration for byte +texts. The safe mode expects symbols already encoded as the dense range +`0..alphabet_size`: + +```rust +use caps_sa::{ExtMemOpts, PackedPrefixSeedPolicy}; + +let opts = ExtMemOpts::default() + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); +``` + +Each selected suffix receives a segment-bounded `u64` key. The key decides +order for prefixes that differ within its fixed depth; equal-key groups retain +the complete LCP merge comparator. The mode requires: + +- symbol type exactly `u8`; +- `max_context == usize::MAX`; +- a `LimitProvider::boundary_rank()` declaration; and +- an alphabet with room for one reserved boundary code. + +`PlainText` and `SegmentedText` declare `BoundaryRank::ShorterFirst`. A custom +STAR-compatible provider whose `boundary_order` places an ending suffix above +a continuing suffix declares `BoundaryRank::LongerFirst`: + +```rust +fn boundary_rank(&self) -> Option { + Some(caps_sa::BoundaryRank::LongerFirst) +} +``` + +`boundary_rank()` describes comparator semantics; it does not activate the +optimization. Unsupported builds fall back without changing output. + +For a gapped byte alphabet, +`PackedPrefixSeedPolicy::remap(max_extra_bytes)` permits an order-preserving +dense copy only if its exact text-length allocation fits the supplied budget. +Allocation failure or an insufficient budget falls back to comparison sort. +Dense inputs avoid the copy under both policies. + +Phase 1 additionally holds one `(u64, I)` record per selected suffix in each +active subarray. For the 6.18-billion-position ruSTAR `u64` construction at +8,192 partitions and 32 workers, that is about 11.5 MiB per active worker; +measured peak RSS increased by 366–376 MiB (about 4.1%). + +`ExtMemOpts::from_env()` recognizes `CAPS_SA_PACKED_PREFIX_SEED` for the +dense-only policy and `CAPS_SA_PACKED_PREFIX_REMAP_BYTES` for an explicit +remap budget. The latter takes precedence when both are set. + #### Geometric LCP memoization Memoization is opt-in and applies only to phase-4 partition merges: diff --git a/docs/src/content/docs/reference/performance.md b/docs/src/content/docs/reference/performance.md index 2bc7f1e..69815c3 100644 --- a/docs/src/content/docs/reference/performance.md +++ b/docs/src/content/docs/reference/performance.md @@ -1,6 +1,6 @@ --- title: Performance -description: Production-shaped caps-sa 0.7 measurements and the standard upstream comparison. +description: Production-shaped caps-sa measurements and the standard upstream comparison. --- Unless stated otherwise, numbers are suffix-array construction time and output @@ -35,6 +35,34 @@ On the focused chromosome-21 backbone plus every GENCODE-derived junction flank, stable measurements improved from 11.97–12.00 seconds to 7.77–7.80 seconds. All 359,616,038 emitted positions matched the reference. +## Optional packed-prefix phase-1 seed + +The opt-in packed-prefix seed was measured against the released v0.7.0 +baseline on the same complete fixture, with one warm-up followed by three +interleaved measured runs. Values below are medians: + +| Configuration | Build | User CPU | Peak RSS | Phase 1 | Phase 4 | +| --- | ---: | ---: | ---: | ---: | ---: | +| Direct LCP | 198.924 s | 5,153.25 s | 9,142,300 KiB | 49.038 s | 144.923 s | +| Geometric memo only | 171.205 s | 4,684.95 s | 9,161,688 KiB | 49.034 s | 117.021 s | +| Packed seed only | 162.650 s | 3,987.38 s | 9,527,240 KiB | 12.204 s | 145.187 s | +| Packed seed + geometric memo | **134.618 s** | **3,507.54 s** | 9,536,692 KiB | **12.179 s** | **117.265 s** | + +The seed reduces phase 1 by 75.1% and improves the memoized ruSTAR +configuration by 21.4%. Memoization changes phase 4 by 19.25% without the seed +and 19.23% with it, confirming that the policies compose rather than competing +for the same work. Every run emitted 6,176,694,310 positions with ordered hash +`e81c8f9881e322148741a23c92ae2000`. + +Peak RSS increased by 366–376 MiB (about 4.1%), matching the bounded +`(u64, u64)` key records held by the 32 active phase-1 tasks. The dense ruSTAR +alphabet required no ranked-text copy. + +On chr21 without annotations, the direct build improved from a 1.291-second +median to 0.934 seconds (27.7%). On the chr21 backbone plus every annotation- +derived flank it improved from 7.709 to 6.025 seconds (21.8%), and the latter +matched all 359,616,038 reference positions exactly. + ## Geometric memoization The table above has geometric memoization enabled on both sides; do not add its diff --git a/examples/rustar_segmented_bench.rs b/examples/rustar_segmented_bench.rs index 3a55ff2..6e7d3bb 100644 --- a/examples/rustar_segmented_bench.rs +++ b/examples/rustar_segmented_bench.rs @@ -20,7 +20,8 @@ //! //! ```text //! cargo run --release --example rustar_segmented_bench -- FIXTURE_DIR \ -//! [--threads N] [--repeat N] [--in-mem] [--work-dir DIR] +//! [--threads N] [--repeat N] [--in-mem] [--work-dir DIR] \ +//! [--packed-prefix-seed | --no-packed-prefix-seed] [--no-boundary-rank] //! ``` //! //! Reports the wall time of the caps-sa build alone, plus the emitted @@ -37,8 +38,8 @@ use std::process; use std::time::Instant; use caps_sa::{ - ExtMemOpts, LimitProvider, Opts, SegmentedText, build_ext_mem_for_filter_with, - build_in_memory_for_positions_with, + ExtMemOpts, LimitProvider, Opts, PackedPrefixSeedPolicy, SegmentedText, + build_ext_mem_for_filter_with, build_in_memory_for_positions_with, }; /// rustar-aligner's `StarSegmentedText`: `SegmentedText` limits with @@ -46,6 +47,7 @@ use caps_sa::{ /// ascending position on ties). struct StarSegmentedText { inner: SegmentedText, + boundary_rank: Option, } impl LimitProvider for StarSegmentedText { @@ -58,6 +60,13 @@ impl LimitProvider for StarSegmentedText { fn boundary_order(&self, p_a: usize, lim_a: usize, p_b: usize, lim_b: usize) -> Ordering { lim_b.cmp(&lim_a).then(p_a.cmp(&p_b)) } + + /// Declare that this comparator's boundary convention is representable by + /// packed keys. Activation remains a separate `ExtMemOpts` policy below. + #[inline] + fn boundary_rank(&self) -> Option { + self.boundary_rank + } } struct Args { @@ -66,10 +75,14 @@ struct Args { repeat: usize, in_mem: bool, work_dir: Option, + packed_prefix_seed: Option, + no_boundary_rank: bool, } const USAGE: &str = "usage: rustar_segmented_bench FIXTURE_DIR [--threads N] \ - [--repeat N] [--in-mem] [--work-dir DIR]"; + [--repeat N] [--in-mem] [--work-dir DIR] \ + [--packed-prefix-seed | --no-packed-prefix-seed] \ + [--no-boundary-rank]"; fn usage_error(message: &str) -> ! { eprintln!("error: {message}\n{USAGE}"); @@ -89,6 +102,8 @@ fn parse_args() -> Args { let mut repeat = 1usize; let mut in_mem = false; let mut work_dir = None; + let mut packed_prefix_seed = None; + let mut no_boundary_rank = false; let mut i = 1; while i < argv.len() { match argv[i].as_str() { @@ -119,6 +134,24 @@ fn parse_args() -> Args { work_dir = Some(PathBuf::from(option_value(&argv, i, "--work-dir"))); i += 2; } + "--packed-prefix-seed" => { + if packed_prefix_seed.is_some() { + usage_error("packed-prefix policy may be specified only once"); + } + packed_prefix_seed = Some(true); + i += 1; + } + "--no-packed-prefix-seed" => { + if packed_prefix_seed.is_some() { + usage_error("packed-prefix policy may be specified only once"); + } + packed_prefix_seed = Some(false); + i += 1; + } + "--no-boundary-rank" => { + no_boundary_rank = true; + i += 1; + } "--help" | "-h" => { eprintln!("{USAGE}"); process::exit(0); @@ -135,12 +168,17 @@ fn parse_args() -> Args { if positional.len() != 1 { usage_error("expected exactly one fixture directory"); } + if in_mem && (packed_prefix_seed.is_some() || no_boundary_rank) { + usage_error("packed-prefix controls apply only to the external-memory path"); + } Args { fixture: PathBuf::from(&positional[0]), threads, repeat, in_mem, work_dir, + packed_prefix_seed, + no_boundary_rank, } } @@ -191,6 +229,7 @@ fn main() { let lp = StarSegmentedText { inner: SegmentedText::from_ends(n, ends), + boundary_rank: (!args.no_boundary_rank).then_some(caps_sa::BoundaryRank::LongerFirst), }; for round in 0..args.repeat { @@ -216,6 +255,14 @@ fn main() { } } else { let mut opts = ExtMemOpts::from_env(); + if let Some(enabled) = args.packed_prefix_seed { + let policy = if enabled { + PackedPrefixSeedPolicy::DenseAlphabetOnly + } else { + PackedPrefixSeedPolicy::Disabled + }; + opts = opts.packed_prefix_seed(policy); + } if let Some(dir) = &args.work_dir { opts = opts.work_dir(dir); } diff --git a/src/ext_mem.rs b/src/ext_mem.rs index 2a4ebbb..087ad6f 100644 --- a/src/ext_mem.rs +++ b/src/ext_mem.rs @@ -35,7 +35,8 @@ use crate::lcp::{LcpDispatch, Symbol}; use crate::lcp_memo::{ GeometricMemo, GeometricMemoizationConfig, LcpMemoizationPolicy, MemoConfig, MemoStats, }; -use crate::limits::{LimitProvider, PlainText}; +use crate::limits::{BoundaryRank, LimitProvider, PlainText}; +use crate::pack::{PackedPrefixSeedPolicy, Packer}; use crate::sample_sort; /// Emit a phase-timing line to stderr if `CAPS_SA_PROFILE` is set in @@ -86,6 +87,11 @@ pub struct ExtMemOpts { /// GRCh38 32-thread benchmark because of channel coordination and /// backpressure, so it is opt-in. pub ordered_phase4_emit: bool, + /// Policy for seeding external-memory phase 1 with fixed-depth packed + /// prefix keys. Disabled by default; pre-encoded dense byte texts can opt + /// into [`PackedPrefixSeedPolicy::DenseAlphabetOnly`] without allocating a + /// second text-sized buffer. + pub packed_prefix_seed: PackedPrefixSeedPolicy, /// Policy for reusing exact long-LCP intervals during phase-4 partition /// merges. Disabled by default; callers with long repeated contexts can /// opt into [`LcpMemoizationPolicy::Geometric`]. @@ -104,6 +110,7 @@ impl Default for ExtMemOpts { work_dir: std::env::temp_dir(), physical_file_count: 0, ordered_phase4_emit: false, + packed_prefix_seed: PackedPrefixSeedPolicy::Disabled, lcp_memoization: LcpMemoizationPolicy::Disabled, collect_lcp_memoization_stats: false, } @@ -128,6 +135,10 @@ impl ExtMemOpts { /// - `CAPS_SA_N_PHYS`: physical backing-file count /// - `CAPS_SA_MAX_CONTEXT`: LCP comparison cap /// - `CAPS_SA_ORDERED_PHASE4=1|true|yes|on`: bounded ordered phase-4 emit + /// - `CAPS_SA_PACKED_PREFIX_SEED=1|true|yes|on`: packed phase-1 keys for a + /// dense byte alphabet + /// - `CAPS_SA_PACKED_PREFIX_REMAP_BYTES`: packed phase-1 keys with this + /// maximum ranked-text allocation /// - `CAPS_SA_GEOMETRIC_MEMO=1|true|yes|on`: geometric LCP memoization /// - `CAPS_SA_MEMO_PROBE`: ordinary symbols compared before table lookup /// - `CAPS_SA_MEMO_MIN_LCP`: minimum exact LCP admitted to the table @@ -135,9 +146,11 @@ impl ExtMemOpts { /// - `CAPS_SA_MEMO_CAPACITY`: maximum entries per partition table /// - `CAPS_SA_MEMO_STATS=1|true|yes|on`: detailed memoization counters /// - /// Invalid and zero-valued numeric overrides are ignored, preserving the - /// corresponding defaults. Memoization tuning variables take effect only - /// when `CAPS_SA_GEOMETRIC_MEMO` enables the policy. + /// Invalid numeric overrides are ignored. Zero-valued count and + /// memoization overrides preserve their defaults; a zero remap budget is + /// valid and permits only already-dense inputs. Memoization tuning + /// variables take effect only when `CAPS_SA_GEOMETRIC_MEMO` enables the + /// policy. pub fn from_env() -> Self { let mut opts = Self::default(); if let Some(dir) = @@ -157,6 +170,11 @@ impl ExtMemOpts { if read_env_bool("CAPS_SA_ORDERED_PHASE4") { opts.ordered_phase4_emit = true; } + if let Some(max_extra_bytes) = read_env_usize("CAPS_SA_PACKED_PREFIX_REMAP_BYTES") { + opts.packed_prefix_seed = PackedPrefixSeedPolicy::remap(max_extra_bytes); + } else if read_env_bool("CAPS_SA_PACKED_PREFIX_SEED") { + opts.packed_prefix_seed = PackedPrefixSeedPolicy::DenseAlphabetOnly; + } if read_env_bool("CAPS_SA_GEOMETRIC_MEMO") { let mut config = GeometricMemoizationConfig::default(); if let Some(v) = read_env_nonzero_usize("CAPS_SA_MEMO_PROBE") { @@ -207,6 +225,12 @@ impl ExtMemOpts { self } + /// Builder-style setter for [`Self::packed_prefix_seed`]. + pub fn packed_prefix_seed(mut self, packed_prefix_seed: PackedPrefixSeedPolicy) -> Self { + self.packed_prefix_seed = packed_prefix_seed; + self + } + /// Builder-style setter for [`Self::lcp_memoization`]. pub fn lcp_memoization(mut self, lcp_memoization: impl Into) -> Self { self.lcp_memoization = lcp_memoization.into(); @@ -1467,6 +1491,27 @@ where (1..p).map(|i| sample[(i * m / p).min(m - 1)]).collect() } +/// Prepare the optional packed seed only after all cheap eligibility checks. +fn packed_prefix_seed_params( + text: &[S], + lp: &L, + opts: &ExtMemOpts, +) -> Option<(Packer, BoundaryRank)> { + if opts.packed_prefix_seed == PackedPrefixSeedPolicy::Disabled + || opts.max_context != usize::MAX + || std::any::TypeId::of::() != std::any::TypeId::of::() + { + return None; + } + + // `boundary_rank` is a semantic capability declaration, not the enable + // switch. Ask only after the caller has opted in and all cheap generic + // eligibility checks pass; then scan the alphabet at most once. + let rank = lp.boundary_rank()?; + let packer = crate::pack::seed_params(text, opts.packed_prefix_seed)?; + Some((packer, rank)) +} + /// Sort each phase-1 subarray and distribute its sorted pieces directly to /// the final partition buckets. /// @@ -1497,6 +1542,14 @@ where let chunk_size = n.div_ceil(p); let partition_buckets: Vec> = (0..p).map(|j| Mutex::new(mk_bucket(j))).collect(); let task_local_sort = p >= rayon::current_num_threads().max(1); + let packed_seed = packed_prefix_seed_params(text, lp, opts); + if opts.packed_prefix_seed != PackedPrefixSeedPolicy::Disabled { + profile_log(if packed_seed.is_some() { + "phase1 packed-prefix seed active" + } else { + "phase1 packed-prefix seed unavailable; using comparison sort" + }); + } (0..p).into_par_iter().try_for_each(|i| -> io::Result<()> { let start = (i * chunk_size).min(n); @@ -1511,7 +1564,22 @@ where let mut sa_w = vec![I::zero(); len]; let mut lcp_arr = vec![I::zero(); len]; let mut lcp_w = vec![I::zero(); len]; - if task_local_sort { + if let Some((packer, rank)) = packed_seed.as_ref() { + crate::pack::seed_subarray( + text, + lp, + packer, + *rank, + &mut sa, + &mut lcp_arr, + &mut sa_w, + &mut lcp_w, + dispatch, + task_local_sort, + ); + // Sorted by key, with the merge kernel run only inside equal-key + // groups. + } else if task_local_sort { sample_sort::merge_sort_task_local( text, lp, @@ -2346,6 +2414,11 @@ mod tests { // all touched variables before releasing it. unsafe { std::env::set_var(key, value) }; } + + fn remove(&self, key: &'static str) { + // Serialized and restored under ENV_LOCK, as in `set`. + unsafe { std::env::remove_var(key) }; + } } impl Drop for EnvGuard { @@ -2406,6 +2479,7 @@ mod tests { #[test] fn memoization_policy_defaults_to_disabled() { let opts = ExtMemOpts::default(); + assert_eq!(opts.packed_prefix_seed, PackedPrefixSeedPolicy::Disabled); assert_eq!(opts.lcp_memoization, LcpMemoizationPolicy::Disabled); assert!(!opts.collect_lcp_memoization_stats); @@ -2420,6 +2494,105 @@ mod tests { ); } + #[test] + fn packed_seed_policy_matches_direct_external_output() { + let dense: Vec = (0..4096).map(|i| ((i * 17 + i / 13) % 6) as u8).collect(); + let direct = ext_mem_sa(&dense, 17); + + let dir = tempdir().unwrap(); + let opts = ExtMemOpts::default() + .subproblem_count(17) + .physical_file_count(1) + .work_dir(dir.path()) + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); + let mut seeded = Vec::with_capacity(dense.len()); + build_ext_mem(&dense, &opts, |pos| { + seeded.push(pos); + Ok(()) + }) + .unwrap(); + assert_eq!(seeded, direct); + + let gapped: Vec = dense.iter().map(|&b| b.saturating_mul(17)).collect(); + let direct = ext_mem_sa(&gapped, 17); + let dir = tempdir().unwrap(); + let opts = ExtMemOpts::default() + .subproblem_count(17) + .physical_file_count(1) + .work_dir(dir.path()) + .packed_prefix_seed(PackedPrefixSeedPolicy::remap(gapped.len())); + let mut seeded = Vec::with_capacity(gapped.len()); + build_ext_mem(&gapped, &opts, |pos| { + seeded.push(pos); + Ok(()) + }) + .unwrap(); + assert_eq!(seeded, direct); + } + + #[test] + fn packed_seed_eligibility_precedes_provider_query() { + struct UnexpectedRank(usize); + impl LimitProvider for UnexpectedRank { + fn lim_at(&self, p: usize) -> usize { + self.0 - p + } + + fn boundary_rank(&self) -> Option { + panic!("ineligible seed must not query boundary_rank") + } + } + + let text = [0u8, 1, 2, 3]; + let lp = UnexpectedRank(text.len()); + assert!(packed_prefix_seed_params(&text, &lp, &ExtMemOpts::default()).is_none()); + + let finite = ExtMemOpts::default() + .max_context(3) + .packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); + assert!(packed_prefix_seed_params(&text, &lp, &finite).is_none()); + + let signed = [0i8, 1, 2, 3]; + let enabled = + ExtMemOpts::default().packed_prefix_seed(PackedPrefixSeedPolicy::DenseAlphabetOnly); + assert!(packed_prefix_seed_params(&signed, &lp, &enabled).is_none()); + + struct NoRank(usize); + impl LimitProvider for NoRank { + fn lim_at(&self, p: usize) -> usize { + self.0 - p + } + } + assert!(packed_prefix_seed_params(&text, &NoRank(text.len()), &enabled).is_none()); + } + + #[test] + fn from_env_parses_packed_seed_policy() { + let _lock = ENV_LOCK.lock().unwrap(); + let keys = [ + "CAPS_SA_PACKED_PREFIX_SEED", + "CAPS_SA_PACKED_PREFIX_REMAP_BYTES", + ]; + let env = EnvGuard::capture(&keys); + for &key in &keys { + env.remove(key); + } + + env.set("CAPS_SA_PACKED_PREFIX_SEED", "true"); + assert_eq!( + ExtMemOpts::from_env().packed_prefix_seed, + PackedPrefixSeedPolicy::DenseAlphabetOnly + ); + + env.set("CAPS_SA_PACKED_PREFIX_REMAP_BYTES", "12345"); + assert_eq!( + ExtMemOpts::from_env().packed_prefix_seed, + PackedPrefixSeedPolicy::Remap { + max_extra_bytes: 12_345, + } + ); + } + #[test] fn geometric_policy_matches_direct_output() { let mut text = Vec::new(); diff --git a/src/lib.rs b/src/lib.rs index 06fc661..b862839 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -22,6 +22,7 @@ mod ext_mem; mod lcp; mod lcp_memo; mod limits; +mod pack; mod sample_sort; pub use ext_mem::{ @@ -36,7 +37,8 @@ pub use ext_mem::{ }; pub use lcp::{LcpDispatch, Symbol, lcp, lcp_scalar, lcp_u8, suffix_cmp}; pub use lcp_memo::{GeometricMemoizationConfig, LcpMemoizationPolicy}; -pub use limits::{LimitProvider, PlainText, SegmentedText}; +pub use limits::{BoundaryRank, LimitProvider, PlainText, SegmentedText}; +pub use pack::PackedPrefixSeedPolicy; pub use sample_sort::{ Opts, build_in_memory, build_in_memory_for_positions, build_in_memory_for_positions_with, build_in_memory_for_positions_with_opts, build_in_memory_with, build_in_memory_with_opts, diff --git a/src/limits.rs b/src/limits.rs index e4e5ead..2faa318 100644 --- a/src/limits.rs +++ b/src/limits.rs @@ -73,6 +73,53 @@ pub trait LimitProvider: Sync { let _ = (p_a, p_b); lim_a.cmp(&lim_b) } + + /// Which side of a real symbol this provider's boundary convention puts + /// the end of a suffix on, when that convention is expressible. + /// + /// Answering `Some` makes the provider eligible for the optional + /// [`PackedPrefixSeedPolicy`][crate::PackedPrefixSeedPolicy]. Activation + /// remains a separate [`ExtMemOpts`][crate::ExtMemOpts] choice. The key + /// packs `min(k, lim_at(p))` symbols and pads the rest with a reserved + /// sentinel placed according to this answer, so key order agrees with + /// [`boundary_order`][Self::boundary_order] whenever the key decides at + /// all. Keys that tie still defer to `boundary_order` itself, which is + /// what lets a convention with a position tie-break (STAR's + /// `lim_b.cmp(&lim_a).then(p_a.cmp(&p_b))`) work: the key never has to + /// express the tie-break, only never to contradict it. + /// + /// The contract is exactly this: for suffixes `a` and `b` whose shared + /// prefix ends because one of them reached its limit, + /// `boundary_order(a, .., b, ..)` must be `Less` iff `a` is the one that + /// ended, under [`BoundaryRank::ShorterFirst`], and `Greater` iff `a` is + /// the one that ended, under [`BoundaryRank::LongerFirst`]. + /// + /// The default is `None`, which declines the seed even when a caller opts + /// in. Answer it only if your `boundary_order` decides purely by which + /// suffix ended first, with at most a tie-break between suffixes that end + /// at the same offset. + #[inline] + fn boundary_rank(&self) -> Option { + None + } +} + +/// How a provider's [`boundary_order`][LimitProvider::boundary_order] ranks a +/// suffix that ends at its segment boundary against one that keeps going. +/// +/// This is the one fact a fixed-depth packed key needs in order to represent +/// a segmented comparator: a key pads a short suffix out to full width, and +/// the padding symbol has to fall on the correct side of every real symbol. +#[derive(Copy, Clone, Debug, PartialEq, Eq)] +pub enum BoundaryRank { + /// The suffix that ends first is smaller, the standard generalized-SA + /// convention and the default `boundary_order`. Keys pad with a symbol + /// below every real one. + ShorterFirst, + /// The suffix that ends first is *larger*, equivalently the longer one is + /// smaller. STAR's spacer-as-largest ordering. Keys pad with a symbol + /// above every real one. + LongerFirst, } /// Default provider for non-segmented texts: `lim_at(p) = n - p`. @@ -99,6 +146,13 @@ impl LimitProvider for PlainText { fn lim_at(&self, p: usize) -> usize { self.n - p } + + /// End of text behaves like a boundary the shorter suffix hits first, + /// which is the lexicographic convention this provider implements. + #[inline] + fn boundary_rank(&self) -> Option { + Some(BoundaryRank::ShorterFirst) + } } /// Provider for texts partitioned into segments at known cumulative @@ -272,6 +326,13 @@ impl LimitProvider for SegmentedText { self.n - p } } + + /// `SegmentedText` keeps the default `boundary_order`, which is + /// shorter-is-smaller. + #[inline] + fn boundary_rank(&self) -> Option { + Some(BoundaryRank::ShorterFirst) + } } #[cfg(test)] diff --git a/src/pack.rs b/src/pack.rs new file mode 100644 index 0000000..0f5710d --- /dev/null +++ b/src/pack.rs @@ -0,0 +1,585 @@ +//! Packed fixed-depth keys, and the phase-1 subarray seed built on them. +//! +//! Phase 1 sorts each subarray from singletons, which is the case the merge +//! kernel handles worst: every leaf merge starts at `m = 0` and orders two +//! suffixes by scanning the text at two random addresses. Sorting by a packed +//! key first resolves the leading `k` symbols with no text access at all (16 +//! symbols for a 6-letter genome alphabet at 4 bits each), and yields the LCP +//! between adjacent runs for free from `(key_a ^ key_b).leading_zeros()`. Only +//! suffixes agreeing through the whole key reach the merge kernel, on the +//! short slice they occupy. +//! +//! This is the bounded-memory counterpart to prefix doubling. Doubling itself +//! is not available here: it needs a rank for every position in the text, +//! which is exactly the memory the external-memory path refuses to spend. A +//! fixed-depth key needs one `u64` per record of the subarray in flight. +//! +//! The key is segment-aware. It packs `min(k, lim_at(p))` symbols, so it never +//! reads into the next segment, and pads with a reserved sentinel placed on +//! the side the provider's [`BoundaryRank`] demands. That is what lets a +//! splice-junction index, where the segment boundary rather than the text end +//! terminates most suffixes, take the path at all. + +use crate::Index; +use crate::lcp::{LcpDispatch, Symbol}; +use crate::limits::{BoundaryRank, LimitProvider}; +use crate::sample_sort; +use rayon::prelude::*; + +/// Whether the external-memory phase-1 sort may start from fixed-depth packed +/// prefix keys. +/// +/// The seed is deliberately opt-in. It adds one `(u64, I)` key record per +/// selected suffix in every active phase-1 task, and it is useful only when +/// the text is byte-valued, comparisons are unbounded, and the +/// [`LimitProvider`] declares a representable [`BoundaryRank`]. Ineligible +/// builds fall back to the ordinary LCP merge-sort without changing output. +#[non_exhaustive] +#[derive(Copy, Clone, Debug, Default, PartialEq, Eq)] +pub enum PackedPrefixSeedPolicy { + /// Use the ordinary comparison-based phase-1 sort. + #[default] + Disabled, + /// Use packed keys only when the byte values occurring in the text are + /// already the dense range `0..alphabet_size`. + /// + /// This mode never allocates a second text-sized buffer. It is the + /// recommended mode for pre-encoded genomic text such as ruSTAR's + /// `A=0, C=1, G=2, T=3, N=4, spacer=5` representation. + DenseAlphabetOnly, + /// Permit an order-preserving dense copy for a gapped byte alphabet, but + /// only when its exact text-length allocation fits `max_extra_bytes`. + /// + /// Dense inputs still avoid the copy. If the budget is too small or the + /// allocation cannot be reserved, construction falls back to the ordinary + /// phase-1 sort. + Remap { max_extra_bytes: usize }, +} + +impl PackedPrefixSeedPolicy { + /// Opt into the allocation-free dense-alphabet path. + pub const fn dense_alphabet_only() -> Self { + Self::DenseAlphabetOnly + } + + /// Permit a dense ranked-text copy up to the given allocation budget. + pub const fn remap(max_extra_bytes: usize) -> Self { + Self::Remap { max_extra_bytes } + } +} + +/// An order-preserving remap of the bytes that actually occur in a text onto +/// a dense code range, plus the resulting key geometry. +/// +/// The field width is driven by how many *distinct* symbols a text uses, not +/// by the largest byte value in it, and the difference is not academic. A raw +/// FASTA uses six symbols, but the largest is `'T'` (84), so packing raw bytes +/// forces 8-bit fields and fits only 8 symbols per key. Ranking those six +/// bytes to `0..6` gives 4-bit fields and 16 symbols per key. A rustar-shaped +/// text is already dense (`0..=5`), so it pays nothing for the map. +/// +/// The map is monotone by construction, since codes are assigned in ascending +/// byte order. That is what keeps a packed key order-preserving: `key_a < +/// key_b` still implies `suffix_a < suffix_b`. +pub(crate) struct Packer { + /// The text with every byte replaced by its code, when the identity map + /// does not already do that. Materializing it once removes a dependent + /// table load from the packing loop. `None` when the text is already + /// dense, so the common pre-coded genomic input pays no extra memory. + ranked: Option>, + /// Bits per packed field. + bits: u32, + /// Symbols per `u64` key. + k: usize, + /// Number of distinct codes in use. Codes are `0..alphabet`; `alphabet` + /// itself is the boundary sentinel when it fits the field. + alphabet: u32, +} + +impl Packer { + /// Build the map for `text`. The field is sized to hold `alphabet`, not + /// `alphabet - 1`, because every key this module builds reserves one code + /// for the boundary sentinel. + fn new(text: &[u8], policy: PackedPrefixSeedPolicy) -> Option { + // Which bytes occur? One parallel pass, folded into a 256-entry set. + let present = text + .par_chunks(1 << 16) + .map(|c| { + let mut seen = [false; 256]; + for &b in c { + seen[b as usize] = true; + } + seen + }) + .reduce( + || [false; 256], + |mut a, b| { + for i in 0..256 { + a[i] |= b[i]; + } + a + }, + ); + + let mut code = [0u8; 256]; + let mut next = 0u16; + let mut identity = true; + for (b, &seen) in present.iter().enumerate() { + if seen { + code[b] = next as u8; + identity &= next as usize == b; + next += 1; + } + } + let bits: u32 = match next { + 0..=1 => 1, + 2..=3 => 2, + 4..=15 => 4, + _ => 8, + }; + let ranked = if identity { + None + } else { + let PackedPrefixSeedPolicy::Remap { max_extra_bytes } = policy else { + return None; + }; + if text.len() > max_extra_bytes { + return None; + } + let mut out = Vec::new(); + out.try_reserve_exact(text.len()).ok()?; + out.resize(text.len(), 0u8); + out.par_chunks_mut(1 << 16) + .zip(text.par_chunks(1 << 16)) + .for_each(|(dst, src)| { + for (d, &s) in dst.iter_mut().zip(src) { + *d = code[s as usize]; + } + }); + Some(out) + }; + Some(Self { + ranked, + bits, + k: 64 / bits as usize, + alphabet: next as u32, + }) + } + + /// Whether a boundary sentinel fits alongside the alphabet in one field. + /// + /// With 8-bit fields and 256 distinct symbols there is no spare code, so + /// keys are unavailable and the caller must fall back. + #[inline] + fn has_sentinel(&self) -> bool { + (self.alphabet as u64) < (1u64 << self.bits) + } + + /// Pack the `min(k, lim)` symbols at `text[p..]`, padding the rest with a + /// boundary sentinel placed according to `rank`. + /// + /// Never reads past `p + lim`, so a key cannot see into the next segment. + /// Under `ShorterFirst` the sentinel is code `0` and every real code is + /// shifted up by one, so a padded field is strictly below any real symbol. + /// Under `LongerFirst` the sentinel is `alphabet`, strictly above every + /// real code, and no shift is needed. + #[inline] + fn key_at_bounded(&self, text: &[u8], p: usize, lim: usize, rank: BoundaryRank) -> u64 { + debug_assert!(self.has_sentinel()); + let src = self.ranked.as_deref().unwrap_or(text); + let take = self.k.min(lim).min(src.len() - p); + let (bias, pad) = match rank { + BoundaryRank::ShorterFirst => (1u64, 0u64), + BoundaryRank::LongerFirst => (0u64, self.alphabet as u64), + }; + let mut key = 0u64; + for &c in &src[p..p + take] { + key = (key << self.bits) | (c as u64 + bias); + } + for _ in take..self.k { + key = (key << self.bits) | pad; + } + key + } + + /// Symbols the two keys share, from their first differing field. + #[inline] + fn shared_fields(&self, key_a: u64, key_b: u64) -> usize { + ((key_a ^ key_b).leading_zeros() / self.bits) as usize + } +} + +/// The alphabet map for `text`, or `None` when a packed key cannot represent +/// this text's order. +/// +/// Computed once per build and handed to [`seed_subarray`], which would +/// otherwise re-scan the whole text for every subarray. +pub(crate) fn seed_params(text: &[S], policy: PackedPrefixSeedPolicy) -> Option { + if policy == PackedPrefixSeedPolicy::Disabled { + return None; + } + // Exactly `u8`, not merely one byte wide. `Symbol` is implemented for + // `i8` too, and a packed key orders its fields as unsigned: `-1` has byte + // `0xFF` and would sort above `1`, inverting the text's real order. + if std::any::TypeId::of::() != std::any::TypeId::of::() { + return None; + } + // SAFETY: `S` is `u8` (just checked by `TypeId`, and `Symbol: 'static` so + // the comparison is exact), so a byte view over the same memory is valid + // for reads of the same length. + let bytes: &[u8] = + unsafe { std::slice::from_raw_parts(text.as_ptr() as *const u8, text.len()) }; + let packer = Packer::new(bytes, policy)?; + packer.has_sentinel().then_some(packer) +} + +/// Sort `sa` into suffix order and fill `lcp`, using a packed fixed-depth key +/// so that most of the ordering costs no text access at all. +/// +/// `sa_w` and `lcp_w` are the caller's existing merge scratch buffers. +#[allow(clippy::too_many_arguments)] +pub(crate) fn seed_subarray( + text: &[S], + lp: &L, + packer: &Packer, + rank: BoundaryRank, + sa: &mut [I], + lcp: &mut [I], + sa_w: &mut [I], + lcp_w: &mut [I], + dispatch: LcpDispatch, + task_local: bool, +) { + let len = sa.len(); + if len < 2 { + if len == 1 { + lcp[0] = I::zero(); + } + return; + } + // SAFETY: a `Packer` is prepared only when `S` is exactly `u8`. + let bytes: &[u8] = + unsafe { std::slice::from_raw_parts(text.as_ptr() as *const u8, text.len()) }; + + let mut keyed: Vec<(u64, I)> = sa + .iter() + .map(|&p| { + let pu = p.to_usize(); + (packer.key_at_bounded(bytes, pu, lp.lim_at(pu), rank), p) + }) + .collect(); + keyed.sort_unstable_by_key(|e| e.0); + for (slot, e) in sa.iter_mut().zip(keyed.iter()) { + *slot = e.1; + } + + // Walk equal-key runs. Between runs the LCP falls straight out of the key + // difference; inside one it needs the merge kernel. + let mut i = 0usize; + while i < len { + let mut j = i + 1; + while j < len && keyed[j].0 == keyed[i].0 { + j += 1; + } + if j - i > 1 { + // Mirror the caller's nesting choice. Phase 1 runs one task per + // subarray once `p` reaches the worker count, and a rayon-splitting + // merge inside such a task only adds scheduling on top of + // parallelism that already saturates. + let sort = if task_local { + sample_sort::merge_sort_task_local + } else { + sample_sort::merge_sort + }; + sort( + text, + lp, + &mut sa[i..j], + &mut sa_w[i..j], + &mut lcp[i..j], + &mut lcp_w[i..j], + usize::MAX, + dispatch, + ); + } else { + lcp[i] = I::zero(); + } + // Boundary entry: LCP against the last element of the previous run. + if i > 0 { + let a = sa[i - 1].to_usize(); + let b = sa[i].to_usize(); + let shared = packer.shared_fields(keyed[i - 1].0, keyed[i].0); + // Cap by both suffixes' limits: a sentinel field can agree with a + // real symbol's field past the end of the shorter suffix, so the + // raw count can overstate the true LCP, and a wrong LCP would + // silently corrupt the order at the next merge level. + lcp[i] = I::from_usize(shared.min(lp.lim_at(a)).min(lp.lim_at(b))); + } + i = j; + } + lcp[0] = I::zero(); +} + +#[cfg(test)] +mod tests { + use super::*; + use crate::limits::{PlainText, SegmentedText}; + + /// STAR's convention: the suffix that hits its boundary first is larger. + struct StarSegmented { + inner: SegmentedText, + } + + impl LimitProvider for StarSegmented { + fn lim_at(&self, p: usize) -> usize { + self.inner.lim_at(p) + } + + fn boundary_order( + &self, + p_a: usize, + lim_a: usize, + p_b: usize, + lim_b: usize, + ) -> std::cmp::Ordering { + lim_b.cmp(&lim_a).then(p_a.cmp(&p_b)) + } + + fn boundary_rank(&self) -> Option { + Some(BoundaryRank::LongerFirst) + } + } + + fn lcg(state: &mut u64) -> u64 { + *state = state + .wrapping_mul(6364136223846793005) + .wrapping_add(1442695040888963407); + *state >> 33 + } + + fn naive_lcp(text: &[u8], p: usize, q: usize, cap: usize) -> usize { + (0..cap) + .position(|i| text[p + i] != text[q + i]) + .unwrap_or(cap) + } + + /// Sort `positions` with the seed, then check the result is a permutation + /// in non-decreasing suffix order with an exact adjacent-LCP array under + /// the provider's own comparator. + /// + /// The assertion is the *property*, not equality with a canonical answer: + /// `SegmentedText`'s default `boundary_order` returns `Equal` for suffixes + /// that end together with equal content, so their relative order is + /// genuinely free and a stable-sort oracle is not a valid reference. + fn check_sorted_with_policy( + text: &[u8], + lp: &L, + positions: Vec, + policy: PackedPrefixSeedPolicy, + ) { + let packer = seed_params(text, policy); + assert!(packer.is_some(), "packer should be available for this text"); + let rank = lp + .boundary_rank() + .expect("test provider should declare a boundary rank"); + let len = positions.len(); + let mut sa = positions.clone(); + let mut lcp = vec![0u32; len]; + let mut sa_w = vec![0u32; len]; + let mut lcp_w = vec![0u32; len]; + let dispatch = LcpDispatch::detect(); + seed_subarray( + text, + lp, + packer.as_ref().unwrap(), + rank, + &mut sa, + &mut lcp, + &mut sa_w, + &mut lcp_w, + dispatch, + true, + ); + + let mut seen = sa.clone(); + seen.sort_unstable(); + let mut want = positions; + want.sort_unstable(); + assert_eq!(seen, want, "seed must permute its input"); + + assert_eq!(lcp[0], 0); + for (i, w) in sa.windows(2).enumerate() { + let (a, b) = (w[0] as usize, w[1] as usize); + assert!( + dispatch.suffix_cmp_with(text, lp, a, b, usize::MAX).is_le(), + "adjacent pair out of order: {a} then {b}", + ); + let cap = lp.lim_at(a).min(lp.lim_at(b)); + let want_lcp = naive_lcp(text, a, b, cap) as u32; + assert_eq!(lcp[i + 1], want_lcp, "wrong adjacent LCP for {a} then {b}",); + } + } + + fn check_sorted(text: &[u8], lp: &L, positions: Vec) { + check_sorted_with_policy( + text, + lp, + positions, + PackedPrefixSeedPolicy::DenseAlphabetOnly, + ); + } + + #[test] + fn seed_orders_a_plain_text() { + let mut state = 12345u64; + let text: Vec = (0..4096).map(|_| (lcg(&mut state) % 4) as u8).collect(); + let lp = PlainText::new(text.len()); + check_sorted(&text, &lp, (0..text.len() as u32).collect()); + } + + #[test] + fn seed_orders_a_segmented_text_shorter_first() { + let mut state = 999u64; + // rustar's alphabet: bases 0..=3, N = 4, spacer = 5. + let text: Vec = (0..8192) + .map(|i| { + if i % 617 < 9 { + 5 + } else { + (lcg(&mut state) % 5) as u8 + } + }) + .collect(); + let ends = spacer_ends(&text); + let lp = SegmentedText::from_ends(text.len(), ends); + let acgt: Vec = (0..text.len() as u32) + .filter(|&p| text[p as usize] < 4) + .collect(); + check_sorted(&text, &lp, acgt); + } + + #[test] + fn seed_orders_a_segmented_text_star_order() { + let mut state = 4242u64; + let text: Vec = (0..8192) + .map(|i| { + if i % 411 < 7 { + 5 + } else { + (lcg(&mut state) % 5) as u8 + } + }) + .collect(); + let ends = spacer_ends(&text); + let lp = StarSegmented { + inner: SegmentedText::from_ends(text.len(), ends), + }; + let acgt: Vec = (0..text.len() as u32) + .filter(|&p| text[p as usize] < 4) + .collect(); + check_sorted(&text, &lp, acgt); + } + + #[test] + fn seed_policy_bounds_ranked_text_allocation() { + let text = vec![0u8, 2, 0, 2, 0, 2, 0, 2]; + assert!( + seed_params(&text, PackedPrefixSeedPolicy::Disabled).is_none(), + "disabled policy must not prepare a seed", + ); + assert!( + seed_params(&text, PackedPrefixSeedPolicy::DenseAlphabetOnly).is_none(), + "dense-only policy must decline a gapped alphabet", + ); + assert!( + seed_params( + &text, + PackedPrefixSeedPolicy::remap(text.len().saturating_sub(1)), + ) + .is_none(), + "a remap over budget must decline", + ); + assert!( + seed_params(&text, PackedPrefixSeedPolicy::remap(text.len())).is_some(), + "an exactly budgeted remap should be available", + ); + + check_sorted_with_policy( + &text, + &PlainText::new(text.len()), + (0..text.len() as u32).collect(), + PackedPrefixSeedPolicy::remap(text.len()), + ); + } + + #[test] + fn seed_requires_a_reserved_boundary_code() { + let with_sentinel: Vec = (0..=254).collect(); + let packer = seed_params(&with_sentinel, PackedPrefixSeedPolicy::DenseAlphabetOnly) + .expect("255 byte values leave one u8 code for the boundary"); + assert_eq!(packer.bits, 8); + assert_eq!(packer.k, 8); + + let all_bytes: Vec = (0..=255).collect(); + assert!( + seed_params(&all_bytes, PackedPrefixSeedPolicy::DenseAlphabetOnly).is_none(), + "all 256 byte values leave no boundary code and must fall back", + ); + } + + #[test] + fn seed_lcps_are_exact_at_key_depths_and_boundary_ties() { + // Three isolated pairs share k-1, k, and k+1 symbols, and the last pair + // is identical through a simultaneous segment boundary. With alphabet + // 0..=2 plus one reserved boundary code, fields are two bits and k=32. + let k = 32usize; + let segments = [ + [vec![0; k - 1], vec![1]].concat(), + [vec![0; k - 1], vec![2]].concat(), + [vec![1; k], vec![0]].concat(), + [vec![1; k], vec![2]].concat(), + [vec![2; k + 1], vec![0]].concat(), + [vec![2; k + 1], vec![1]].concat(), + vec![0; k + 2], + vec![0; k + 2], + ]; + let lengths: Vec = segments.iter().map(Vec::len).collect(); + let mut starts = Vec::with_capacity(segments.len()); + let mut text = Vec::new(); + for segment in &segments { + starts.push(text.len() as u32); + text.extend_from_slice(segment); + } + let packer = seed_params(&text, PackedPrefixSeedPolicy::DenseAlphabetOnly).unwrap(); + assert_eq!(packer.k, k); + + let shorter_first = SegmentedText::from_lengths(text.len(), &lengths); + check_sorted(&text, &shorter_first, starts.clone()); + + let longer_first = StarSegmented { + inner: SegmentedText::from_lengths(text.len(), &lengths), + }; + check_sorted(&text, &longer_first, starts); + } + + /// Segment ends for a spacer-separated text: one past each maximal + /// non-spacer run, closing at the text length. + fn spacer_ends(text: &[u8]) -> Vec { + let mut ends = Vec::new(); + let mut in_run = false; + for (i, &b) in text.iter().enumerate() { + if b == 5 { + if in_run { + ends.push(i as u64); + in_run = false; + } + } else { + in_run = true; + } + } + if ends.last() != Some(&(text.len() as u64)) { + ends.push(text.len() as u64); + } + ends + } +}