Skip to main content

hermes_core/segment/builder/
graph_bisection.rs

1//! Recursive Graph Bisection (BP) for BMP document ordering.
2//!
3//! Based on Dhulipala et al. (KDD 2016) and Mackenzie et al. — the same
4//! algorithm used in Lucene and PISA for document reordering.
5//!
6//! Directly optimizes log-gap cost: docs sharing dimensions end up in the
7//! same BMP blocks, producing tight upper bounds and effective pruning.
8//!
9//! Memory is budgeted: CSR terms plus roughly 32 bytes/document of graph
10//! scratch, with lazily initialized direct-index term-degree arrays.
11
12#[cfg(feature = "native")]
13use rayon::prelude::*;
14use rustc_hash::FxHashMap;
15
16/// Per-partition left/right term degrees with direct compact-term indexing.
17///
18/// Recursive BP used to zero two `num_terms`-long vectors at every node. At
19/// 100k vocabulary terms and hundreds of thousands of fine partitions, that
20/// turns into a large amount of memory traffic unrelated to actual postings.
21/// A one-bit initialization map lets us retain array-speed lookups while only
22/// touching degree slots present in the current partition.
23struct TermDegrees {
24    values: Vec<std::mem::MaybeUninit<[u32; 2]>>,
25    initialized: Vec<u64>,
26}
27
28impl TermDegrees {
29    fn new(num_terms: usize) -> Self {
30        let mut values = Vec::with_capacity(num_terms);
31        values.resize_with(num_terms, std::mem::MaybeUninit::uninit);
32        Self {
33            values,
34            initialized: vec![0; num_terms.div_ceil(64)],
35        }
36    }
37
38    #[inline]
39    fn entry_mut(&mut self, term: usize) -> &mut [u32; 2] {
40        let word = term / 64;
41        let mask = 1u64 << (term % 64);
42        if self.initialized[word] & mask == 0 {
43            self.values[term].write([0, 0]);
44            self.initialized[word] |= mask;
45        }
46        // SAFETY: the bit above is set only after writing this exact slot.
47        unsafe { self.values[term].assume_init_mut() }
48    }
49
50    #[inline]
51    fn get(&self, term: usize) -> [u32; 2] {
52        let word = term / 64;
53        let mask = 1u64 << (term % 64);
54        if self.initialized[word] & mask == 0 {
55            return [0, 0];
56        }
57        // SAFETY: an initialized bit is published only after the slot write;
58        // scoring reads degrees after construction, with no concurrent writes.
59        unsafe { *self.values[term].assume_init_ref() }
60    }
61}
62
63// ── Forward index (CSR) ──────────────────────────────────────────────────
64
65/// Forward index in CSR format: doc `d`'s terms are `terms[offsets[d]..offsets[d+1]]`.
66///
67/// Term IDs are remapped to compact range `0..num_terms` for flat-array degree tracking.
68pub(crate) struct ForwardIndex {
69    terms: Vec<u32>,
70    /// u64, not u32: a 58M-doc / ~85-dims-per-doc reorder pass carries ~4.9B
71    /// postings — u32 prefix sums wrapped and the CSR carving panicked
72    /// (prod 2026-07-14, "mid > len"). The old 8 GB memory budget masked it
73    /// by dropping dims below the u32 limit.
74    offsets: Vec<u64>,
75    pub num_terms: usize,
76}
77
78/// Build CSR offsets (prefix sums) from per-entity counts. u64 output — the
79/// sum of counts legitimately exceeds u32::MAX on large reorder passes.
80fn build_csr_offsets(counts: &[u32]) -> Vec<u64> {
81    let mut offsets = Vec::with_capacity(counts.len() + 1);
82    offsets.push(0u64);
83    for &c in counts {
84        offsets.push(offsets.last().unwrap() + c as u64);
85    }
86    offsets
87}
88
89impl ForwardIndex {
90    #[inline]
91    pub fn num_docs(&self) -> usize {
92        if self.offsets.is_empty() {
93            0
94        } else {
95            self.offsets.len() - 1
96        }
97    }
98
99    #[inline]
100    fn doc_terms(&self, doc: usize) -> &[u32] {
101        let start = self.offsets[doc] as usize;
102        let end = self.offsets[doc + 1] as usize;
103        &self.terms[start..end]
104    }
105
106    /// Total postings in the forward index.
107    pub fn total_postings(&self) -> u64 {
108        self.offsets.last().copied().unwrap_or(0)
109    }
110}
111
112/// Build virtual→real and real→virtual vid maps from a BMP doc map.
113///
114/// A virtual slot is real iff its doc-map entry is not the `u32::MAX` padding
115/// sentinel. Realness must come from the doc map itself: block-copy merged
116/// segments carry each source's tail padding as *interior* padding, so
117/// `vid < num_real_docs` does NOT identify real docs there.
118///
119/// Returns `(virtual_to_real, real_to_virtual)` where `virtual_to_real[vid]`
120/// is the dense real index or `u32::MAX` for padding.
121pub(crate) fn build_vid_maps(bmp: &crate::segment::reader::bmp::BmpIndex) -> (Vec<u32>, Vec<u32>) {
122    let ids = bmp.doc_map_ids_slice();
123    let num_virtual = bmp.num_virtual_docs as usize;
124    let mut virtual_to_real = vec![u32::MAX; num_virtual];
125    let mut real_to_virtual = Vec::with_capacity(bmp.num_real_docs() as usize);
126    for (vid, (slot, chunk)) in virtual_to_real
127        .iter_mut()
128        .zip(ids.as_chunks::<4>().0)
129        .enumerate()
130    {
131        let doc_id = u32::from_le_bytes(*chunk);
132        if doc_id != u32::MAX {
133            *slot = real_to_virtual.len() as u32;
134            real_to_virtual.push(vid as u32);
135        }
136    }
137    if real_to_virtual.len() != bmp.num_real_docs() as usize {
138        log::warn!(
139            "[reorder] BMP doc map has {} real slots but footer says num_real_docs={} — trusting the doc map",
140            real_to_virtual.len(),
141            bmp.num_real_docs(),
142        );
143    }
144    (virtual_to_real, real_to_virtual)
145}
146
147/// One (source, block) unit of forward-index construction. Because
148/// [`build_vid_maps`] assigns real ids in ascending vid order, a block's real
149/// docs form the contiguous per-source range
150/// `real_start..real_start + real_len` — blocks can be processed in parallel
151/// with disjoint output slices.
152struct BlockJob {
153    src: u32,
154    block_id: u32,
155    /// Per-source real index of the block's first real doc.
156    real_start: u32,
157    /// Number of real (non-padding) docs in the block.
158    real_len: u32,
159}
160
161/// Enumerate jobs in (source, block) order — cumulative `real_len` tiles the
162/// global real-id space `0..total_docs` exactly.
163fn build_block_jobs(
164    bmps: &[&crate::segment::reader::bmp::BmpIndex],
165    vid_maps: &[(Vec<u32>, Vec<u32>)],
166) -> Vec<BlockJob> {
167    let total_blocks: usize = bmps.iter().map(|b| b.num_blocks as usize).sum();
168    let mut jobs = Vec::with_capacity(total_blocks);
169    for (src, (bmp, (v2r, _))) in bmps.iter().zip(vid_maps).enumerate() {
170        let block_size = bmp.bmp_block_size as usize;
171        let mut real_cursor = 0u32;
172        for block_id in 0..bmp.num_blocks as usize {
173            let vid_start = block_id * block_size;
174            let vid_end = ((block_id + 1) * block_size).min(v2r.len());
175            let real_len = v2r[vid_start..vid_end]
176                .iter()
177                .filter(|&&r| r != u32::MAX)
178                .count() as u32;
179            jobs.push(BlockJob {
180                src: src as u32,
181                block_id: block_id as u32,
182                real_start: real_cursor,
183                real_len,
184            });
185            real_cursor += real_len;
186        }
187    }
188    jobs
189}
190
191/// Merge the smaller df map into the larger (rayon reduce combiner).
192#[cfg(feature = "native")]
193fn merge_df_maps(
194    mut a: FxHashMap<u32, usize>,
195    mut b: FxHashMap<u32, usize>,
196) -> FxHashMap<u32, usize> {
197    if a.len() < b.len() {
198        std::mem::swap(&mut a, &mut b);
199    }
200    for (k, v) in b {
201        *a.entry(k).or_insert(0) += v;
202    }
203    a
204}
205
206/// Build forward index from BmpIndex sources (single or multi-source).
207///
208/// Documents are identified by dense *real* indices assigned sequentially
209/// across sources: source 0 gets 0..n0, source 1 gets n0..n0+n1, etc., where
210/// each n is the source's real (non-padding) doc count derived from its doc
211/// map via [`build_vid_maps`]. Returns `(forward_index, per_source_real_doc_counts)`.
212///
213/// Filters dims with doc_freq outside `[min_doc_freq, max_doc_freq]`.
214/// If the estimated forward index memory exceeds `memory_budget_bytes`, the
215/// highest-frequency dims are dropped to stay within budget. This prevents OOM
216/// for huge segments at the cost of slightly reduced reorder quality.
217///
218/// Remaps term IDs to compact range for flat-array degree tracking.
219#[cfg(test)]
220pub(crate) fn build_forward_index_from_bmps(
221    bmps: &[&crate::segment::reader::bmp::BmpIndex],
222    min_doc_freq: usize,
223    max_doc_freq: usize,
224    memory_budget_bytes: usize,
225) -> (ForwardIndex, Vec<usize>) {
226    let vid_maps: Vec<(Vec<u32>, Vec<u32>)> = bmps.iter().map(|b| build_vid_maps(b)).collect();
227    build_forward_index_from_bmps_with_maps(
228        bmps,
229        &vid_maps,
230        min_doc_freq,
231        max_doc_freq,
232        memory_budget_bytes,
233    )
234}
235
236/// Variant for reorder callers that already need the virtual/real maps during
237/// output encoding. Reusing them avoids a second full document-map scan and a
238/// duplicate real-to-virtual allocation on very large segments.
239pub(crate) fn build_forward_index_from_bmps_with_maps(
240    bmps: &[&crate::segment::reader::bmp::BmpIndex],
241    vid_maps: &[(Vec<u32>, Vec<u32>)],
242    min_doc_freq: usize,
243    max_doc_freq: usize,
244    memory_budget_bytes: usize,
245) -> (ForwardIndex, Vec<usize>) {
246    debug_assert_eq!(bmps.len(), vid_maps.len());
247    let source_doc_counts: Vec<usize> = vid_maps.iter().map(|(_, r2v)| r2v.len()).collect();
248    let total_docs: usize = source_doc_counts.iter().sum();
249
250    if total_docs == 0 {
251        return (
252            ForwardIndex {
253                terms: Vec::new(),
254                offsets: Vec::new(),
255                num_terms: 0,
256            },
257            source_doc_counts,
258        );
259    }
260
261    // Job list: one entry per (source, block). Real ids are assigned in
262    // ascending vid order (see build_vid_maps), so each block owns a
263    // contiguous real-id range — every phase below can process blocks in
264    // parallel, writing disjoint slices.
265    let jobs = build_block_jobs(bmps, vid_maps);
266
267    // Phase 1: count doc freq per dimension across all sources + assign compact IDs
268    let count_block_df = |job: &BlockJob, acc: &mut FxHashMap<u32, usize>| {
269        let bmp = bmps[job.src as usize];
270        let (v2r, _) = &vid_maps[job.src as usize];
271        let block_size = bmp.bmp_block_size as usize;
272        for (dim_id, postings) in bmp.iter_block_terms(job.block_id) {
273            let mut n = 0usize;
274            for p in postings {
275                let vid = job.block_id as usize * block_size + p.local_slot as usize;
276                if v2r[vid] != u32::MAX && p.impact > 0 {
277                    n += 1;
278                }
279            }
280            if n > 0 {
281                *acc.entry(dim_id).or_insert(0) += n;
282            }
283        }
284    };
285    #[cfg(feature = "native")]
286    let dim_df: FxHashMap<u32, usize> = jobs
287        .par_iter()
288        .fold(FxHashMap::default, |mut acc, job| {
289            count_block_df(job, &mut acc);
290            acc
291        })
292        .reduce(FxHashMap::default, merge_df_maps);
293    #[cfg(not(feature = "native"))]
294    let dim_df: FxHashMap<u32, usize> = {
295        let mut acc = FxHashMap::default();
296        for job in &jobs {
297            count_block_df(job, &mut acc);
298        }
299        acc
300    };
301
302    // Filter dims by [min_doc_freq, max_doc_freq] range
303    let mut eligible: Vec<(u32, usize)> = dim_df
304        .iter()
305        .filter(|&(_, df)| *df >= min_doc_freq && *df <= max_doc_freq)
306        .map(|(&dim_id, &df)| (dim_id, df))
307        .collect();
308    drop(dim_df);
309
310    // Memory budget: estimate forward index + bisection scratch.
311    // Peak ≈ 4*total_postings (terms array) + 32*total_docs (u64 offsets,
312    //   counts, docs, gains, indices, new_left, new_right) + 8*num_terms
313    //   (degree arrays).
314    let total_postings_est = eligible
315        .iter()
316        .fold(0usize, |total, (_, df)| total.saturating_add(*df));
317    let entity_scratch_bytes = total_docs.saturating_mul(32);
318    let estimated_bytes = total_postings_est
319        .saturating_mul(4)
320        .saturating_add(entity_scratch_bytes)
321        .saturating_add(eligible.len().saturating_mul(8));
322
323    if estimated_bytes > memory_budget_bytes && !eligible.is_empty() {
324        // Sort by df ascending — keep discriminative low-df dims first,
325        // drop highest-df dims which contribute the most postings.
326        eligible.sort_by_key(|&(_, df)| df);
327
328        // Account for each retained term together with its postings. The old
329        // calculation charged the eight-byte degree slot for every candidate
330        // before deciding how many to retain; a large rare-term vocabulary
331        // could therefore make the target zero even when a useful subset fit.
332        let mut used_bytes = entity_scratch_bytes;
333        let mut cum = 0usize;
334        let mut keep_count = 0;
335        for &(_, df) in &eligible {
336            let term_bytes = df.saturating_mul(4).saturating_add(8);
337            if term_bytes > memory_budget_bytes.saturating_sub(used_bytes) {
338                break;
339            }
340            used_bytes = used_bytes.saturating_add(term_bytes);
341            cum = cum.saturating_add(df);
342            keep_count += 1;
343        }
344
345        let dropped = eligible.len() - keep_count;
346        eligible.truncate(keep_count);
347
348        log::warn!(
349            "[reorder] memory budget {:.0} MB: estimated {:.0} MB, dropped {} highest-df dims, keeping {} ({} postings)",
350            memory_budget_bytes as f64 / (1024.0 * 1024.0),
351            estimated_bytes as f64 / (1024.0 * 1024.0),
352            dropped,
353            keep_count,
354            cum,
355        );
356    }
357
358    if eligible.is_empty() {
359        // The caller emits an identity permutation when there is no graph
360        // signal. Avoid allocating per-document counts and u64 offsets only
361        // to discover that the terms array is empty, especially when the
362        // configured budget is below the fixed document scratch cost.
363        return (
364            ForwardIndex {
365                terms: Vec::new(),
366                offsets: Vec::new(),
367                num_terms: 0,
368            },
369            source_doc_counts,
370        );
371    }
372
373    let mut term_remap: FxHashMap<u32, u32> = FxHashMap::default();
374    for &(dim_id, _) in &eligible {
375        let compact_id = term_remap.len() as u32;
376        term_remap.insert(dim_id, compact_id);
377    }
378    let num_active_terms = term_remap.len();
379    drop(eligible);
380
381    // Phase 2: count terms per doc (filtered) — per-block disjoint slices
382    let mut counts = vec![0u32; total_docs];
383    let fill_block_counts = |job: &BlockJob, out: &mut [u32]| {
384        let bmp = bmps[job.src as usize];
385        let (v2r, _) = &vid_maps[job.src as usize];
386        let block_size = bmp.bmp_block_size as usize;
387        for (dim_id, postings) in bmp.iter_block_terms(job.block_id) {
388            if !term_remap.contains_key(&dim_id) {
389                continue;
390            }
391            for p in postings {
392                let vid = job.block_id as usize * block_size + p.local_slot as usize;
393                let real = v2r[vid];
394                if real != u32::MAX && p.impact > 0 {
395                    out[(real - job.real_start) as usize] += 1;
396                }
397            }
398        }
399    };
400    {
401        let mut slices: Vec<(&BlockJob, &mut [u32])> = Vec::with_capacity(jobs.len());
402        let mut rest: &mut [u32] = &mut counts;
403        for job in &jobs {
404            let (head, tail) = rest.split_at_mut(job.real_len as usize);
405            slices.push((job, head));
406            rest = tail;
407        }
408        #[cfg(feature = "native")]
409        slices
410            .into_par_iter()
411            .for_each(|(job, out)| fill_block_counts(job, out));
412        #[cfg(not(feature = "native"))]
413        for (job, out) in slices {
414            fill_block_counts(job, out);
415        }
416    }
417
418    // Phase 3: build CSR offsets (u64 — sums exceed u32::MAX at scale)
419    let offsets = build_csr_offsets(&counts);
420    let total = *offsets.last().unwrap() as usize;
421    drop(counts);
422
423    // Phase 4: fill terms (compact IDs) — each block writes the contiguous
424    // terms range covering its real docs; per-doc write cursors are local.
425    let mut terms = vec![0u32; total];
426    let fill_block_terms = |job: &BlockJob, global_real_start: usize, out: &mut [u32]| {
427        let bmp = bmps[job.src as usize];
428        let (v2r, _) = &vid_maps[job.src as usize];
429        let block_size = bmp.bmp_block_size as usize;
430        // local_slot is u8, so a block never holds more than 256 real docs
431        assert!(job.real_len as usize <= 256, "BMP block exceeds 256 docs");
432        let mut cursor = [0u32; 256];
433        let base = offsets[global_real_start] as usize;
434        for (dim_id, postings) in bmp.iter_block_terms(job.block_id) {
435            let Some(&compact) = term_remap.get(&dim_id) else {
436                continue;
437            };
438            for p in postings {
439                let vid = job.block_id as usize * block_size + p.local_slot as usize;
440                let real = v2r[vid];
441                if real != u32::MAX && p.impact > 0 {
442                    let local = (real - job.real_start) as usize;
443                    let pos =
444                        offsets[global_real_start + local] as usize - base + cursor[local] as usize;
445                    out[pos] = compact;
446                    cursor[local] += 1;
447                }
448            }
449        }
450    };
451    {
452        let mut slices: Vec<(&BlockJob, usize, &mut [u32])> = Vec::with_capacity(jobs.len());
453        let mut rest: &mut [u32] = &mut terms;
454        let mut global_real = 0usize;
455        for job in &jobs {
456            let len =
457                (offsets[global_real + job.real_len as usize] - offsets[global_real]) as usize;
458            let (head, tail) = rest.split_at_mut(len);
459            slices.push((job, global_real, head));
460            rest = tail;
461            global_real += job.real_len as usize;
462        }
463        #[cfg(feature = "native")]
464        slices
465            .into_par_iter()
466            .for_each(|(job, g, out)| fill_block_terms(job, g, out));
467        #[cfg(not(feature = "native"))]
468        for (job, g, out) in slices {
469            fill_block_terms(job, g, out);
470        }
471    }
472
473    (
474        ForwardIndex {
475            terms,
476            offsets,
477            num_terms: num_active_terms,
478        },
479        source_doc_counts,
480    )
481}
482
483/// Build a forward index over BLOCKS (one entity per block, its terms = the
484/// block's header dim list). Used by block-level reorder: BP over blocks is
485/// ~block_size× cheaper than over records and only needs to decide superblock
486/// assignment. Blocks are numbered globally across sources in source order.
487///
488/// Dims appearing in fewer than 2 blocks carry no clustering signal and are
489/// dropped; the memory budget applies as in the record-level builder.
490pub(crate) fn build_forward_index_from_blocks(
491    bmps: &[&crate::segment::reader::bmp::BmpIndex],
492    memory_budget_bytes: usize,
493) -> ForwardIndex {
494    let total_blocks: usize = bmps.iter().map(|b| b.num_blocks as usize).sum();
495    if total_blocks == 0 {
496        return ForwardIndex {
497            terms: Vec::new(),
498            offsets: Vec::new(),
499            num_terms: 0,
500        };
501    }
502
503    // (source, block) pairs in global block order — the parallel unit.
504    let blocks: Vec<(u32, u32)> = bmps
505        .iter()
506        .enumerate()
507        .flat_map(|(src, bmp)| (0..bmp.num_blocks).map(move |b| (src as u32, b)))
508        .collect();
509
510    // Phase 1: dim → number of blocks containing it (block-level df)
511    let count_block_bf = |&(src, block_id): &(u32, u32), acc: &mut FxHashMap<u32, usize>| {
512        for (dim_id, _) in bmps[src as usize].iter_block_terms(block_id) {
513            *acc.entry(dim_id).or_insert(0) += 1;
514        }
515    };
516    #[cfg(feature = "native")]
517    let dim_bf: FxHashMap<u32, usize> = blocks
518        .par_iter()
519        .fold(FxHashMap::default, |mut acc, b| {
520            count_block_bf(b, &mut acc);
521            acc
522        })
523        .reduce(FxHashMap::default, merge_df_maps);
524    #[cfg(not(feature = "native"))]
525    let dim_bf: FxHashMap<u32, usize> = {
526        let mut acc = FxHashMap::default();
527        for b in &blocks {
528            count_block_bf(b, &mut acc);
529        }
530        acc
531    };
532
533    let max_bf = (total_blocks as f64 * 0.9) as usize;
534    let mut eligible: Vec<(u32, usize)> = dim_bf
535        .iter()
536        .filter(|&(_, bf)| *bf >= 2 && *bf <= max_bf.max(2))
537        .map(|(&dim_id, &bf)| (dim_id, bf))
538        .collect();
539    drop(dim_bf);
540
541    let total_postings_est = eligible
542        .iter()
543        .fold(0usize, |total, (_, bf)| total.saturating_add(*bf));
544    let entity_scratch_bytes = total_blocks.saturating_mul(32);
545    let estimated_bytes = total_postings_est
546        .saturating_mul(4)
547        .saturating_add(entity_scratch_bytes)
548        .saturating_add(eligible.len().saturating_mul(8));
549    if estimated_bytes > memory_budget_bytes && !eligible.is_empty() {
550        eligible.sort_by_key(|&(_, bf)| bf);
551        let mut used_bytes = entity_scratch_bytes;
552        let mut cum = 0usize;
553        let mut keep = 0;
554        for &(_, bf) in &eligible {
555            let term_bytes = bf.saturating_mul(4).saturating_add(8);
556            if term_bytes > memory_budget_bytes.saturating_sub(used_bytes) {
557                break;
558            }
559            used_bytes = used_bytes.saturating_add(term_bytes);
560            cum = cum.saturating_add(bf);
561            keep += 1;
562        }
563        log::warn!(
564            "[reorder] block-level fwd index over budget — dropped {} highest-bf dims",
565            eligible.len() - keep,
566        );
567        eligible.truncate(keep);
568    }
569
570    if eligible.is_empty() {
571        return ForwardIndex {
572            terms: Vec::new(),
573            offsets: Vec::new(),
574            num_terms: 0,
575        };
576    }
577
578    let mut term_remap: FxHashMap<u32, u32> = FxHashMap::default();
579    for &(dim_id, _) in &eligible {
580        let compact = term_remap.len() as u32;
581        term_remap.insert(dim_id, compact);
582    }
583    let num_terms = term_remap.len();
584    drop(eligible);
585
586    // Phase 2+3: counts and CSR fill — one entity per block, so each block
587    // maps to a single count cell and a contiguous terms range.
588    let count_remapped = |&(src, block_id): &(u32, u32)| -> u32 {
589        bmps[src as usize]
590            .iter_block_terms(block_id)
591            .filter(|(dim_id, _)| term_remap.contains_key(dim_id))
592            .count() as u32
593    };
594    #[cfg(feature = "native")]
595    let counts: Vec<u32> = blocks.par_iter().map(count_remapped).collect();
596    #[cfg(not(feature = "native"))]
597    let counts: Vec<u32> = blocks.iter().map(count_remapped).collect();
598
599    let offsets = build_csr_offsets(&counts);
600    let total = *offsets.last().unwrap() as usize;
601    drop(counts);
602
603    let mut terms = vec![0u32; total];
604    let fill_block = |&(src, block_id): &(u32, u32), out: &mut [u32]| {
605        let mut n = 0usize;
606        for (dim_id, _) in bmps[src as usize].iter_block_terms(block_id) {
607            if let Some(&compact) = term_remap.get(&dim_id) {
608                out[n] = compact;
609                n += 1;
610            }
611        }
612    };
613    {
614        let mut slices: Vec<(&(u32, u32), &mut [u32])> = Vec::with_capacity(blocks.len());
615        let mut rest: &mut [u32] = &mut terms;
616        for (gb, b) in blocks.iter().enumerate() {
617            let len = (offsets[gb + 1] - offsets[gb]) as usize;
618            let (head, tail) = rest.split_at_mut(len);
619            slices.push((b, head));
620            rest = tail;
621        }
622        #[cfg(feature = "native")]
623        slices
624            .into_par_iter()
625            .for_each(|(b, out)| fill_block(b, out));
626        #[cfg(not(feature = "native"))]
627        for (b, out) in slices {
628            fill_block(b, out);
629        }
630    }
631
632    ForwardIndex {
633        terms,
634        offsets,
635        num_terms,
636    }
637}
638
639// ── Recursive Graph Bisection ────────────────────────────────────────────
640
641/// CPU/depth budget for a BP pass. BP is an anytime algorithm: stopping at
642/// any depth or deadline still yields a valid permutation, and because the
643/// output layout becomes the next pass's input order, repeated budgeted
644/// passes warm-start and deepen (top levels converge in ~0 swaps, the budget
645/// flows to deeper levels).
646#[derive(Clone, Copy, Debug, Default)]
647pub struct BpBudget {
648    /// Stop recursion at partitions of at most this many docs instead of
649    /// descending to block granularity. `None` = full depth. Capping at
650    /// superblock granularity (superblock_size × block_size docs) keeps most
651    /// of the superblock-pruning win at ~⅓ less depth.
652    pub min_partition_docs: Option<usize>,
653    /// Wall-clock cap for the whole BP computation. The pass ends cleanly at
654    /// the deadline with whatever depth it reached (`converged = false`).
655    /// Ignored on wasm (no monotonic clock).
656    pub time_budget: Option<std::time::Duration>,
657}
658
659impl BpBudget {
660    /// Unbudgeted: full depth, no deadline.
661    pub fn full() -> Self {
662        Self::default()
663    }
664}
665
666/// Recursive graph bisection. Returns `(perm, converged)` where
667/// `perm[new_pos] = old_index` and `converged` is false iff the wall-clock
668/// budget ended the pass before it finished (a depth cap alone is a chosen
669/// target, not an interruption — it reports converged).
670///
671/// `min_partition_size` should be the BMP block_size (64).
672/// `max_iters` controls convergence (20 is standard).
673///
674/// Term IDs in the forward index must be compact (0..num_terms) so we can
675/// use flat arrays for O(1) degree lookups instead of hash maps.
676pub(crate) fn graph_bisection(
677    fwd: &ForwardIndex,
678    min_partition_size: usize,
679    max_iters: usize,
680    budget: BpBudget,
681) -> (Vec<u32>, bool) {
682    let n = fwd.num_docs();
683    if n == 0 {
684        return (Vec::new(), true);
685    }
686
687    let effective_min_partition = budget
688        .min_partition_docs
689        .unwrap_or(0)
690        .max(min_partition_size);
691
692    let mut docs: Vec<u32> = (0..n as u32).collect();
693    let depth = if effective_min_partition > 0 {
694        ((n as f64) / (effective_min_partition as f64))
695            .log2()
696            .ceil() as usize
697    } else {
698        0
699    };
700    let log_table = build_log_table(4096);
701
702    log::debug!(
703        "BP graph_bisection: n={}, min_partition={}, max_iters={}, depth=~{}, time_budget={:?}",
704        n,
705        effective_min_partition,
706        max_iters,
707        depth,
708        budget.time_budget,
709    );
710
711    #[cfg(feature = "native")]
712    let deadline = budget.time_budget.map(|duration| {
713        let now = std::time::Instant::now();
714        now.checked_add(duration).unwrap_or(now)
715    });
716    #[cfg(not(feature = "native"))]
717    let deadline: Option<()> = None;
718    #[cfg(not(feature = "native"))]
719    let _ = deadline;
720
721    let exhausted = std::sync::atomic::AtomicBool::new(false);
722    #[cfg(feature = "native")]
723    bisect(
724        &mut docs,
725        fwd,
726        effective_min_partition,
727        max_iters,
728        &log_table,
729        deadline,
730        &exhausted,
731    );
732    #[cfg(not(feature = "native"))]
733    bisect(
734        &mut docs,
735        fwd,
736        effective_min_partition,
737        max_iters,
738        &log_table,
739        None,
740        &exhausted,
741    );
742
743    let converged = !exhausted.load(std::sync::atomic::Ordering::Relaxed);
744    if !converged {
745        log::info!(
746            "BP graph_bisection: wall-clock budget {:?} exhausted at n={} — emitting partial (still valid) permutation",
747            budget.time_budget,
748            n,
749        );
750    }
751    (docs, converged)
752}
753
754/// Recursive bisection of a document slice.
755///
756/// Uses flat `Vec<u32>` degree arrays indexed by compact term_id for cache-friendly
757/// O(1) lookups (vs FxHashMap which has poor cache locality at scale).
758///
759/// Gain computation is parallelized via rayon for large partitions (n > 4096).
760/// Adaptive iteration count reduces work at top levels where coarse splits
761/// converge faster and dominate total runtime.
762fn bisect(
763    docs: &mut [u32],
764    fwd: &ForwardIndex,
765    min_partition_size: usize,
766    max_iters: usize,
767    log_table: &[f32],
768    #[cfg(feature = "native")] deadline: Option<std::time::Instant>,
769    #[cfg(not(feature = "native"))] deadline: Option<()>,
770    exhausted: &std::sync::atomic::AtomicBool,
771) {
772    let n = docs.len();
773    if n <= min_partition_size {
774        return;
775    }
776    // Anytime cutoff: leave this subtree in its current (valid) order.
777    if exhausted.load(std::sync::atomic::Ordering::Relaxed) {
778        return;
779    }
780    #[cfg(feature = "native")]
781    if let Some(dl) = deadline
782        && std::time::Instant::now() >= dl
783    {
784        exhausted.store(true, std::sync::atomic::Ordering::Relaxed);
785        return;
786    }
787    #[cfg(not(feature = "native"))]
788    let _ = deadline;
789
790    let mid = n / 2;
791    let nt = fwd.num_terms;
792
793    // Adaptive iteration count: large partitions converge faster with
794    // coarse splits, so fewer refinement passes suffice. The fine-grained
795    // clustering is handled by deeper recursion levels with full iterations.
796    let effective_iters = if n > 100_000 {
797        max_iters.min(12)
798    } else {
799        max_iters
800    };
801
802    // Compact term IDs permit direct indexing. Slots are initialized lazily so
803    // deep partitions do not zero the full vocabulary on every recursive node.
804    let mut degrees = TermDegrees::new(nt);
805
806    for (i, &doc) in docs.iter().enumerate() {
807        let side = usize::from(i >= mid);
808        for &term in fwd.doc_terms(doc as usize) {
809            degrees.entry_mut(term as usize)[side] += 1;
810        }
811    }
812
813    // Scratch buffers reused across iterations
814    let mut gains: Vec<f32> = vec![0.0; n];
815    let mut indices: Vec<usize> = (0..n).collect();
816    let mut new_left: Vec<u32> = Vec::with_capacity(mid);
817    let mut new_right: Vec<u32> = Vec::with_capacity(n - mid);
818
819    for iter in 0..effective_iters {
820        // Anytime cutoff between refinement passes: keep the current split.
821        #[cfg(feature = "native")]
822        if let Some(dl) = deadline
823            && std::time::Instant::now() >= dl
824        {
825            exhausted.store(true, std::sync::atomic::Ordering::Relaxed);
826            break;
827        }
828        // Compute gain for each document (approx_1 from Dhulipala et al.)
829        // Parallelized for large partitions where per-doc work dominates.
830        compute_gains(docs, fwd, mid, &degrees, log_table, &mut gains);
831
832        // Partition: the `mid` LOWEST keys (strongest left affinity) go left
833        indices.clear();
834        indices.extend(0..n);
835        indices.select_nth_unstable_by(mid, |&a, &b| {
836            gains[a].total_cmp(&gains[b]).then_with(|| a.cmp(&b))
837        });
838
839        // Apply partition, update degree arrays for swapped docs
840        new_left.clear();
841        new_right.clear();
842        let mut swap_count: usize = 0;
843
844        for (rank, &idx) in indices.iter().enumerate() {
845            let doc = docs[idx];
846            let was_left = idx < mid;
847            let now_left = rank < mid;
848
849            if now_left {
850                new_left.push(doc);
851            } else {
852                new_right.push(doc);
853            }
854
855            if was_left != now_left {
856                swap_count += 1;
857                for &term in fwd.doc_terms(doc as usize) {
858                    let degree = degrees.entry_mut(term as usize);
859                    if was_left {
860                        degree[0] -= 1;
861                        degree[1] += 1;
862                    } else {
863                        degree[1] -= 1;
864                        degree[0] += 1;
865                    }
866                }
867            }
868        }
869
870        docs[..mid].copy_from_slice(&new_left);
871        docs[mid..].copy_from_slice(&new_right);
872
873        if swap_count == 0 {
874            break;
875        }
876
877        // Early termination: if < 0.5% of docs swapped, partition is stable
878        if iter > 2 && swap_count < n / 200 {
879            break;
880        }
881
882        // Cooling: break early if gains are negligible
883        if iter > 5 {
884            let max_abs_gain = gains
885                .iter()
886                .copied()
887                .fold(0.0f32, |max_gain, gain| max_gain.max(gain.abs()));
888            if max_abs_gain < 0.001 {
889                break;
890            }
891        }
892    }
893
894    // Drop scratch before recursion to free memory for sub-problems
895    drop(degrees);
896    drop(gains);
897    drop(indices);
898    drop(new_left);
899    drop(new_right);
900
901    let (left, right) = docs.split_at_mut(mid);
902    #[cfg(feature = "native")]
903    rayon::join(
904        || {
905            bisect(
906                left,
907                fwd,
908                min_partition_size,
909                max_iters,
910                log_table,
911                deadline,
912                exhausted,
913            )
914        },
915        || {
916            bisect(
917                right,
918                fwd,
919                min_partition_size,
920                max_iters,
921                log_table,
922                deadline,
923                exhausted,
924            )
925        },
926    );
927    #[cfg(not(feature = "native"))]
928    {
929        bisect(
930            left,
931            fwd,
932            min_partition_size,
933            max_iters,
934            log_table,
935            deadline,
936            exhausted,
937        );
938        bisect(
939            right,
940            fwd,
941            min_partition_size,
942            max_iters,
943            log_table,
944            deadline,
945            exhausted,
946        );
947    }
948}
949
950/// Compute gains for all documents, parallelized via rayon for large partitions.
951///
952/// Each doc's gain is independent: iterate its terms, accumulate the log-gap
953/// cost delta of moving it to the other side. Read-only access to degree arrays
954/// makes this embarrassingly parallel.
955#[inline(never)]
956fn compute_gains(
957    docs: &[u32],
958    fwd: &ForwardIndex,
959    mid: usize,
960    degrees: &TermDegrees,
961    log_table: &[f32],
962    gains: &mut [f32],
963) {
964    // Single coherent key: HIGH = belongs in the RIGHT half.
965    // Left docs get +approx_one(from=left, to=right) — a misplaced left doc
966    // (terms concentrated right) scores high. Right docs get
967    // -approx_one(from=right, to=left) — a misplaced right doc scores low.
968    // This matches the reference two-sided formulation (compute_gains_left /
969    // compute_gains_right with negation); ranking both halves by raw
970    // "move gain" instead made both sides' misplaced docs rank identically,
971    // so the partition step could never exchange them.
972    let gain_for_doc = |i: usize| -> f32 {
973        let doc = docs[i] as usize;
974        let in_left = i < mid;
975        let mut g = 0.0f32;
976        for &term in fwd.doc_terms(doc) {
977            let [left, right] = degrees.get(term as usize);
978            let (from, to) = if in_left {
979                (left, right)
980            } else {
981                (right, left)
982            };
983            let move_gain = fast_log2_lookup(to as usize + 2, log_table)
984                - fast_log2_lookup(from as usize, log_table)
985                - std::f32::consts::LOG2_E / (1.0 + to as f32);
986            g += if in_left { move_gain } else { -move_gain };
987        }
988        g
989    };
990
991    #[cfg(feature = "native")]
992    {
993        if docs.len() > 4096 {
994            gains
995                .par_iter_mut()
996                .enumerate()
997                .for_each(|(i, gain)| *gain = gain_for_doc(i));
998        } else {
999            for (i, gain) in gains.iter_mut().enumerate().take(docs.len()) {
1000                *gain = gain_for_doc(i);
1001            }
1002        }
1003    }
1004    #[cfg(not(feature = "native"))]
1005    {
1006        for (i, gain) in gains.iter_mut().enumerate().take(docs.len()) {
1007            *gain = gain_for_doc(i);
1008        }
1009    }
1010}
1011
1012// ── Helpers ──────────────────────────────────────────────────────────────
1013
1014/// Build precomputed log2 table for values 0..size.
1015fn build_log_table(size: usize) -> Vec<f32> {
1016    let mut table = vec![0.0f32; size];
1017    // log2(0) is undefined; use a large negative value
1018    table[0] = -10.0;
1019    for (i, entry) in table.iter_mut().enumerate().skip(1) {
1020        *entry = (i as f32).log2();
1021    }
1022    table
1023}
1024
1025/// Fast log2 with precomputed table lookup.
1026#[inline]
1027fn fast_log2_lookup(val: usize, table: &[f32]) -> f32 {
1028    if val < table.len() {
1029        table[val]
1030    } else {
1031        (val as f32).log2()
1032    }
1033}
1034
1035#[cfg(test)]
1036mod tests {
1037    use super::*;
1038
1039    #[test]
1040    fn lazy_term_degrees_initialize_only_on_first_write() {
1041        let mut degrees = TermDegrees::new(130);
1042        assert_eq!(degrees.get(65), [0, 0]);
1043        degrees.entry_mut(65)[0] += 3;
1044        degrees.entry_mut(65)[1] += 2;
1045        assert_eq!(degrees.get(65), [3, 2]);
1046        assert_eq!(degrees.get(64), [0, 0]);
1047        assert_eq!(
1048            degrees
1049                .initialized
1050                .iter()
1051                .map(|w| w.count_ones())
1052                .sum::<u32>(),
1053            1
1054        );
1055    }
1056
1057    /// Regression: CSR offsets were u32 and wrapped past 4.29B postings —
1058    /// a 58M-doc / ~85-dims-per-doc prod reorder pass (~4.9B postings)
1059    /// panicked with "mid > len" in the terms carving. The old 8 GB memory
1060    /// budget masked the overflow by dropping dims; raising the budget
1061    /// exposed it. Offsets must be u64.
1062    #[test]
1063    fn test_csr_offsets_do_not_wrap_past_u32() {
1064        let counts = [1_500_000_000u32; 3]; // 4.5B total > u32::MAX
1065        let offsets = build_csr_offsets(&counts);
1066        assert_eq!(
1067            offsets,
1068            vec![0, 1_500_000_000, 3_000_000_000, 4_500_000_000]
1069        );
1070        assert!(*offsets.last().unwrap() > u32::MAX as u64);
1071    }
1072
1073    /// Build a simple forward index from (doc_id, terms) pairs.
1074    fn make_fwd(docs: &[&[u32]], num_terms: usize) -> ForwardIndex {
1075        let mut terms = Vec::new();
1076        let mut offsets = vec![0u64];
1077        for doc_terms in docs {
1078            terms.extend_from_slice(doc_terms);
1079            offsets.push(terms.len() as u64);
1080        }
1081        ForwardIndex {
1082            terms,
1083            offsets,
1084            num_terms,
1085        }
1086    }
1087
1088    #[test]
1089    fn test_bp_empty() {
1090        let fwd = ForwardIndex {
1091            terms: Vec::new(),
1092            offsets: Vec::new(),
1093            num_terms: 0,
1094        };
1095        let (perm, _) = graph_bisection(&fwd, 4, 20, BpBudget::full());
1096        assert!(perm.is_empty());
1097    }
1098
1099    #[test]
1100    fn test_bp_small() {
1101        // 4 docs, min_partition_size=4 → no bisection, identity
1102        let fwd = make_fwd(&[&[0, 1], &[0, 2], &[1, 3], &[2, 3]], 4);
1103        let (perm, _) = graph_bisection(&fwd, 4, 20, BpBudget::full());
1104        assert_eq!(perm.len(), 4);
1105        // All docs present
1106        let mut sorted = perm.clone();
1107        sorted.sort();
1108        assert_eq!(sorted, vec![0, 1, 2, 3]);
1109    }
1110
1111    #[test]
1112    fn test_bp_clusters() {
1113        // 8 docs in 2 clear clusters:
1114        // Cluster A (docs 0-3): share terms 0, 1
1115        // Cluster B (docs 4-7): share terms 2, 3
1116        let fwd = make_fwd(
1117            &[
1118                &[0, 1],
1119                &[0, 1],
1120                &[0, 1],
1121                &[0, 1],
1122                &[2, 3],
1123                &[2, 3],
1124                &[2, 3],
1125                &[2, 3],
1126            ],
1127            4,
1128        );
1129        let (perm, _) = graph_bisection(&fwd, 4, 20, BpBudget::full());
1130        assert_eq!(perm.len(), 8);
1131
1132        // After bisection, docs from same cluster should be in same half
1133        let left: Vec<u32> = perm[..4].to_vec();
1134
1135        // Either all of cluster A is in left and B in right, or vice versa
1136        let a_in_left = left.iter().filter(|&&d| d < 4).count();
1137        let b_in_left = left.iter().filter(|&&d| d >= 4).count();
1138        assert!(
1139            (a_in_left == 4 && b_in_left == 0) || (a_in_left == 0 && b_in_left == 4),
1140            "Clusters should be separated: a_left={}, b_left={}",
1141            a_in_left,
1142            b_in_left,
1143        );
1144    }
1145
1146    #[test]
1147    fn test_bp_permutation_valid() {
1148        // 16 docs with mixed terms: terms range from 0..4 and 10..18
1149        let docs: Vec<Vec<u32>> = (0..16).map(|i| vec![i / 4, 10 + i / 2]).collect();
1150        let doc_refs: Vec<&[u32]> = docs.iter().map(|v| v.as_slice()).collect();
1151        let fwd = make_fwd(&doc_refs, 18); // max term = 10 + 15/2 = 17, so need 18
1152        let (perm, _) = graph_bisection(&fwd, 4, 20, BpBudget::full());
1153
1154        assert_eq!(perm.len(), 16);
1155        // Must be a valid permutation
1156        let mut sorted = perm.clone();
1157        sorted.sort();
1158        let expected: Vec<u32> = (0..16).collect();
1159        assert_eq!(sorted, expected);
1160    }
1161
1162    /// Depth-capped BP: with min_partition_docs above the cluster size, only
1163    /// the top-level split happens — clusters still separate (coarse
1164    /// clustering), the permutation stays valid, and the pass converges
1165    /// (a depth cap is a chosen target, not an interruption).
1166    #[test]
1167    fn test_bp_depth_cap_separates_clusters_and_converges() {
1168        // Mostly-separated clusters with one misplaced doc per half — the
1169        // top-level swap pass must exchange docs 3 and 4.
1170        let fwd = make_fwd(
1171            &[
1172                &[0, 1],
1173                &[0, 1],
1174                &[0, 1],
1175                &[2, 3],
1176                &[0, 1],
1177                &[2, 3],
1178                &[2, 3],
1179                &[2, 3],
1180            ],
1181            4,
1182        );
1183        let budget = BpBudget {
1184            min_partition_docs: Some(4),
1185            time_budget: None,
1186        };
1187        let (perm, converged) = graph_bisection(&fwd, 2, 20, budget);
1188        assert!(converged, "depth cap must report converged");
1189        assert_eq!(perm.len(), 8);
1190        let mut sorted = perm.clone();
1191        sorted.sort();
1192        assert_eq!(
1193            sorted,
1194            (0..8).collect::<Vec<u32>>(),
1195            "must stay a valid permutation"
1196        );
1197        // Top-level split separates the clusters (docs {0,1,2,4} share terms
1198        // 0/1; docs {3,5,6,7} share terms 2/3)
1199        let cluster_a = [0u32, 1, 2, 4];
1200        let a_in_left = perm[..4].iter().filter(|d| cluster_a.contains(d)).count();
1201        assert!(
1202            a_in_left == 4 || a_in_left == 0,
1203            "clusters should separate at the top level: {:?}",
1204            perm
1205        );
1206    }
1207
1208    /// Zero wall-clock budget: the pass ends immediately, reports
1209    /// converged=false, and still emits a valid (identity) permutation.
1210    #[test]
1211    fn test_bp_zero_time_budget_emits_valid_partial_permutation() {
1212        let docs: Vec<Vec<u32>> = (0..64).map(|i| vec![i % 4]).collect();
1213        let doc_refs: Vec<&[u32]> = docs.iter().map(|v| v.as_slice()).collect();
1214        let fwd = make_fwd(&doc_refs, 4);
1215        let budget = BpBudget {
1216            min_partition_docs: None,
1217            time_budget: Some(std::time::Duration::ZERO),
1218        };
1219        let (perm, converged) = graph_bisection(&fwd, 4, 20, budget);
1220        assert!(!converged, "zero budget must report unconverged");
1221        assert_eq!(perm.len(), 64);
1222        let mut sorted = perm.clone();
1223        sorted.sort();
1224        assert_eq!(sorted, (0..64).collect::<Vec<u32>>());
1225    }
1226
1227    #[test]
1228    fn test_fast_log2() {
1229        let table = build_log_table(4096);
1230        assert!((table[1] - 0.0).abs() < 0.001);
1231        assert!((table[2] - 1.0).abs() < 0.001);
1232        assert!((table[4] - 2.0).abs() < 0.001);
1233        assert!((table[1024] - 10.0).abs() < 0.001);
1234        // Fallback for values beyond table
1235        let val = fast_log2_lookup(8192, &table);
1236        assert!((val - 13.0).abs() < 0.001);
1237    }
1238}