[MRG] add max_containment to MinHash class.

src/sourmash/cli/compare.py

@@ -27,6 +27,10 @@ def subparser(subparsers):
         '--containment', action='store_true',
         help='calculate containment instead of similarity'
     )
+    subparser.add_argument(
+        '--max-containment', action='store_true',
+        help='calculate max containment instead of similarity'
+    )
     subparser.add_argument(
         '--from-file',
         help='a text file containing a list of files to load signatures from'

src/sourmash/cli/search.py

@@ -34,12 +34,16 @@ def subparser(subparsers):
     )
     subparser.add_argument(
         '--containment', action='store_true',
-        help='evaluate containment rather than similarity'
+        help='score based on containment rather than similarity'
+    )
+    subparser.add_argument(
+        '--max-containment', action='store_true',
+        help='score based on max containment rather than similarity'
     )
     subparser.add_argument(
         '--ignore-abundance', action='store_true',
         help='do NOT use k-mer abundances if present; note: has no effect if '
-        '--containment is specified'
+        '--containment or --max-containment is specified'
     )
     subparser.add_argument(
         '--scaled', metavar='FLOAT', type=float, default=0,

src/sourmash/commands.py

@@ -7,7 +7,8 @@
 import sys
 
 import screed
-from .compare import compare_all_pairs, compare_serial_containment
+from .compare import (compare_all_pairs, compare_serial_containment,
+                      compare_serial_max_containment)
 from . import MinHash
 from .sbtmh import load_sbt_index, create_sbt_index
 from . import signature as sig
@@ -91,16 +92,24 @@ def compare(args):
         error('cannot mix scaled signatures with bounded signatures')
         sys.exit(-1)
 
+    is_containment = False
+    if args.containment or args.max_containment:
+        is_containment = True
+
+        if args.containment and args.max_containment:
+            notify("ERROR: cannot specify both --containment and --max-containment!")
+            sys.exit(-1)
+
     # complain if --containment and not is_scaled
-    if args.containment and not is_scaled:
-        error('must use scaled signatures with --containment option')
+    if is_containment and not is_scaled:
+        error('must use scaled signatures with --containment and --max-containment')
         sys.exit(-1)
 
     # notify about implicit --ignore-abundance:
-    if args.containment:
+    if is_containment:
         track_abundances = any(( s.minhash.track_abundance for s in siglist ))
         if track_abundances:
-            notify('NOTE: --containment means signature abundances are flattened.')
+            notify('NOTE: --containment and --max-containment means signature abundances are ignored.')
 
     # if using --scaled, downsample appropriately
     printed_scaled_msg = False
@@ -127,6 +136,8 @@ def compare(args):
     labeltext = [str(item) for item in siglist]
     if args.containment:
         similarity = compare_serial_containment(siglist)
+    elif args.max_containment:
+        similarity = compare_serial_max_containment(siglist)
     else:
         similarity = compare_all_pairs(siglist, args.ignore_abundance,
                                        n_jobs=args.processes)
@@ -437,8 +448,14 @@ def search(args):
         query.minhash = query.minhash.downsample(scaled=args.scaled)
 
     # set up the search databases
+    is_containment = args.containment or args.max_containment
+    if is_containment:
+        if args.containment and args.max_containment:
+            notify("ERROR: cannot specify both --containment and --max-containment!")
+            sys.exit(-1)
+
     databases = sourmash_args.load_dbs_and_sigs(args.databases, query,
-                                                not args.containment)
+                                                not is_containment)
 
     # forcibly ignore abundances if query has no abundances
     if not query.minhash.track_abundance:
@@ -450,8 +467,11 @@ def search(args):
 
     # do the actual search
     results = search_databases(query, databases,
-                               args.threshold, args.containment,
-                               args.best_only, args.ignore_abundance,
+                               threshold=args.threshold,
+                               do_containment=args.containment,
+                               do_max_containment=args.max_containment,
+                               best_only=args.best_only,
+                               ignore_abundance=args.ignore_abundance,
                                unload_data=True)
 
     n_matches = len(results)
@@ -477,7 +497,8 @@ def search(args):
         notify("** reporting only one match because --best-only was set")
 
     if args.output:
-        fieldnames = ['similarity', 'name', 'filename', 'md5']
+        fieldnames = ['similarity', 'name', 'filename', 'md5',
+                      'query_filename', 'query_name', 'query_md5']
 
         with FileOutput(args.output, 'wt') as fp:
             w = csv.DictWriter(fp, fieldnames=fieldnames)
@@ -486,6 +507,7 @@ def search(args):
             for sr in results:
                 d = dict(sr._asdict())
                 del d['match']
+                del d['query']
                 w.writerow(d)
 
     # save matching signatures upon request

src/sourmash/compare.py

@@ -42,7 +42,7 @@ def compare_serial(siglist, ignore_abundance, downsample=False):
 def compare_serial_containment(siglist, downsample=False):
     """Compare all combinations of signatures and return a matrix
     of containments. Processes combinations serially on a single
-    process. Best to use when there is few signatures.
+    process. Best to only use when there are few signatures.
 
     :param list siglist: list of signatures to compare
     :param boolean downsample by scaled if True
@@ -61,6 +61,28 @@ def compare_serial_containment(siglist, downsample=False):
     return containments
 
 
+def compare_serial_max_containment(siglist, downsample=False):
+    """Compare all combinations of signatures and return a matrix
+    of max_containments. Processes combinations serially on a single
+    process. Best to only use when there are few signatures.
+
+    :param list siglist: list of signatures to compare
+    :param boolean downsample by scaled if True
+    :return: np.array similarity matrix
+    """
+    import numpy as np
+
+    n = len(siglist)
+
+    containments = np.ones((n, n))
+    for i in range(n):
+        for j in range(n):
+            containments[i][j] = siglist[j].max_containment(siglist[i],
+                                                         downsample=downsample)
+
+    return containments
+
+
 def similarity_args_unpack(args, ignore_abundance, downsample):
     """Helper function to unpack the arguments. Written to use in pool.imap
     as it can only be given one argument."""

src/sourmash/index.py

@@ -38,7 +38,9 @@ def find(self, search_fn, *args, **kwargs):
                 matches.append(node)
         return matches
 
-    def search(self, query, *args, **kwargs):
+    def search(self, query, threshold=None,
+               do_containment=False, do_max_containment=False,
+               ignore_abundance=False, **kwargs):
         """Return set of matches with similarity above 'threshold'.
 
         Results will be sorted by similarity, highest to lowest.
@@ -55,16 +57,18 @@ def search(self, query, *args, **kwargs):
         """
 
         # check arguments
-        if 'threshold' not in kwargs:
+        if threshold is None:
             raise TypeError("'search' requires 'threshold'")
-        threshold = kwargs['threshold']
+        threshold = float(threshold)
 
-        do_containment = kwargs.get('do_containment', False)
-        ignore_abundance = kwargs.get('ignore_abundance', False)
+        if do_containment and do_max_containment:
+            raise TypeError("'do_containment' and 'do_max_containment' cannot both be True")
 
         # configure search - containment? ignore abundance?
         if do_containment:
             query_match = lambda x: query.contained_by(x, downsample=True)
+        elif do_max_containment:
+            query_match = lambda x: query.max_containment(x, downsample=True)
         else:
             query_match = lambda x: query.similarity(
                 x, downsample=True, ignore_abundance=ignore_abundance)
@@ -73,9 +77,9 @@ def search(self, query, *args, **kwargs):
         matches = []
 
         for ss in self.signatures():
-            similarity = query_match(ss)
-            if similarity >= threshold:
-                matches.append((similarity, ss, self.filename))
+            score = query_match(ss)
+            if score >= threshold:
+                matches.append((score, ss, self.filename))
 
         # sort!
         matches.sort(key=lambda x: -x[0])

src/sourmash/lca/lca_db.py

@@ -300,7 +300,8 @@ def save(self, db_name):
 
             json.dump(save_d, fp)
 
-    def search(self, query, *args, **kwargs):
+    def search(self, query, threshold=None, do_containment=False,
+               do_max_containment=False, ignore_abundance=False, **kwargs):
         """Return set of matches with similarity above 'threshold'.
 
         Results will be sorted by similarity, highest to lowest.
@@ -319,18 +320,18 @@ def search(self, query, *args, **kwargs):
             return []
 
         # check arguments
-        if 'threshold' not in kwargs:
+        if threshold is None:
             raise TypeError("'search' requires 'threshold'")
-        threshold = kwargs['threshold']
-        do_containment = kwargs.get('do_containment', False)
-        ignore_abundance = kwargs.get('ignore_abundance', False)
+        threshold = float(threshold)
+
         mh = query.minhash
         if ignore_abundance:
             mh.track_abundance = False
 
         # find all the matches, then sort & return.
         results = []
-        for x in self._find_signatures(mh, threshold, do_containment):
+        for x in self._find_signatures(mh, threshold, do_containment,
+                                       do_max_containment):
             (score, match, filename) = x
             results.append((score, match, filename))
 
@@ -455,7 +456,8 @@ def _signatures(self):
         return sigd
 
     def _find_signatures(self, minhash, threshold, containment=False,
-                       ignore_scaled=False):
+                         max_containment=False,
+                         ignore_scaled=False):
         """
         Do a Jaccard similarity or containment search, yield results.
 
@@ -467,7 +469,7 @@ def _find_signatures(self, minhash, threshold, containment=False,
         if self.scaled > minhash.scaled:
             minhash = minhash.downsample(scaled=self.scaled)
         elif self.scaled < minhash.scaled and not ignore_scaled:
-            # note that containment can be calculated w/o matching scaled.
+            # note that containment cannot be calculated w/o matching scaled.
             raise ValueError("lca db scaled is {} vs query {}; must downsample".format(self.scaled, minhash.scaled))
 
         query_mins = set(minhash.hashes)
@@ -494,6 +496,9 @@ def _find_signatures(self, minhash, threshold, containment=False,
             # calculate the containment or similarity
             if containment:
                 score = count / len(query_mins)
+            elif max_containment:
+                denom = min((len(query_mins), match_size))
+                score = count / denom
             else:
                 # query_mins is size of query signature
                 # match_size is size of match signature

src/sourmash/minhash.py

@@ -543,6 +543,18 @@ def contained_by(self, other, downsample=False):
 
         return self.count_common(other, downsample) / len(self)
 
+    def max_containment(self, other, downsample=False):
+        """\
+        Calculate maximum containment.
+        """
+        if not (self.scaled and other.scaled):
+            raise TypeError("can only calculate containment for scaled MinHashes")
+        min_denom = min((len(self), len(other)))
+        if not min_denom:
+            return 0.0
+
+        return self.count_common(other, downsample) / min_denom
+
     def __add__(self, other):
         if not isinstance(other, MinHash):
             raise TypeError("can only add MinHash objects to MinHash objects!")

src/sourmash/sbt.py

@@ -337,7 +337,10 @@ def find(self, search_fn, *args, **kwargs):
 
         return matches
 
-    def search(self, query, *args, **kwargs):
+    def search(self, query, threshold=None,
+               ignore_abundance=False, do_containment=False,
+               do_max_containment=False, best_only=False,
+               unload_data=False, **kwargs):
         """Return set of matches with similarity above 'threshold'.
 
         Results will be sorted by similarity, highest to lowest.
@@ -350,15 +353,17 @@ def search(self, query, *args, **kwargs):
           * ignore_abundance: default False. If True, and query signature
             and database support k-mer abundances, ignore those abundances.
         """
-        from .sbtmh import search_minhashes, search_minhashes_containment
+        from .sbtmh import (search_minhashes, search_minhashes_containment,
+                            search_minhashes_max_containment)
         from .sbtmh import SearchMinHashesFindBest
         from .signature import SourmashSignature
 
-        threshold = kwargs['threshold']
-        ignore_abundance = kwargs.get('ignore_abundance', False)
-        do_containment = kwargs.get('do_containment', False)
-        best_only = kwargs.get('best_only', False)
-        unload_data = kwargs.get('unload_data', False)
+        if threshold is None:
+            raise TypeError("'search' requires 'threshold'")
+        threshold = float(threshold)
+
+        if do_containment and do_max_containment:
+            raise TypeError("'do_containment' and 'do_max_containment' cannot both be True")
 
         # figure out scaled value of tree, downsample query if needed.
         leaf = next(iter(self.leaves()))
@@ -378,13 +383,21 @@ def search(self, query, *args, **kwargs):
         if do_containment:
             search_fn = search_minhashes_containment
             query_match = lambda x: tree_query.contained_by(x, downsample=True)
+        elif do_max_containment:
+            search_fn = search_minhashes_max_containment
+            query_match = lambda x: tree_query.max_containment(x,
+                                                               downsample=True)
 
         if best_only:            # this needs to be reset for each SBT
             search_fn = SearchMinHashesFindBest().search
 
         # now, search!
         results = []
-        for leaf in self.find(search_fn, tree_query, threshold, unload_data=unload_data):
+
+        # here, self.find is used only to find candidate nodes;
+        for leaf in self.find(search_fn, tree_query, threshold,
+                              unload_data=unload_data):
+            # the actual calculation of node match is done here:
             similarity = query_match(leaf.data)
 
             # tree search should always/only return matches above threshold