STPD dev notes

These are ongoing notes and logs on improving the STPD (Becker et al. 2025) and a new variant of it that I’ll call the jump index.

1 Results Link to heading

twice the same table; bottom one in millions

short repetitive sequence, 1% error rate, alphabet size 2 and 4

long repetitive sequence, 0.1% error rate, alphabet size 2 and 4

Datasets from pizza&chili repetitive corpus:

• influenza: 78k copies of influenza; probably with very high variance
• dna.001.1: is 100 synthetic copies of a 1 MB string, each with 0.1% independent mutations
• einstein.de.txt: revisions of the German Einstein Wikipedia article
• Escherichia_Coli: 23 copies of E. Coli.

Variables:

• \(n\): size of input
• \(k\): the number of copies
• \(r\): BWT-runs
• \(e = |\mathsf{CDAWG}|\) is the number of edges in the CDAWG and is copied from Cleary et al. (2024).
  • My numbers in the CDAWG-e column are just slightly (<100) off. Maybe because I don’t append a sentinel?
• \(|\mathsf{stpd}|\): STPD size
  • equal to number of CDAWG nodes with a non-path incoming edge?
• src: number of text positions from which we jump
• (src, c): number of (text positions, character) for which there is a jump
• #(links): total number of links in Jump Index.
  • TODO: Why not equal to |CDAWG|/2?

Observations:

• For all but influenza, we have roughly \(r = n/(k/2) = 2\cdot n/k\).
• the STPD size is always 60% to 75% of \(r\)
• Number of links for pos and lex is very close to \(r\)
• Number of links for colex is much worse, but why? Especially when there are many repeats.

2 Faster construction by skipping runs Link to heading

before: 60s

after: 35s

Idea: if an SA interval \([l, r)\) is contained inside a BWT run, we can skip processing it completely.

This results in close to \(4r\) samples in total. Probably because each ‘critical’ interval is processed once for each character somehow?

• After reusing BWT between variants: 29s
• After switching to non-precomputed RMQ: still 29s
• After switching to S=256 instead of S=64: still 29s
• using par_iter for permuted_pi: 22s
• using S=128: 21s
• build RMQ using par_iter for block mins: 20s

New overall results: (67s)

TODO 2.1 We loose 1 or 2 samples in the new setting. Why? Link to heading

TODO 2.2 More fine-grained pruning to directly construct the final result without duplicates? Link to heading

TODO 2.3 Construction direction via run-length BWT? Link to heading

TODO 2.4 PHF on (pos, c) Link to heading

TODO 2.5 compressed targets Link to heading

Instead of storing all jump targets, store a list of all target locations (as in STPD). But now partition it by the last character. Then store one list per final character (sorted by text pos, for lg(n/r) bits each) and the index in the list (lg r bits each). This should help compression when targets are reused.

TODO 2.6 Prefix table up to k=10 Link to heading

TODO 3 Queries Link to heading

• 1-PHF to map \((i,a)\) to its range

TODO 4 suffix links Link to heading

• How good are our suffix links??? How many do we need to traverse to get to the next position on average?

5 Ideas Link to heading

• max LCP length relates to max run length: both increase for more repetitive texts
• pBWT: highlight jump-targets
  • Either binary search for them in \(O(l)\) (\(l = \Theta(n/r)\) copies), or store target in \(\lg l\) bits.
• STPD, say of size \(s\):
  • \(s\) \(\lg n\) bit samples, shuffled
  • can we compress them?
    • typically, we will jump from pos \(i\) in one copy to pos \(\approx i\) in another copy? So then a binary search over \(\ell\) things is enough?
    • Consecutive STPD samples are in similar text positions?
  • Jump targets will usually point to roughly the same position in the text
• Compressing the Jump Index:
  • MPHF for \((i, a)\)
  • LCP \(\geq 255\) can be stored in a secondary structure with 32-bit values, so 8 bits is enough for the primary structure
• Use a prefix lookup for the first ~10 bp = ~20 bits

6 Analysing and reducing space usage Link to heading

6.1 Movi 1 Link to heading

Movi (Zakeri et al. 2024, 2025) is an engineered version of a run-length FM-index. Similar to the r-index, but avoiding all the rank and select structures and instead inlining just the right data so that most LF-mappings can be done with only a single cache miss.

Movi1 uses 16 bytes per run.

6.2 Movi 2 Link to heading

Movi2 splits some runs:

• Runs longer than \(2^\ell\) are split into two (+5% for HPRCv2), so the length can be stored in \(\ell\) bits
• In a run of A, if we instead want to append a C, the best match is in the run of C’s preceding or succeeding the current run. Most of the time, this choice is the same for all elements of the run, and a single bit up/down pointer for each of \(\sigma-1\) of the non-run characters suffices. If this is not the case, the run is split, for 10% extra runs on HPRCv2.

It then stores:

• 2-bit character
• \(\ell=11\)-bit run length
• the \(\log_2(n)\)-bit position where the entire run maps to:
  • a \(\log_2( r)=36\)-bit run-index \(\xi\)
  • an 11-bit index in the run
• \(\sigma-1=3\) bits indicating whether to round up or down for each non-run characters.

For a total of 63 bits or 8 bytes, with the bulk being the \(r\log_2n\) bits for the target indices.

6.2.1 Block-based run-indices Link to heading

The largest component are the run-indices \(\chi\). But these can be compressed further, since \(\xi\) is increasing over all runs for the same character. An elias-fano structure would be space efficient, but a practical solution is to only store \(\xi\) values (for all characters) every \(b\) rows and then linearly interpolate between them, only storing the delta to the prediction. In this mode, \(\ell=10\) and \(\Delta=23\) (I think), for a total of 6 bytes.

Instead of storing \(\Delta\), a further down-sampling followed by a linear scan is also possible, bringing the space per run to just 3 bytes. (I do not yet fully understand the details here.)

6.3 STPD Link to heading

Suppose we have \(k\) similar sequences of total length \(n=k\cdot m\). In practice, the STPD has roughly \(s=\frac 23 r\) samples (see tables above).

Each sample is a text position, for a total space of \(s\log_2 n \approx r\log_2 n\), i.e., similar to the move-index. In fact, 6 bytes (48 bits) per sample is plenty, matching the block-based movi2 (but not the sampled version).

Do note though that this excludes the relative-Lempel-Ziv compressed text, as well as a 2-4 bit/element RMQ structure that selects the best match in the interval of the STPD (sparse prefix array) that matches the pattern prefix:

• For colex-, this cam be omitted and we take the start of the interval instead.
• For pos-, we could also iterate over the entire matching interval. Not sure if that would be slow.
• For lex-, the RMQ structure is needed.

Unlike for movi, where the target runs are increasing and can thus be compressed, here we must sort the \(s\) integers in colex order, meaning that (a priory) we cannot compress them via something like Elias-Fano coding.

But there is hope!

6.4 Compressing STPD: STPD for pBWT Link to heading

Consider the setting where the \(k\) sequences only differ by say 1% of substitutions, and nothing else. Then we can put them below each other (as in the pBWT), and for each column where the pattern matches, we could report the first input sequence in which it matches.

Thus, we get something like this:

1
2
3
4
5

Pat: ABCDE
T1   ABCDEFGHIJX
T2   ABCDEFGHIJY
T3   ABCDEZGHIJX
T4   ABCDEZGHIXX

Here, the pattern matches a prefix of \(T_1\). Then, if we extend the pattern by Z, it does not match \(T_1\) anymore, and instead we have to jump to \(T_3\). Thus, the Z in \(T_3\) would be an STPD-sample (in the concatenation of the texts).

Thus, in every column where “things happen” we get one or more STPD samples. But in this setting, we don’t have to binary search over all samples in \(O(\log s)\)! Instead, we only have to binary search the samples in the current column, which is more like \(O(\log s/m)\), or either way at most \(O(\log k)\).

Thus, we could ideally use some perfect hash function to map active columns (with at least one sample) to a variable-length list of samples in around 2 bits per active column. Or we could store in \(\log_2 k\) bits the number of STPD samples in each active column, and then distribute columns over multiple PHFs based on the number of samples.

Then, all that’s needed is to store the subset of the \(k\) sequences that has an STPD sample in this column, which is only \(\log_2 k\) bits (around 9 for HPRCv2) rather than \(\log_2 n\) (around 44).

Thus, this scheme needs only \(m\log k + 2m +s\log_2 k\) bits.

A big restriction, though, is that this does not allow jumping between different columns: once the pattern in anchored to some column, it can only switch between rows.

Nevertheless, it feels like this should provide some possible ways to optimize the generic STPD as well: most of the time, jumping happens between ’equivalent positions’ in different genomes, and if we have some course multiple-sequence-alignment or so between them, then it might be sufficient to only store low bits most of the time to indicate ‘jump to the equivalent location in the input sequence with ID \(x\)’.

6.5 Jump Index Link to heading

• HPRCv2: 466 copies of Gbp => 40.5 bits of text
• At 8 bytes/run and including rc and 15% extra runs, movi2 has 47.3 GB.
  • 5.1 G runs
  • 2.55 G runs for 466 copies of fwd-only
  • 1.72 G runs for a single copy of a human genome

Initial space we need per pointer:

• source: 41 bits
  • hide in MPHF on (pos, c) => 2 bits
  • EF to store all lengths?
  • Or EF for source locations + EF for lengths?
  • <1 byte total?
• target: 41 bits = 5 bytes
• LCP len: 16 bits (probably 8 + fallback) = 2 bytes

6.5.1 Compact link Link to heading

Current

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

pub struct Link {
    /// 5 bytes?
    pub source: usize,
    /// 2 bits?
    pub c: char,
    /// 2 bytes?
    pub lcp: usize,
    /// 5 bytes?
    pub target: usize,
}

• naive: 24 bytes
• packed: 100 bits = 12.5 bytes
  • 41 bit src
  • 2 bit char
  • 16 bit LCP
  • 41 bit target
• EF: we have ~2G runs, so we encode the top 31 bits with only 2 bits => 71 bits each

6.5.2 Mutable Link to heading

• Giulio’s mutable rank/select B-tree on the high part of the EF
  • maybe 20% overhead; that’s all fine
  • Does it support insertions though… It’ll need to be modified for this
• B-tree on the low parts of the EF values, so insertion is fast
  • Needs special support for indexing
  • Split so that low parts are 64 bits??? 9-bit LCP?

TODO 6.6 Link to heading

• Number of edges in graph on STPD samples?
  • (this is larger than just the number of links)
  • equivalent to CDAWG? Smaller?
• Number of suffix links?

7 CDAWG Link to heading

What is the relation between suffixient sets/STPD and the CDAWG?

• STPD is a path decomposition of the CDAWG
  • Size is the number of nodes with >1 incoming edge, ie at least one incoming edge not in the path decomposition. Up to 10x smaller than the CDAWG, see table at the top.
• Jump Index has size equal to the number of non-path edges in the CDAWG?
  • In the table, the size of the jump index varies between 20% and 50% of the CDAWG. I’m not yet sure why it’s not always at least 50%, ie where at most half the edges are heavy-path and the other (at least) half light edges needing a jump pointer.
  • Hmm; the JI jumps to a different location depending on the length of the substring already matched, whereas the CDAWG is independent of the already-matched prefix and so is more like COLEX-? Or is this not an issue because such cases have different CDAWG nodes?
• CDAWG is bad (\(O(n)\)) for AAAAAAAA#. But maybe fine in practice?
• CDAWG nodes are maximal repeats \(S\), and have an equivalence class of right-maximal strings \(\{S[0..), S[1..), S[2..), \dots, S[k..)\}\) such that all of those (but the first) are not left-maximal, but \(S[k+1\dots)\) is left-maximal. (Belazzougui and Cunial 2017)
  • All strings in the equivalence class have the same range-length in the BWT and for \(i\geq 1\), all \(S[i..)\) the range is (contained in) a single run with character \(S[i-1]\). OTOH, the BWT interval for \(S[0..)\) contains at least two characters.
  • All paths to a node in the CDAWG are a suffix of one another.
• The in-neighbours of a node exactly partition the range of ‘active lengths’.

CDAWG:

• A node for every maximal repeat, that is: a substring \(S\) that is repeated and can be extended both left and right in at least 2 ways. Or equivalently, \(aS\) and \(Sa\) occur stricly less than \(S\) for all symbols \(a\).
• Space: number of right-extensions of maximal repeats

“Composite Repetition-Aware Data Structures” (Belazzougui et al. 2015)

• (google scholar citations to this paper give a good source of papers)
• introduces \(e\), the number of arcs in the CDAWG
• Index combining run-length BWT with CDAWG

“Representing the Suffix Tree with the Cdawg” (Belazzougui and Cunial 2017)

• Heavy-path decomposition of CDAWG; supports additional suffix-tree operations

7.1 Further papers Link to heading

“Suffix Trees, Dawgs and Cdawgs for Forward and Backward Tries” (Inenaga 2020)

• Takes as input not a string, but a suffix tree or so.

“The Cdawg Index and Pattern Matching on Grammar-Compressed Strings” (Cleary et al. 2024)

• Builds a CDAWG on grammars; 4x smaller than the original CDAWG

“Indexing Highly Repetitive String Collections, Part I: Repetitiveness Measures” (Navarro 2021b)

“Indexing Highly Repetitive String Collections, Part Ii: Compressed Indexes” (Navarro 2021a)

8 Fast construction Link to heading

• Goal: \(O( r)\) space; \(O(n)\) is wasteful
  • BWT+LCP in \(O( r)\) is possible (?)
• Skip constructing SA on repetitive parts?
  • Simply dynamically add them and don’t process them at all?
  • For each mutation find the longest extension to left/right to make it unique. Then only process that window?

9 HPRCv2 Link to heading

Get the AGC (assembled genomes compressor) version as described here (3.3 GB .agc):

1

wget https://s3-us-west-2.amazonaws.com/human-pangenomics/submissions/B4174A5F-F20E-4DCF-8470-F8A907B640BC--HPRCv2_0.6.1_pr_agc_submission/HPRC_r2_assemblies_0.6.1.agc

Or get all of them as alignments against a reference from here (5.2 GB .paf.gz).

Rust AGC api: https://docs.rs/ragc-core/0.1.1/ragc_core/

10 Questions Link to heading

• Why is the number of jumps so close to \(r\)
• data layout
  • dynamic vs static implementation?
• lg(k) instead of lg(n) bits per jump?
• suffix links?

10.1 Notes Link to heading

• BML: faster MEM finding
• moni align
• phi (and inverse) also give LCP values
• classic: find all occurrences of all MEMs, but not what is needed in practice
  • 100 matches to 3 locations is not super useful
  • Only care about distinct locations in pangenome graphs
    • Paper Joni + Paten, WABI 2025: tag array
    • only store distinct mapping locations for runs
• Latin 2025 paper Gonzalo, Joni, Travis
• Jannick Olbrich: MSA
  • sample runs; how many columns you move left after LF mapping
• Heng:
  • vg-surject: maps alignment to graph back to standard reference
  • tag array of vg-surject values
• Goal: mix tag arrays with suffixient sets
  • Problem: these don’t give BWT interval
  • Can only enumerate the interval; not “RMQ” query it
  • not count query, but interval query
• lg(n) vs lg(h)
• movi: r lg n, but r lg(n/r)
  • move structure: row-major
  • compression: col-major
• paper from Roberto’s festschrift
• SMEM: MEM aligned by pBWT columns
• suffixient set for findings MEMs; then r-index to locate all occurrences?
• SPIRE paper:
  • Heng: grlBWT; LMS grammar
  • Ben: PFP parsing
• z-fast trie vs centroid decomposition
• WABI 22 ??

Demo:

• patterns up to length 20
  • obvious that most stuff is useless

10.2 Jump Index Link to heading

• dynamic construction
  • quickly discard existing regions
• dynamic predecessor structure is all that’s needed
  • log(u/n) ideally
• similar to ukkonen’s SA construction
• Store lg(h) instead of lg(n), relative to MSA
  • all relative to reference?
  • break each pointer into 2 pieces:
    • this genome lg(h)
    • roughly this position lg(w)
  • if mapping is monotone; use checkpoints
• Can enumerate all occurrences via DFS
• Associate TAG array information with jump links?
  • But we still
• Can we do counts?
• Draw as automaton
  • What happens when we add a copy?
  • How many states do we gain?

11 Shrinking practical memory usage Link to heading

• SA: 24GB = 8B * 3G, because libsais needs i64 for \(2^{31.5}\) elements
• PLCP: 24GB
• Shrink SA to u32 for 12GB
• Shrink PLCP to u16 for 6GB
• Build LCP, also 6GB
• Peak so far: 48GB

New:

• Text: 1 B -> 3 GB
• BWT: 1 B -> 3 GB
• SA: 32 bit = 4 B -> 12 GB
• LCP: 32 bit = 4 B -> 12 GB
• Dropped range boundaries
• RMQ via segtree instead of sparse table
• SA RMQ: 1.5 b -> 500 MB
• LCP RMQ: 1.5 b -> 500 MB

Total: 31GB

• 32 byte links: ~23GB for ~1/3rd of the data

• 16 byte links: ~24GB for 560/750 intervals = 3/4th of the data

• Replace LCP by Compact PLCP version

  • 12.4 GB => 0.8 GB

• 20 GB now: 3+3+12+0.5+0.5+0.8

  • 40 GB left for links :)

• collecting links is SLOW because of RMQ requiring up to 128 select queries…

  • 36 GB before dedup; ~56 GB total
  • 27 GB after dedup; ~47 GB total

• POS- for human genome in 9.5min wall time (576s) on 6 cores / 12 threads.

Max LCP: 3111690 => 22 bits

• links can fit in 12 bytes, but that makes sorting them much much slower.
• EF compresses links from 27 GB to 13 GB!
  • MPHF would make it even smaller but need something to counter false positives.

References Link to heading