Odd behavior in external benchmarks

[feldroop]

I need to figure out all of the things that seem to be going wrong here: RagnarGrootKoerkamp/quadrank#1

[RagnarGrootKoerkamp]

Pushed support for MemSize and MemDbg to my fork https://github.com/RagnarGrootKoerkamp/genedex.

In the smallest configuration (CondensedTextWithRankSupport<i32, Block512>) and with alphabet::ascii_dna for a text of length 537M, this prints the below.
i.e.: 200 of 214MB is in interleaved_blocks.

The minimal required space is 537/4 = 135MB, so you're at 58% overhead, or 3.19bits per base.

214.5 MB 100.00% ⏺: genedex::FmIndex<i32, genedex::text_with_rank_support::condensed::CondensedTextWithRankSupport<i32, genedex::text_with_rank_support::block::Block512>>
  316  B   0.00% ├╴alphabet: genedex::alphabet::Alphabet
  280  B   0.00% │ ├╴io_to_dense_representation_table: alloc::vec::Vec<u8>
   28  B   0.00% │ ├╴dense_to_io_representation_table: alloc::vec::Vec<u8>
    8  B   0.00% │ ╰╴num_io_symbols_not_searcheable: usize
   72  B   0.00% ├╴count: alloc::vec::Vec<usize>
212.4 MB  99.02% ├╴text_with_rank_support: genedex::text_with_rank_support::condensed::CondensedTextWithRankSupport<i32, genedex::text_with_rank_support::block::Block512>
201.7 MB  94.05% │ ├╴interleaved_blocks: alloc::vec::Vec<genedex::text_with_rank_support::block::Block512>
10.51 MB   4.90% │ ├╴interleaved_block_offsets: alloc::vec::Vec<u16>
164.2 kB   0.08% │ ├╴interleaved_superblock_offsets: alloc::vec::Vec<i32>
    8  B   0.00% │ ├╴text_len: usize
    8  B   0.00% │ ╰╴alphabet_size: usize
2.102 MB   0.98% ├╴suffix_array: genedex::sampled_suffix_array::SampledSuffixArray<i32>
2.101 MB   0.98% │ ├╴suffix_array_data: alloc::vec::Vec<u32>
  100  B   0.00% │ ├╴text_border_lookup: std::collections::hash::map::HashMap<usize, i32>
    8  B   0.00% │ ├╴sampling_rate: usize
    0  B   0.00% │ ╰╴_compression_marker: core::marker::PhantomData<i32>
   88  B   0.00% ├╴text_ids: genedex::text_id_search_tree::TexdIdSearchTree
   32  B   0.00% │ ├╴nodes: alloc::vec::Vec<genedex::text_id_search_tree::Node>
   56  B   0.00% │ ╰╴sentinel_indices: alloc::vec::Vec<usize>
  200  B   0.00% ╰╴lookup_tables: genedex::lookup_table::LookupTables<i32>
   32  B   0.00%   ├╴factors: alloc::vec::Vec<usize>
  160  B   0.00%   ├╴tables: alloc::vec::Vec<genedex::lookup_table::LookupTable<i32>>
    8  B   0.00%   ╰╴num_symbols: usize

[RagnarGrootKoerkamp]

I think you're pretty good here already. Maybe previous comparisons where with the sampled suffix array included, but I've disabled that since it would otherwise take a lot of space and I have not yet implemented lookup so far, only count.

Also given that you have a num_dense_symbols of 5, 25% overhead (2.5bits/symbol) is expected, and you're currently another ~25% above that, which is quite allright I'd say. (I'm at 33% above the baseline at 2.67bits/bp, but I did make some effort to avoid a size-5 alphabet.)

Also note that in the readme, you're at 5.2GB for a human genome of ~3.2Gbp or 800MB packed. That comes out at 6.2bits/bp, which is a bit higher, but most of that is probably the sampled suffix array.

[feldroop]

Wow, what an awesome tool!

The condensed layout is a reordering of storing the fixed-size bit-compressed chars in a bitarray + rank support (interleaved (super)block offsets). Since 537 * 3 / 8 = 201, the size of the interleaved blocks is as expected.

I need to compare this with the other implementations. For condensed 64, only the rank support should be larger. For the flat layout, the memory usage of the main array should change to 537 * 5 / 8 = 336.

[feldroop]

Hey @RagnarGrootKoerkamp,

I was wondering, how did you manage to get rid of the sentinel? (For locate it's pretty straightforward, but I'm unsure about count)

Even if we leave multiple texts/chromosomes aside, there still is the issue with matches wrapping around the end of the text.

The only idea I had is if we assume a maximum query length of k, then we could store the suffix array positions of the last k suffixes and check if the final interval of count contains such a position.

But I don't see a way without this assumption.

[RagnarGrootKoerkamp]

Ah good point. I'm not super familiar with the details and I guess I kinda cheated my way through this 😅.

The way that BWA solves this is by explicitly storing the index of the $ in the BWT, and then adjusting interval endpoints by +-1 when they are after that position (or something along those lines).

Another option that I kinda like but is also iffy is to just pad the sequence with some random data, or just many As, as you're unlikely to have that as an interesting query.

[feldroop]

Ah, that's super interesting! So you don't store the sentinel in the occ table, but post-process the interval after each step to make it equivalent to what you'd get if the sentinel was there?

I think that could be generalized for m texts using one or two binary searches in an array of length m in each step. Do you know if BWA does that, too? I'm a little scared of the performance impact, but I might try it at some point. Smaller≈faster, after all 😄.

Also, what do you think about handling N's and other degenerate IUPAC chars that can be found in many reference genomes? Those were my original motivation for not caring about the sentinel, even though the dream of a 2 bit DNA alphabet is so sweet.

The padding approach might be the easiest and most efficient, but I agree it's iffy. My goal with this library is to have a general purpose implementation of the FM-Index, which is robust and provides all necessary features for typical applications (potentially even outside of bioinformatics). So I think the padding approach is not the way I want to go.

[RagnarGrootKoerkamp]

Ah, that's super interesting! So you don't store the sentinel in the occ table, but post-process the interval after each step to make it equivalent to what you'd get if the sentinel was there?

Exactly. See here for how BWA does this:
https://github.com/lh3/bwa/blob/b92993c1161e73167181558856567ef2f367e3f0/bwt.c#L114

I think that could be generalized for m texts using one or two binary searches in an array of length m in each step. Do you know if BWA does that, too? I'm a little scared of the performance impact, but I might try it at some point. Smaller≈faster, after all 😄.

Hmm yeah; this might require some custom solution, but you'd want to avoid that. A hacky solution again would be to do $AAAAAA or $<fixed pseudo random string>, so that all $ are clustered in the same position in the BWT.

Also, what do you think about handling N's and other degenerate IUPAC chars that can be found in many reference genomes? Those were my original motivation for not caring about the sentinel, even though the dream of a 2 bit DNA alphabet is so sweet.

Personally I don't really believe in IUPAC :P I think for FM-index specifically, it makes more sense to just random/deterministically project them to ACGT. Either way you're probably going to query an FM-index while allowing some small number of errors.

The padding approach might be the easiest and most efficient, but I agree it's iffy. My goal with this library is to have a general purpose implementation of the FM-Index, which is robust and provides all necessary features for typical applications (potentially even outside of bioinformatics). So I think the padding approach is not the way I want to go.

Yeah, that's very fair. I agree robustness should be very high on the list of design principles.