[MRG] re-implement the actual gather protocol with a cleaner interface.

[ctb]

Note: This is a PR into #1370; review & merge can wait until that PR is merged.

This PR trials an actual gather protocol at the class level, as a potential alternate fix to #1263. See also comments here and above.

The core code is located here, which provides a CounterGather class that collects and prioritizes matches for gather and also supports cross-database gather.

The key piece of the refactoring is that CounterGather now provides two methods, peek(query) and consume(...). peek(query) provides the best containment result from this counter, but does not adjust any of the internal information; consume(...) is used to remove a match (potentially from a different CounterGather object).

Below is some code implementing multi-database gather that (when shoehorned into the current gather code :) passes all of the tests:

    def gather(self, query, threshold_bp):
        results = []

        best_result = None
        best_intersect_mh = None

        # find the best score across multiple counters, without consuming       
        for counter in self.counters:
            result = counter.peek(query.minhash, query.minhash.scaled,
                                  threshold_bp)
            if result:
                (sr, intersect_mh) = result

                if best_result is None or sr.score > best_result.score:
                    best_result = sr
                    best_intersect_mh = intersect_mh

        if best_result:
            # remove the best result from each counter                          
            for counter in self.counters:
                counter.consume(best_intersect_mh)

            # and done!                                                         
            return [best_result]
        return []

TODO:

• test flatten and downsampling
• test consume in more detail, maybe

[codecov]

Codecov Report

Merging #1489 (af52c41) into add/prefetch_cli (903239b) will increase coverage by 0.23%.
The diff coverage is 100.00%.

@@                 Coverage Diff                  @@
##           add/prefetch_cli    #1489      +/-   ##
====================================================
+ Coverage             94.89%   95.12%   +0.23%     
====================================================
  Files                    99       99              
  Lines                 16662    17044     +382     
  Branches               1545     1562      +17     
====================================================
+ Hits                  15811    16213     +402     
+ Misses                  616      603      -13     
+ Partials                235      228       -7     

Flag	Coverage Δ
python	95.12% <100.00%> (+0.23%)	⬆️

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

Impacted Files	Coverage Δ
tests/test_sourmash.py	99.72% <ø> (-0.01%)	⬇️
src/sourmash/commands.py	85.21% <100.00%> (+0.12%)	⬆️
src/sourmash/index.py	95.22% <100.00%> (+4.43%)	⬆️
src/sourmash/search.py	96.42% <100.00%> (+0.48%)	⬆️
tests/test_index.py	100.00% <100.00%> (ø)
src/sourmash/minhash.py	93.33% <0.00%> (+0.63%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 903239b...af52c41. Read the comment docs.

[ctb]

aaaaaand just because, I implemented Index.counter_gather(...) and completely swapped out the gather implementation in favor of the new approach.

As a bonus (?), I may have discovered a bug in multigather.

[ctb]

Note: there are some potential cleanups and refactorings in gather_databases but I'm going to rest on my laurels at this point, and await feedback before continuing.

[ctb]

Some benchmarking for y'all using @bluegenes new zipfile databases -- running on farm head node 😏

These use the code above, so prefetch and CounterGather and (because zipfiles) --linear.

tl;dr:

• searching 35k GTDB reps is 4 minutes and 85 MB.
• searching 300k GTDB all is 24 minutes and 4 GB.
• it's almost all user space, not I/O (presumably zip is doing some caching here?)

against ~35k GTDB reps database

% /usr/bin/time -v sourmash gather signatures/7a65a45d2739f50e5ecf920825482629.sig gtdb-r202.genomic-reps.k31.zip

== This is sourmash version 4.0.1.dev37+gad03e1ec. ==
== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==

select query k=31 automatically.
loaded query: GCF_002950215.1 Shigella flexn... (k=31, DNA)
loaded 1 databases.

Using EXPERIMENTAL feature: prefetch enabled!

overlap     p_query p_match avg_abund
---------   ------- ------- ---------
4.4 Mbp      100.0%  100.0%       1.1    GCF_002950215.1 Shigella flexneri 2a ...

found 1 matches total;
the recovered matches hit 100.0% of the query

        Command being timed: "sourmash gather signatures/7a65a45d2739f50e5ecf920825482629.sig gtdb-r202.genomic-reps.k31.zip"
        User time (seconds): 230.93
        System time (seconds): 2.18
        Percent of CPU this job got: 100%
        Elapsed (wall clock) time (h:mm:ss or m:ss): 3:51.08
        Average shared text size (kbytes): 0
        Average unshared data size (kbytes): 0
        Average stack size (kbytes): 0
        Average total size (kbytes): 0
        Maximum resident set size (kbytes): 85020
        Average resident set size (kbytes): 0
        Major (requiring I/O) page faults: 19
        Minor (reclaiming a frame) page faults: 33958
        Voluntary context switches: 1474
        Involuntary context switches: 301629
        Swaps: 0
        File system inputs: 2725792
        File system outputs: 8
        Socket messages sent: 0
        Socket messages received: 0
        Signals delivered: 0
        Page size (bytes): 4096
        Exit status: 0

against ~300k GTDB all database

% /usr/bin/time -v sourmash gather signatures/7a65a45d2739f50e5ecf920825482629.sig gtdb-r202.genomic.k31.zip

== This is sourmash version 4.0.1.dev37+gad03e1ec. ==
== Please cite Brown and Irber (2016), doi:10.21105/joss.00027. ==

select query k=31 automatically.
loaded query: GCF_002950215.1 Shigella flexn... (k=31, DNA)
loaded 1 databases.

Using EXPERIMENTAL feature: prefetch enabled!

overlap     p_query p_match avg_abund
---------   ------- ------- ---------
4.4 Mbp      100.0%  100.0%       1.1    GCF_002950215.1 Shigella flexneri 2a ...

found 1 matches total;
the recovered matches hit 100.0% of the query

        Command being timed: "sourmash gather signatures/7a65a45d2739f50e5ecf920825482629.sig gtdb-r202.genomic.k31.zip"
        User time (seconds): 1414.67
        System time (seconds): 14.18
        Percent of CPU this job got: 100%
        Elapsed (wall clock) time (h:mm:ss or m:ss): 23:46.91
        Average shared text size (kbytes): 0
        Average unshared data size (kbytes): 0
        Average stack size (kbytes): 0
        Average total size (kbytes): 0
        Maximum resident set size (kbytes): 4016316
        Average resident set size (kbytes): 0
        Major (requiring I/O) page faults: 20
        Minor (reclaiming a frame) page faults: 1611759
        Voluntary context switches: 2750
        Involuntary context switches: 546447
        Swaps: 0
        File system inputs: 16377032
        File system outputs: 8
        Socket messages sent: 0
        Socket messages received: 0
        Signals delivered: 0
        Page size (bytes): 4096
        Exit status: 0

[luizirber]

Nice!

A question about the usage of gather: We always pass a query argument, but the current use case doesn't modify the query (other than removing the best match). Because of this we can "optimize" the counter, which can be calculated only once and updated (because the updated query is always a subset of the original query).

Are we trying to cover use cases where the query can change in other ways, especially when it becomes something that is not a subset of the original query anymore? Technically this is supported with the current API (even tho no one does it), but more broadly it precludes the optimization I described (because the counter might be inconsistent and need to be recalculated).

[ctb]

Nice!

👍

A question about the usage of gather: We always pass a query argument, but the current use case doesn't modify the query (other than removing the best match). Because of this we can "optimize" the counter, which can be calculated only once and updated (because the updated query is always a subset of the original query).

Yes, and I should have an assert statement in there to make sure it's a subset of the original too..

In various refactorings (there were several that you didn't see 😆 ), I made use of the query in various ways, but I think we might be back to only making use of scaled, which can change (upwards), I think? It's not entirely clear to me even after staring at the code a whole bunch over the last few days. I would want to write a bunch of tests before committing to anything!

I like your idea of optimizing the query in that way.

Are we trying to cover use cases where the query can change in other ways, especially when it becomes something that is not a subset of the original query anymore? Technically this is supported with the current API (even tho no one does it), but more broadly it precludes the optimization I described (because the counter might be inconsistent and need to be recalculated).

Nope! I am not trying to cover use cases beyond what is already there, and definitely want to keep it as subset of original query (which I think makes perfect sense given min-set-cov understanding).

I actually want to tighten restrictions more, but this code is tricky and it took a while to get it all working right. If and when you think we should move forward with this approach, I can start writing clearer specifications and simpler tests on the Python side to tie the code down more. So, just say the word :).

[ctb]

ok @luizirber I think the tests are pretty thorough now, although I'd like to do a bit more white-box testing of consume if I can figure out a good way to do it.

[ctb]

done. planning to wait to ask for detailed review/merge until #1370 is in, however; off the top of my head, I don't see any point in making #1370 a bigger PR :).

[luizirber]

This looks great =]

I think it can be merged into #1370 now, even if it makes it longer to review again.