Jump Index 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)
• phi using lg(h); cool trick by Jounan(?)
• 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
• Michael Sansung: papers on how r changes with point mutations
• 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.
• Idea to reduce memory usage: store LCP lengths separately so that large values don’t take all space.
• 2nd idea by Nicola: only store the LCP when needed at all.
• 3rd idea by Ruben: verify all length \(2^i\) prefixes of the pattern.

11.1 Memory usage with into_par_iter for merging Link to heading

• EF size: 20.3 GB
• during partitioning: 40-41 GB
• EF size: 21.2 GB
• merging: 45..46..47..48..49..52..54..5541
• EF size: 13.3 GB
• extend: 52 .. 52 .. 52
• final EF size: 13.6GB
• memory usage still 52 GB???

11.2 Memory usage with into_iter for merging Link to heading

• collecting: 29 .. 30 .. 40
• EF size: 20.3 GB
• during partitioning: 40-41 GB
• EF size: 21.2 GB
• merging: 41 .. 44 .. 46 .. 47 .. 52
• EF size: 13.3 GB
• extend: 52 .. 52 .. 52
• final EF size: 13.6GB
• memory usage still 48 GB???
• after printing LCP stats: still 48 GB when computing JI stats

11.3 Memory usage with drain fix?? Link to heading

• collecting: 29 .. 30 .. 40
• EF size: 20.3 GB
• during partitioning: 40-41 GB
• EF size: 21.2 GB
• merging: 41 .. 44 .. 46 .. 47 .. 52
• EF size: 13.3 GB
• extend: 52 .. 52 .. 52
• final EF size: 13.6GB
• memory usage still 48 GB???
• after printing LCP stats: still 48 GB when computing JI stats

THE ISSUE WAS THE ALLOCATOR NOT FREEING MEMORY FROM SMALL ALLOCATIONS.

11.4 LCP stats Link to heading

• max 3111690 => 22 bits
• average: 17 => 5+2 = 7 bits with EF
• Number of (pos,c) with >0 link: 815311808
• Number of (pos,c) with >1 link: 275454818 (33.7%)
• Number of (pos,c) with >2 link: 59460310
• Number of (pos,c) with >3 link: 13060635
• Number of (pos,c) with >4 link: 2907064

11.5 2-bit alphabet Link to heading

• 1 EF on (src, c, lcp, target)

  final EF size: 12.538 GB : 1697342725 * (2 + 56) bits

11.6 Compressing LCP values Link to heading

• 1 EF on (src, c, target)

  compact EF size: 6.250 GB : 1345582601 * (2 + 35) bits

• 1 EF on (src, c, lcp, target) with additional targets for each c

  LCP EF size: 2.687 GB : 351760124 * (2 + 59) bits

• Total: 8.937 GB

11.7 Dropping LCP values Link to heading

• 1 EF on (src, c, target)

  compact EF without LCP: 7.816 GB : 1697342725 * (2 + 34) bits

11.8 MPHF for (source, c) Link to heading

• PHF mapping (source, c) to idx phf: 0.449 GB
• Bitvector mapping index to offset bits: 0.264 GB
• List of target values targets: 6.899 GB
• Total: 7.612 GB

11.9 Highly repetitive 8-bit einstein.de.txt: Link to heading

• EF on (s,c,l,t): 24.2 MB
• EF on (s,c,t) + EF on (s,c,l,t): 12.8 + 5.1 = 17.9 MB
• EF on (s,c,t): 16.0 MB
• PHF + bits + ts: 0.8 + 0.4 + 11.9 = 13.2 MB

TODO 12 Link to heading

• Drop RMQs; find them on-the-fly during DFS
• Try separately encoding STPD (jump targets) and JI
• Store only LCP when there are >1 links in a state; might give 4x saving for that part.
• Dynamic push-only EF structure that auto-adapts?

Turns out the issue is that allocators don’t always return all memory…

NEXT UP:

• Compress links by only storing ‘double’ LCP values option 1:

  • EF<SmallLink>: source, char, target
  • Ef<LargeLink>: source, char, LCP, target

  option 2:

  • EF<source>
  • EF<SmallLink>: source_idx, char, target
  • Ef<LargeLink>: source_idx, char, LCP, target

• Store STPD and LPlF array for suffix links?

13 Pseudo matching lengths Link to heading

I did some benchmark for Pseudo Matching Lengths as in the Mov2 paper, which if I understand correctly comes down to: extend from the start of the pattern as much as possible until there is no match anymore. Then repeatedly continue that from the first character that wasn’t matched yet. So this is like some kind of greedy matches but not necessarily MEMs because we start somewhere in the middle of each MEM.

Movi2 reports 4000 reads/s PML throughput for simulated (pacbio?) reads of average length 2500, so 10M bases/s, or 100ns/base, which indeed is close to a cache miss. I’m not sure what the default error rate is here.

I’m taking random substrings of length uniform in 1 to 5000 and adding 1% or 0.1% random mutations. Then I get the following, where ‘part len’ is the average length of the pseudo matches and ‘jump dist’ is the average distance between consecutive jumps.

1
2
3
4
5
6
7
8
9

Influenza:
1.00% errors: Avg   5787 reads/sec. Avg part len:   86.2, avg jump dist:   10.3, avg jumps/part:   8.38
0.10% errors: Avg  39635 reads/sec. Avg part len:  696.3, avg jump dist:   72.7, avg jumps/part:   9.58
E.coli:
1.00% errors: Avg  16360 reads/sec. Avg part len:   83.2, avg jump dist:    9.2, avg jumps/part:   9.07
0.10% errors: Avg  94055 reads/sec. Avg part len:  691.1, avg jump dist:   72.6, avg jumps/part:   9.52
Human genome:
1.00% errors: Avg  14234 reads/sec. Avg part len:   87.3, avg jump dist:    7.4, avg jumps/part:  11.84
0.10% errors: Avg  81313 reads/sec. Avg part len:  696.3, avg jump dist:   57.6, avg jumps/part:  12.09

So with 10x fewer errors everything is much faster, as can be expected since once we anchor on the right location there is no jumping around anymore. That also gives much longer matches and correspondingly a larger distance between jumps. In particular the number of jumps needed to anchor a pseudo match/part is roughly constant at 9 either way. Which probably also means that a prefix lookup table on the first ~10 characters (TODO) will avoid most of those. Also note that the throughput for the human genome is similar to the other cases, even though this is a 30x larger dataset, so this is promising for the extension to full HPRCv2.

14 June 15 Link to heading

14.1 Plan Link to heading

• Extract link computation to separate function
• Store EF-encoded version on disk
• Build and test POS+ version.
  • Bug where 10..2947 of E.coli is not found?
  • Include empty suffix?
• Store POS+ links on disk as well.
• Make SuffixLink type for suffix links, pointing from target to source.
• Store Links and ReverseLink(reverse text) at the same time
• Use suffix links for fast matching statistics
• Dynamic construction

14.2 Testing POS+ Link to heading

1
2
3
4
5
6
7

01234567890123456
AGAGCGAGAGCGCGC#
               ^ root


sa: [0, 6, 2, 8, 15, 14, 4, 12, 10, 5, 1, 7, 13, 3, 11, 9]
Root anchor: 4

Bug in E.coli dataset for POS+:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

Testing mapping
0: mismatch at 112689514: got 1 wanted 2
link: Link { source: 112689514, c: 2, lcp: 3, target: 112689513 }
2: mismatch at 112689515: got end of text wanted 2
link: Link { source: 112689515, c: 2, lcp: 2, target: 112689366 }
4: mismatch at 112689368: got 3 wanted 0
link: Link { source: 112689368, c: 0, lcp: 4, target: 112688049 }
5: mismatch at 112688050: got 2 wanted 1
link: Link { source: 112688050, c: 1, lcp: 5, target: 112685534 }
9: mismatch at 112685538: got 2 wanted 0
link: Link { source: 112685538, c: 0, lcp: 10, target: 110925697 }
11: mismatch at 110925699: got 3 wanted 1
link: Link { source: 110925699, c: 1, lcp: 11, target: 110289776 }
14: mismatch at 110289779: got 2 wanted 3
link: Link { source: 110289779, c: 3, lcp: 14, target: 107623798 }
426: mismatch at 107624210: got 3 wanted 1
link: Link { source: 107624210, c: 1, lcp: 436, target: 37421115 }
462: mismatch at 37421151: got 1 wanted 2
no link found; next is Link { source: 37421258, c: 2, lcp: 14, target: 15827280 }
prev link: Link { source: 37420687, c: 1, lcp: 12, target: 16407000 }

thread 'main' (867347) panicked at src/jump_index.rs:735:13:
substring 10..2947 not found
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

14.3 Failing tests Link to heading

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

01234567890123456789
CCBDDACADDCCADDACBDD
A    *
B *
D  *
CA     *
CB*
   DA*
   DC     *
   DDC    *
     AD *
   DDACB         *
      CADDA    *
CCBDDACADDCCADDACBDD
01234567890123456789

Bug: CCA and CCB are a BWT-run, but the character preceding CCB wraps around the text!

1
2
3
4
5
6

      012345678901234567890123456789
text: DCACCCABDCDCACCCABDCDCCCCCABAC
pat:  **                 CDCCCCC
       C*       *
               CDC*         *
                            ^ MISSING

• missing link for prefix DCACCCABDCDCACCCA => DCACCCABDCDCACCCC
  • had off-by-one in prefix-run check

1
2
3
4
5
6
7

      012345678901234567890123456789
text: DCCACDCBDADCCACDCBDADCAACACBDA
pat:                  CBDADCAACACB
      **
       C*    *
            CBDADC*         *
                            ^ MISSING

• missing link for prefix DCCACDCBDADCC => DCCACDCBDADCA
  • had off-by-one in prefix-run check

TODO 14.4 Support length-0 pattern Link to heading

For POS+, the length-0 pattern currently matches before instead of after the last character.

• TODO: Include the empty suffix in the suffix array.

14.5 Suffix link debugging Link to heading

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

   0123456789012345
T: AGAGCGAGAGCGCGC#   POS- LINKS
   #              *   0->15 (0)
   C   *              0->4  (0)
   G*                 0->1  (0)
   AGC *              2->4  (2)
    GAGA   *          4->8  (3)
      GC#         *   5->15 (2)
   AGAGCGC     *      6->12 (6)
            GCGC# *  13->15 (4)
T: AGAGCGAGAGCGCGC#
   0123456789012345
   #              *  reverse POS+ LINKS
   CGAG*
   G*
     CGCG    *
           CGCGC
             *
  GAG
    *
  GAGAGCG
        *
T: AGAGCGAGAGCGCGC#
   0123456789012345
                     SUFFIX LINKS
   A
   *                 1->0 (0)
   *GAG              2->0 (2)
   *   C             5->0 (0)
    *  CGAG          5->1 (3)
   *    GAGAGC       6->0 (5)
      *      CGCG   11->3 (3)
            *CGCGC  11->9 (4)
   *              # 16->0 (0)
T: AGAGCGAGAGCGCGC#
   0123456789012345

14.6 Match statistics debugging Link to heading

1
2
3
4
5
6
7

T: AGAGCGAGAGCGCGC#
   0123456789012345
P: AGAGCGAGGCG
           *
    GAGG
   AGG
   G

1
2
3
4
5
6

T: AGAGCGAGAGCGCGC#
   0123456789012345
P:  GAGCGCGAGCGCGC
         *
          GAGCGCGAGCGCGC
                 *

14.7 Match statistics results!!! Link to heading

for 100k reads; avg len 2500

Human genome index: 2 * 12.5 = 25 GB (for fwd links and suffix links)

RLZ:

• as fast for influenza
• 2x faster for e.coli
• 3x faster than MS for human genome

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

influenza:
11.5 for 100000 reads of average length 2500 with 1.00% errors. Avg   8672 reads/sec. Avg part len:   27.7, avg jump dist:   19.5, avg jumps/part:   1.42, avg suffix-link dist:   20.9, avg suffix-links/part:   1.33 (jump lookups/pt:   2.75, suffix link lookups/pt:   1.33, suffix link skips/pt:   2.25)
 2.4 for 100000 reads of average length 2500 with 0.10% errors. Avg  40113 reads/sec. Avg part len:  324.1, avg jump dist:  109.7, avg jumps/part:   2.95, avg suffix-link dist:  206.4, avg suffix-links/part:   1.57 (jump lookups/pt:   4.52, suffix link lookups/pt:   1.57, suffix link skips/pt:   6.80)

e.coli:
MS:  9.798 for 100000 reads of average length 2500 with 1.00% errors. Avg  10206 reads/sec. Avg part len:   20.0, avg jump dist:   17.7, avg jumps/part:   1.14, avg suffix-link dist:   18.9, avg suffix-links/part:   1.06 (jump lookups/pt:   2.20, suffix link lookups/pt:   1.06, suffix link skips/pt:   4.89)
MS:  2.070 for 100000 reads of average length 2500 with 0.10% errors. Avg  48298 reads/sec. Avg part len:  234.6, avg jump dist:  114.4, avg jumps/part:   2.05, avg suffix-link dist:  208.6, avg suffix-links/part:   1.12 (jump lookups/pt:   3.18, suffix link lookups/pt:   1.12, suffix link skips/pt:  24.30)

human:
MS:  16.835 for 100000 reads of average length 2500 with 1.00% errors. Avg   5940 reads/sec. Avg part len:   18.1, avg jump dist:   15.7, avg jumps/part:   1.15, avg suffix-link dist:   17.1, avg suffix-links/part:   1.06 (jump lookups/pt:   2.21, suffix link lookups/pt:   1.06, suffix link skips/pt:   7.01)
MS:  2.908 for 100000 reads of average length 2500 with 0.10% errors. Avg  34391 reads/sec. Avg part len:  209.0, avg jump dist:   98.7, avg jumps/part:   2.12, avg suffix-link dist:  194.7, avg suffix-links/part:   1.07 (jump lookups/pt:   3.19, suffix link lookups/pt:   1.07, suffix link skips/pt:  53.61)

15 Dynamic construction Link to heading

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

0123456789012345
AGAGCGAGAGCGCGC#

         1->0 sufl (0) ** init
*G       0->1 link (0)
         2->0 sufl (0) ** wrong len
A
AG
  *
  AGC    2->4 link (2)
C        0->4 link (0)
         5->0 sufl (0) ** wrong len
         5->1 sufl (3) ** MISSING
 G
 GA
 GAG
 GAGA    4->8 link (3)
AGA      6->0 sufl (3) ** wrong len
AGAG
AGAGC
AGAGCG
AGAGCGC  6->12 link (6)
   GC   11->3  sufl (2) ** wrong len
   GCG
          CGCGC ???
         GCGC   11->9  sufl (4) ** MISSING
         GCGC#  13->15 link (4)
   GC#           5->15 link (2) ** MISSING
#                0->15 link (0)
                16->0  sufl (0)  ** OK
AGAGCGAGAGCGCGC#
0123456789012345

     GAGAGCGC
      ------
AGAGCG
   GC      **

          CGCG
   GC      --
           ***
   GCG

          CGCGC
   GCG     ---
           ****
         GCGC

15.1 Take 2 Link to heading

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34

   0123456789012345
T: AGAGCGAGAGCGCGC#  POS- LINKS   SL         not in the pattern
   G*                0->1  (0)    1->0  (0)  v
   AGC *             2->4  (2)    2->0  (2)
   C   *             0->4  (0)    5->0  (0)  *
    GAGA   *         4->8  (3)    5->1  (3)
   AGAGCGC     *     6->12 (6)    6->0  (5)  *
            GCGC# * 13->15 (4)   11->9  (4)
      GC#         *  5->15 (2)   11->3  (3)  *
   #              *  0->15 (0)   16->0  (0)  *
T: AGAGCGAGAGCGCGC#
   0123456789012345  POS- LINKS   SUFFIX LINKS
  A*                 0->1  (0)    1->0  (0)
  G
   A
   AG  *             2->4  (2)    2->0  (2)
  C    *             0->4  (0)    5->0  (0)
    G
    GA
    GAG    *         4->8  (3)    5->1  (3)
   AG
   AGA
   AGAG
   AGAGC
   AGAGCG      *     6->12 (6)    6->0  (5)  *** Dynamic algo gives len 6 here!
    G  *             2->4  (1)
      GC
      GCG      *     6->12 (3)
            GCGC  * 13->15 (4)   11->9  (4)
      GC          *  5->15 (2)   11->3  (3)
  #               *  0->15 (0)   16->0  (0)

T: AGAGCGAGAGCGCGC#
   0123456789012345

Thinks to track:

• in: longest suffix T[pos-len..pos] = T[i-len..i] that occurred at pos<i

• if T[pos] == T[i]:

  • pos+=1
  • len+=1
  • done

• out: longest suffix T[pos'-len'..pos'] = T[i+1-len'..i+1] that occurred at pos'≤i

15.2 More debugging of dynamic construction Link to heading

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

0123456789
BDDAAABDAC
fwd:
    C    *   4->9 (2)  got len 1
suf:
  DA  BDA    7->2 (2)  got len 1

BD    BD
  A     *
| empty suffix link
   A extend with A
  DA -> could have extended directly
BDDAAABDAC
0123456789



T: BDDAAABDAC
P: DAC
FWD1: [0->3 A 0, 0->9 C 0, 0->1 D 0, 2->3 A 1, 2->8 A 2, 4->6 B 1, 4->9 C 2, 5->6 B 2]
FWD2: [0->3 A 0, 0->9 C 0, 0->1 D 0, 2->3 A 1, 2->8 A 2, 4->6 B 1, 4->9 C 1, 5->6 B 2]
SUF1: [1->0 0, 2->0 0, 2->1 1, 4->0 0, 4->3 2, 6->0 2, 7->2 2, 10->0 0]
SUF2: [1->0 0, 2->0 0, 2->1 1, 4->0 0, 4->3 2, 6->0 2, 8->3 1, 10->0 0]

15.3 IDEA: Store one EF for each character?? Link to heading

15.4 It’s working!!! Link to heading

Influenza:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

=== DYNAMIC CONSTRUCTION ===
   1.0Mbp    0.06s 17.80Mbp/s fwd links    0.1M = every 12.8bp suf links    0.1M = every 12.9bp fwd    2.73MB suf    2.71MB
   2.1Mbp    0.14s 14.80Mbp/s fwd links    0.2M = every 11.7bp suf links    0.2M = every 11.7bp fwd    5.98MB suf    5.97MB
   4.2Mbp    0.27s 15.58Mbp/s fwd links    0.3M = every 13.3bp suf links    0.3M = every 13.3bp fwd   10.51MB suf   10.47MB
   8.4Mbp    0.47s 17.86Mbp/s fwd links    0.5M = every 16.5bp suf links    0.5M = every 16.6bp fwd   16.92MB suf   16.86MB
  16.8Mbp    0.79s 21.16Mbp/s fwd links    0.8M = every 21.4bp suf links    0.8M = every 21.5bp fwd   26.13MB suf   26.03MB
  33.6Mbp    1.21s 27.70Mbp/s fwd links    1.1M = every 30.6bp suf links    1.1M = every 30.7bp fwd   36.55MB suf   36.43MB
  67.1Mbp    1.91s 35.08Mbp/s fwd links    1.6M = every 42.9bp suf links    1.6M = every 43.0bp fwd   52.13MB suf   51.98MB
 134.2Mbp    3.52s 38.11Mbp/s fwd links    2.6M = every 51.4bp suf links    2.6M = every 51.6bp fwd   87.05MB suf   86.72MB
=== DYNAMIC CONSTRUCTION IN 4.230 seconds ===

E.coli:

1
2
3
4
5
6
7
8
9

=== DYNAMIC CONSTRUCTION ===
   1.0Mbp    0.69s  1.51Mbp/s fwd links    0.8M = every  1.4bp suf links    0.8M = every  1.4bp fwd   25.09MB suf   25.09MB
   2.1Mbp    1.58s  1.33Mbp/s fwd links    1.5M = every  1.4bp suf links    1.5M = every  1.4bp fwd   50.17MB suf   50.17MB
   4.2Mbp    3.52s  1.19Mbp/s fwd links    3.0M = every  1.4bp suf links    3.0M = every  1.4bp fwd   99.37MB suf   99.37MB
   8.4Mbp    5.77s  1.45Mbp/s fwd links    4.5M = every  1.8bp suf links    4.5M = every  1.8bp fwd  151.21MB suf  151.20MB
  16.8Mbp   11.35s  1.48Mbp/s fwd links    8.1M = every  2.1bp suf links    8.1M = every  2.1bp fwd  270.47MB suf  270.49MB
  33.6Mbp   14.59s  2.30Mbp/s fwd links   10.1M = every  3.3bp suf links   10.1M = every  3.3bp fwd  336.17MB suf  336.16MB
  67.1Mbp   19.96s  3.36Mbp/s fwd links   13.3M = every  5.1bp suf links   13.3M = every  5.1bp fwd  441.83MB suf  441.62MB
=== DYNAMIC CONSTRUCTION IN 24.007 seconds ===

Einstein.de:

1
2
3
4
5
6
7
8
9

=== DYNAMIC CONSTRUCTION ===
   1.0Mbp    0.01s  77.92Mbp/s fwd links    0.0M = every  69.3bp suf links    0.0M = every  69.3bp fwd    0.50MB suf    0.50MB
   2.1Mbp    0.02s  98.83Mbp/s fwd links    0.0M = every 109.3bp suf links    0.0M = every 109.1bp fwd    0.64MB suf    0.64MB
   4.2Mbp    0.03s 123.52Mbp/s fwd links    0.0M = every 182.5bp suf links    0.0M = every 182.3bp fwd    0.77MB suf    0.77MB
   8.4Mbp    0.06s 136.17Mbp/s fwd links    0.0M = every 263.0bp suf links    0.0M = every 261.4bp fwd    1.06MB suf    1.07MB
  16.8Mbp    0.11s 149.39Mbp/s fwd links    0.0M = every 378.7bp suf links    0.0M = every 374.4bp fwd    1.48MB suf    1.49MB
  33.6Mbp    0.22s 155.96Mbp/s fwd links    0.1M = every 490.0bp suf links    0.1M = every 485.0bp fwd    2.28MB suf    2.31MB
  67.1Mbp    0.41s 164.68Mbp/s fwd links    0.1M = every 671.8bp suf links    0.1M = every 663.6bp fwd    3.33MB suf    3.37MB
=== DYNAMIC CONSTRUCTION IN 0.560 seconds ===

15.5 Optimizing the relative store Link to heading

Starting point:

• fwd: EF + BTreeSet
• suf: EF + BTreeSet

Better:

• fwd: EF + HashMap<(source, c): u64, (len: u32, idx: u32)> + Vec<Vec<(lcp, target): u64>>
• suf: EF + Vec<(source, lcp, target)> with binary search or an S+-tree.

References Link to heading