MRG: Fix scaled determination in fastgather and fastmultigather

src/fastgather.rs

@@ -2,6 +2,7 @@
 use anyhow::Result;
 use sourmash::prelude::Select;
 use sourmash::selection::Selection;
+use sourmash::sketch::minhash::KmerMinHash;
 
 use crate::utils::{
     consume_query_by_gather, load_collection, load_sketches_above_threshold, write_prefetch,
@@ -12,8 +13,7 @@ use crate::utils::{
 pub fn fastgather(
     query_filepath: String,
     against_filepath: String,
-    threshold_bp: usize,
-    scaled: usize,
+    threshold_bp: u64,
     selection: Selection,
     gather_output: Option<String>,
     prefetch_output: Option<String>,
@@ -40,31 +40,35 @@ pub fn fastgather(
     let query_md5 = query_sig.md5sum();
 
     // clone here is necessary b/c we use full query_sig in consume_query_by_gather
-    let query_sig_ds = query_sig.select(&selection)?; // downsample
-    let query_mh = match query_sig_ds.try_into() {
+    let query_sig_ds = query_sig.select(&selection)?; // downsample as needed.
+    let query_mh: KmerMinHash = match query_sig_ds.try_into() {
         Ok(query_mh) => query_mh,
         Err(_) => {
             bail!("No query sketch matching selection parameters.");
         }
     };
+
+    let mut against_selection = selection;
+    let scaled = query_mh.scaled();
+    against_selection.set_scaled(scaled);
+
     // load collection to match against.
     let against_collection = load_collection(
         &against_filepath,
-        &selection,
+        &against_selection,
         ReportType::Against,
         allow_failed_sigpaths,
     )?;
 
     // calculate the minimum number of hashes based on desired threshold
-    let threshold_hashes: u64 = {
-        let x = threshold_bp / scaled;
+    let threshold_hashes = {
+        let x = threshold_bp / scaled as u64;
         if x > 0 {
             x
         } else {
             1
         }
-    }
-    .try_into()?;
+    };
 
     eprintln!(
         "using threshold overlap: {} {}",

src/fastmultigather.rs

@@ -30,7 +30,7 @@ pub fn fastmultigather(
     query_filepath: String,
     against_filepath: String,
     threshold_bp: u32,
-    scaled: Option<usize>,
+    scaled: Option<u32>,
     selection: Selection,
     allow_failed_sigpaths: bool,
     save_matches: bool,
@@ -47,7 +47,7 @@ pub fn fastmultigather(
     )?;
 
     let scaled = match scaled {
-        Some(s) => s as u32,
+        Some(s) => s,
         None => {
             let scaled = *query_collection.max_scaled().expect("no records!?");
             eprintln!(
@@ -59,10 +59,10 @@ pub fn fastmultigather(
     };
 
     let mut against_selection = selection;
-    against_selection.set_scaled(scaled as u32);
+    against_selection.set_scaled(scaled);
 
     let threshold_hashes: u64 = {
-        let x = threshold_bp / scaled;
+        let x = threshold_bp as u64 / scaled as u64;
         if x > 0 {
             x as u64
         } else {

src/fastmultigather_rocksdb.rs

@@ -34,6 +34,22 @@ pub fn fastmultigather_rocksdb(
         allow_failed_sigpaths,
     )?;
 
+    // set scaled from query collection.
+    let scaled = match selection.scaled() {
+        Some(s) => s,
+        None => {
+            let scaled = *query_collection.max_scaled().expect("no records!?");
+            eprintln!(
+                "Setting scaled={} based on max scaled in query collection",
+                scaled
+            );
+            scaled
+        }
+    };
+
+    let mut against_selection = selection;
+    against_selection.set_scaled(scaled);
+
     // set up a multi-producer, single-consumer channel.
     let (send, recv) =
         std::sync::mpsc::sync_channel::<BranchwaterGatherResult>(rayon::current_num_threads());
@@ -55,8 +71,8 @@ pub fn fastmultigather_rocksdb(
     let send = query_collection
         .par_iter()
         .filter_map(|(coll, _idx, record)| {
-            let threshold = threshold_bp / selection.scaled()?;
-            let ksize = selection.ksize()?;
+            let threshold = threshold_bp / against_selection.scaled().expect("scaled is not set!?");
+            let ksize = against_selection.ksize().expect("ksize not set!?");
 
             // query downsampling happens here
             match coll.sig_from_record(record) {
@@ -78,7 +94,7 @@ pub fn fastmultigather_rocksdb(
                             hash_to_color,
                             threshold as usize,
                             &query_mh,
-                            Some(selection.clone()),
+                            Some(against_selection.clone()),
                         );
                         if let Ok(matches) = matches {
                             for match_ in &matches {
@@ -101,8 +117,8 @@ pub fn fastmultigather_rocksdb(
                                     query_filename: query_filename.clone(),
                                     query_name: query_name.clone(),
                                     query_md5: query_md5.clone(),
-                                    query_bp: query_mh.n_unique_kmers() as usize,
-                                    ksize: ksize as usize,
+                                    query_bp: query_mh.n_unique_kmers(),
+                                    ksize: ksize as u16,
                                     moltype: query_mh.hash_function().to_string(),
                                     scaled: query_mh.scaled(),
                                     query_n_hashes: query_mh.size() as u64,
@@ -170,6 +186,11 @@ pub fn fastmultigather_rocksdb(
     let i: usize = processed_sigs.fetch_max(0, atomic::Ordering::SeqCst);
     eprintln!("DONE. Processed {} search sigs", i);
 
+    if i == 0 {
+        eprintln!("ERROR: no search sigs found!?");
+        do_fail = true;
+    }
+
     let skipped_paths = skipped_paths.load(atomic::Ordering::SeqCst);
     let failed_paths = failed_paths.load(atomic::Ordering::SeqCst);
     let failed_gathers = failed_gathers.load(atomic::Ordering::SeqCst);

src/lib.rs

@@ -34,13 +34,13 @@ fn do_manysearch(
     siglist_path: String,
     threshold: f64,
     ksize: u8,
-    scaled: usize,
+    scaled: Option<u32>,
     moltype: String,
     output_path: Option<String>,
     ignore_abundance: Option<bool>,
 ) -> anyhow::Result<u8> {
     let againstfile_path: PathBuf = siglist_path.clone().into();
-    let selection = build_selection(ksize, Some(scaled), &moltype);
+    let selection = build_selection(ksize, scaled, &moltype);
     eprintln!("selection scaled: {:?}", selection.scaled());
     let allow_failed_sigpaths = true;
 
@@ -88,21 +88,20 @@ fn do_manysearch(
 fn do_fastgather(
     query_filename: String,
     siglist_path: String,
-    threshold_bp: usize,
+    threshold_bp: u64,
     ksize: u8,
-    scaled: usize,
+    scaled: Option<u32>,
     moltype: String,
     output_path_prefetch: Option<String>,
     output_path_gather: Option<String>,
 ) -> anyhow::Result<u8> {
-    let selection = build_selection(ksize, Some(scaled), &moltype);
+    let selection = build_selection(ksize, scaled, &moltype);
     let allow_failed_sigpaths = true;
 
     match fastgather::fastgather(
         query_filename,
         siglist_path,
         threshold_bp,
-        scaled,
         selection,
         output_path_prefetch,
         output_path_gather,
@@ -122,9 +121,9 @@ fn do_fastgather(
 fn do_fastmultigather(
     query_filenames: String,
     siglist_path: String,
-    threshold_bp: usize,
+    threshold_bp: u64,
     ksize: u8,
-    scaled: Option<usize>,
+    scaled: Option<u32>,
     moltype: String,
     output_path: Option<String>,
     save_matches: bool,
@@ -191,16 +190,17 @@ fn set_global_thread_pool(num_threads: usize) -> PyResult<usize> {
 }
 
 #[pyfunction]
+#[pyo3(signature = (siglist, ksize, scaled, moltype, output, colors, use_internal_storage))]
 fn do_index(
     siglist: String,
     ksize: u8,
-    scaled: usize,
+    scaled: Option<u32>,
     moltype: String,
     output: String,
     colors: bool,
     use_internal_storage: bool,
 ) -> anyhow::Result<u8> {
-    let selection = build_selection(ksize, Some(scaled), &moltype);
+    let selection = build_selection(ksize, scaled, &moltype);
     let allow_failed_sigpaths = false;
     match index::index(
         siglist,
@@ -238,7 +238,7 @@ fn do_multisearch(
     siglist_path: String,
     threshold: f64,
     ksize: u8,
-    scaled: Option<usize>,
+    scaled: Option<u32>,
     moltype: String,
     estimate_ani: bool,
     output_path: Option<String>,
@@ -272,13 +272,13 @@ fn do_pairwise(
     siglist_path: String,
     threshold: f64,
     ksize: u8,
-    scaled: usize,
+    scaled: Option<u32>,
     moltype: String,
     estimate_ani: bool,
     write_all: bool,
     output_path: Option<String>,
 ) -> anyhow::Result<u8> {
-    let selection = build_selection(ksize, Some(scaled), &moltype);
+    let selection = build_selection(ksize, scaled, &moltype);
     let allow_failed_sigpaths = true;
     match pairwise::pairwise(
         siglist_path,

src/manysearch.rs

@@ -129,7 +129,7 @@ pub fn manysearch(
                                     query_md5: query.md5sum.clone(),
                                     match_name: against_name.clone(),
                                     containment: containment_query_in_target,
-                                    intersect_hashes: overlap as usize,
+                                    intersect_hashes: overlap as u64,
                                     match_md5: Some(against_md5.clone()),
                                     jaccard: Some(jaccard),
                                     max_containment: Some(max_containment),

src/manysearch_rocksdb.rs

@@ -89,7 +89,7 @@ pub fn manysearch_rocksdb(
                                     query_md5: query_sig.md5sum(),
                                     match_name: path.clone(),
                                     containment,
-                                    intersect_hashes: overlap,
+                                    intersect_hashes: overlap as u64,
                                     match_md5: None,
                                     jaccard: None,
                                     max_containment: None,

src/python/sourmash_plugin_branchwater/__init__.py

@@ -70,7 +70,7 @@ def __init__(self, p):
         p.add_argument(
             "-s",
             "--scaled",
-            default=1000,
+            default=None,
             type=int,
             help="scaled factor at which to do comparisons",
         )
@@ -172,9 +172,9 @@ def __init__(self, p):
         p.add_argument(
             "-s",
             "--scaled",
-            default=1000,
+            default=None,
             type=int,
-            help="scaled factor at which to do comparisons (default: 1000)",
+            help="scaled factor at which to do comparisons (default: determined from query)",
         )
         p.add_argument(
             "-m",
@@ -333,9 +333,9 @@ def __init__(self, p):
         p.add_argument(
             "-s",
             "--scaled",
-            default=1000,
+            default=None,
             type=int,
-            help="scaled factor at which to do comparisons",
+            help="scaled factor at which to select sketches (default: chosen as max from collection)",
         )
         p.add_argument(
             "-m",
@@ -512,7 +512,7 @@ def __init__(self, p):
         p.add_argument(
             "-s",
             "--scaled",
-            default=1000,
+            default=None,
             type=int,
             help="scaled factor at which to do comparisons",
         )

src/python/tests/test_fastgather.py

@@ -1325,3 +1325,40 @@ def test_fullres_vs_sourmash_gather(runtmp):
     fg_total_weighted_hashes = set(gather_df["total_weighted_hashes"])
     g_total_weighted_hashes = set(sourmash_gather_df["total_weighted_hashes"])
     assert fg_total_weighted_hashes == g_total_weighted_hashes == set([73489])
+
+
+def test_equal_matches(runtmp):
+    # check that equal matches get returned from fastgather
+    base = sourmash.MinHash(scaled=1, ksize=31, n=0)
+
+    a = base.copy_and_clear()
+    b = base.copy_and_clear()
+    c = base.copy_and_clear()
+
+    a.add_many(range(0, 1000))
+    b.add_many(range(1000, 2000))
+    c.add_many(range(0, 2000))
+
+    ss = sourmash.SourmashSignature(a, name="g_a")
+    sourmash.save_signatures([ss], open(runtmp.output("a.sig"), "wb"))
+    ss = sourmash.SourmashSignature(b, name="g_b")
+    sourmash.save_signatures([ss], open(runtmp.output("b.sig"), "wb"))
+    ss = sourmash.SourmashSignature(c, name="g_mg")
+    sourmash.save_signatures([ss], open(runtmp.output("mg.sig"), "wb"))
+
+    runtmp.sourmash("sig", "cat", "a.sig", "b.sig", "-o", "combined.sig.zip")
+
+    runtmp.sourmash(
+        "scripts",
+        "fastgather",
+        "mg.sig",
+        "combined.sig.zip",
+        "-o",
+        "out.csv",
+        "--threshold-bp",
+        "0",
+    )
+
+    df = pandas.read_csv(runtmp.output("out.csv"))
+    assert len(df) == 2
+    assert set(df["intersect_bp"]) == {1000}