MRG: fix RocksDB-based gather & other rust-based infelicities revealed by plugins

[ctb]

This PR fixes a bug in disk_revindex.rs::RevIndexOps::gather where RocksDB-based gather was
incorrectly subtracting hashes multiple times from the counter in situations of high redundancy.

For example, consider this Venn diagram of the 3-way intersection between three sketches:

When a metagenome contains the union of all three of these sketches, the broken implementation would subtract the 10 at the center multiple times. This was caused by removing hashes from the matches, rather than the intersection, each pass through the counter.

Of note, this made RocksDB-based fastmultigather return incorrect results, ref sourmash-bio/sourmash_plugin_branchwater#322; first discovered in #3138 (comment).

The PR fixes this, and adds a more complete pair of tests, based on test_gather_metagenome_num_results in the Python tests for sourmash.

This PR also adjusts the hash function string for DNA sketches in Rust to be uppercase DNA rather than lowercase dna, ref sourmash-bio/sourmash_plugin_directsketch#49

And remember, it's not just the destination - it's the friends you make along the way, like env_logger.

• Fixes RocksDB-based RevIndex gather functionality is broken. #3139
• Fixes DNA moltype output in manifest is dna not DNA sourmash_plugin_directsketch#49

For consideration:

Right now we are calculating the intersection twice, once in disk_revindex.rs and once in calculate_gather_stats in revindex/mod.rs. This is unnecessary, of course. But the function signature for calculate_gather_stats would need to change to take the intersection as an argument. We could:

• keep calculating it twice, just for simplicity;
• change calculate_gather_stats to take an optional intersection, and calculate it if not provided;
• change calculate_gather_stats to require an intersection.

TODO:

• add at least one explicit test for the moltype fix

Other notes:

• there is a discrepancy between the Python (sourmash gather) and Rust (sourmash_plugin_branchwater results) calculations for remaining_bp. It seems to me like the Python one is definitely wrong; not yet sure about Rust. Viz gather is calculating remaining_bp incorrectly #3194.

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 86.69%. Comparing base (a133e68) to head (231a134).
Report is 78 commits behind head on latest.

Additional details and impacted files

@@            Coverage Diff             @@
##           latest    #3193      +/-   ##
==========================================
+ Coverage   86.67%   86.69%   +0.01%     
==========================================
  Files         136      136              
  Lines       15832    15840       +8     
  Branches     2716     2716              
==========================================
+ Hits        13722    13732      +10     
+ Misses       1800     1798       -2     
  Partials      310      310              

Flag	Coverage Δ
hypothesis-py	25.35% <ø> (ø)
python	92.33% <ø> (ø)
rust	62.26% <100.00%> (+0.16%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[ctb]

The situation so far:

from sourmash-bio/sourmash_plugin_branchwater#322 (comment),

in #3138 (comment), commit 10d5ee8, I add an analog to the Python test test_gather_metagenome. This demonstrates what is hopefully 😭 the same problem - we're getting only 6 matches in the Rust code (instead of 11), and the 6th match starts to diverge from the values we see in the Python implementation.

I've narrowed this down further to running gather on combined.sig against three sketches, the ones for NC_000853.1, NC_011978.1, and NC_009486.1:

signature: NC_009486.1 Thermotoga petrophila RKU-1, complete genome
signature: NC_011978.1 Thermotoga neapolitana DSM 4359, complete genome
signature: NC_000853.1 Thermotoga maritima MSB8 chromosome, complete genome

These three have f_match column with fastmulgather+rocksdb that is broken:

vs

the correct (I'm assuming) set of matches from sourmash gather:

note the last entry, .484 vs .431.

These are the first three matches that have a mutual overlap among all three:

so my spidey sense is tingling about that, but I don't have a clear idea of exactly how it would result in what I'm seeing yet ;).

[ctb]

And I've confirmed with simple Python code that the f_match numbers returned by sourmash gather are correct, at least to my understanding of the gather algorithm 😆

import sourmash

combined_sig = list(sourmash.load_file_as_signatures('combined.sig',
                                                     ksize=21))[0]
combined_mh = combined_sig.minhash
combined_set = set(combined_mh.hashes)

match1_sig = list(sourmash.load_file_as_signatures('match1_NC_000853.1.sig',
                                                   ksize=21))[0]
match1_mh = match1_sig.minhash
match1_set = set(match1_mh.hashes)

match2_sig = list(sourmash.load_file_as_signatures('match2_NC_011978.1.sig',
                                                   ksize=21))[0]

match2_mh = match2_sig.minhash
match2_set = set(match2_mh.hashes)

match5_sig = list(sourmash.load_file_as_signatures('match5_NC_009486.1.sig',
                                                   ksize=21))[0]
match5_mh = match5_sig.minhash
match5_set = set(match5_mh.hashes)

f_match_1 = len(match1_set.intersection(combined_set)) / len(match1_set)
print('1', f_match_1)

combined_set -= match1_set

f_match_2 = len(match2_set.intersection(combined_set)) / len(match2_set)
print('2', f_match_2)

combined_set -= match2_set

f_match_5 = len(match5_set.intersection(combined_set)) / len(match5_set)
print('5', f_match_5)

yields:

1 1.0
2 0.898936170212766
5 0.4842105263157895

which matches the sourmash gather f_match column. Ausgezeichnet!

[ctb]

EUREKA!

[ctb]

tl;dr it looks like this line:

sourmash/src/core/src/index/revindex/disk_revindex.rs

Line 374 in a133e68

match_mh

is double-subtracting hashes that were already removed previously, because they are in the (complete) match but no longer in the intersection with the reduced query.

If you look at the Venn diagram above, there are 10 hashes that are at the intersection of all three - these 10 hashes are being removed from the counter EACH time through the loop.

[ctb]

This is ready for review @luizirber @bluegenes. Note duplicate calculation of intersection; thoughts welcome.

[ctb]

oh, and should we bump the sourmash-rs core version as part of this PR?

[luizirber]

👏

And +1 on doing a r0.14.0 release

[luizirber]

Nice! And with this change there is a block below:

counter.entry(dataset).and_modify(|e| {
    if *e > 0 {
        *e -= 1
    }
});

that can be

counter.entry(dataset).and_modify(|e| { *e -= 1 });

instead, and will panic if e goes below 0.

With this change the error triggers in latest in the index::revindex::test::revindex_load_and_gather_2 test, and doesn't trigger in this PR

[ctb]

added in 231a134, thanks!

[ctb]

And +1 on doing a r0.14.0 release

version bumped in 231a134. I'll merge this once tests pass, but wait a day or two to cut an actual release, in case @bluegenes has comments.