Repository: sourmash-bio/sourmash (language: python)

Issue:
MRG: Calculate all gather stats in rust; use for rocksdb gather
This PR adds a `calculate_gather_stats` function (and criterion benchmark) to calculate all gather statistics directly in rust. It makes use of this function in `disk_revindex.rs`, meaning that gather on a rocksdb database will now produce the full `GatherResult`. 

Note that this PR does not change the python layer to use these stats. This is just a step on the way to that :).

## New functions:

- `Sketch::MinHash::sum_abunds()` - sum abundances in a MinHash
- `Sketch::MinHash::n_unique_kmers()` - len(self) * self.scaled.
- `Sketch::inflate(abunds_from)` - inflate abundances in self using abund values from `abunds_from` MinHash
- `Sketch::inflated_abundances(abunds_from)` - same process as `inflate`, but just return the abundance vector and sum of abundances.
- `ani_utils::ani_from_containment` - a streamlined function to get ANI from containment. I suspect a lot of the time added for ANI calculation (ref: testing over in branchwater) is the way we calculate it in python (including calculating a total number of k-mers). This excludes anything not really needed for basic ANI without confidence intervals.
- `ani_utils::ani_from_containment_ci` - get high and low CI ANI from containment
- `index::calculate_gather_stats` - calculate all gather statistics that can be separated from the gather iteration and return the full `GatherResult`

The "expensive" calculations:
- abundance-weighted values require re-calculating intersections and/or manipulating the abundance vector
  > Thankfully I think the new `inflate`/`inflated_abundances` code is pretty efficient? It uses the `itertools` `merge_join_by` strategy from https://github.com/sourmash-bio/sourmash/compare/lirber/itertools_merge. Thanks @luizirber :)
- abundance `median` and `stddev` currently require cloning abunds, not sure how big of an issue that is.
- ANI with confidence intervals

## Punted Issues
- #3020: calculate all gather stats for all other database types. start with `linear.rs` since it uses `GatherResult` struct already.
- #3021 try moving `calculate_gather_stats` outside of the `gather` functions to enable postprocessing in parallel.
- #3022 : write ffi for the new functions?
  - `sketch::minhash::sum_abunds`
  - `sketch::minhash::n_unique_kmers`,
  - `sketch::minhash::inflate`
  - `sketch::minhash::inflated_abundances`
  - `ani_utils::ani_from_containment`
  - `ani_utils::ani_from_containment_ci`
  - `ani_utils::prob_nothing_in_common`

ref https://github.com/sourmash-bio/sourmash_plugin_branchwater/issues/187.



Discussion / hints:

Fix patch:
diff --git a/Cargo.lock b/Cargo.lock
index 1b11617c52..41a0843ea8 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -52,6 +52,15 @@ version = "1.0.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "41ed9a86bf92ae6580e0a31281f65a1b1d867c0cc68d5346e2ae128dddfa6a7d"
 
+[[package]]
+name = "approx"
+version = "0.5.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cab112f0a86d568ea0e627cc1d6be74a1e9cd55214684db5561995f6dad897c6"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "autocfg"
 version = "1.1.0"
@@ -867,6 +876,16 @@ dependencies = [
  "pkg-config",
 ]
 
+[[package]]
+name = "matrixmultiply"
+version = "0.3.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7574c1cf36da4798ab73da5b215bbf444f50718207754cb522201d78d1cd0ff2"
+dependencies = [
+ "autocfg",
+ "rawpointer",
+]
+
 [[package]]
 name = "md5"
 version = "0.7.0"
@@ -929,6 +948,35 @@ version = "0.0.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a2983372caf4480544083767bf2d27defafe32af49ab4df3a0b7fc90793a3664"
 
+[[package]]
+name = "nalgebra"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d506eb7e08d6329505faa8a3a00a5dcc6de9f76e0c77e4b75763ae3c770831ff"
+dependencies = [
+ "approx",
+ "matrixmultiply",
+ "nalgebra-macros",
+ "num-complex",
+ "num-rational",
+ "num-traits",
+ "rand",
+ "rand_distr",
+ "simba",
+ "typenum",
+]
+
+[[package]]
+name = "nalgebra-macros"
+version = "0.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "01fcc0b8149b4632adc89ac3b7b31a12fb6099a0317a4eb2ebff574ef7de7218"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn 1.0.104",
+]
+
 [[package]]
 name = "needletail"
 version = "0.5.1"
@@ -970,6 +1018,15 @@ dependencies = [
  "minimal-lexical",
 ]
 
+[[package]]
+name = "num-complex"
+version = "0.4.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "23c6602fda94a57c990fe0df199a035d83576b496aa29f4e634a8ac6004e68a6"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "num-integer"
 version = "0.1.44"
@@ -991,6 +1048,17 @@ dependencies = [
  "num-traits",
 ]
 
+[[package]]
+name = "num-rational"
+version = "0.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0638a1c9d0a3c0914158145bc76cff373a75a627e6ecbfb71cbe6f453a5a19b0"
+dependencies = [
+ "autocfg",
+ "num-integer",
+ "num-traits",
+]
+
 [[package]]
 name = "num-traits"
 version = "0.2.15"
@@ -1038,6 +1106,12 @@ dependencies = [
  "syn 2.0.46",
 ]
 
+[[package]]
+name = "paste"
+version = "1.0.14"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "de3145af08024dea9fa9914f381a17b8fc6034dfb00f3a84013f7ff43f29ed4c"
+
 [[package]]
 name = "peeking_take_while"
 version = "0.1.2"
@@ -1246,6 +1320,16 @@ dependencies = [
  "getrandom",
 ]
 
+[[package]]
+name = "rand_distr"
+version = "0.4.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "32cb0b9bc82b0a0876c2dd994a7e7a2683d3e7390ca40e6886785ef0c7e3ee31"
+dependencies = [
+ "num-traits",
+ "rand",
+]
+
 [[package]]
 name = "rand_xorshift"
 version = "0.3.0"
@@ -1255,6 +1339,12 @@ dependencies = [
  "rand_core",
 ]
 
+[[package]]
+name = "rawpointer"
+version = "0.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "60a357793950651c4ed0f3f52338f53b2f809f32d83a07f72909fa13e4c6c1e3"
+
 [[package]]
 name = "rayon"
 version = "1.8.1"
@@ -1354,6 +1444,12 @@ dependencies = [
  "librocksdb-sys",
 ]
 
+[[package]]
+name = "roots"
+version = "0.0.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "082f11ffa03bbef6c2c6ea6bea1acafaade2fd9050ae0234ab44a2153742b058"
+
 [[package]]
 name = "rustc-hash"
 version = "1.1.0"
@@ -1393,6 +1489,15 @@ version = "1.0.9"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "73b4b750c782965c211b42f022f59af1fbceabdd026623714f104152f1ec149f"
 
+[[package]]
+name = "safe_arch"
+version = "0.7.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f398075ce1e6a179b46f51bd88d0598b92b00d3551f1a2d4ac49e771b56ac354"
+dependencies = [
+ "bytemuck",
+]
+
 [[package]]
 name = "safemem"
 version = "0.3.3"
@@ -1469,6 +1574,19 @@ version = "1.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64"
 
+[[package]]
+name = "simba"
+version = "0.6.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f0b7840f121a46d63066ee7a99fc81dcabbc6105e437cae43528cea199b5a05f"
+dependencies = [
+ "approx",
+ "num-complex",
+ "num-traits",
+ "paste",
+ "wide",
+]
+
 [[package]]
 name = "simdutf8"
 version = "0.1.4"
@@ -1505,6 +1623,7 @@ dependencies = [
  "getrandom",
  "getset",
  "histogram",
+ "itertools 0.12.0",
  "log",
  "md5",
  "memmap2",
@@ -1523,8 +1642,11 @@ dependencies = [
  "rkyv",
  "roaring",
  "rocksdb",
+ "roots",
  "serde",
  "serde_json",
+ "statrs",
+ "streaming-stats",
  "tempfile",
  "thiserror",
  "twox-hash",
@@ -1541,6 +1663,28 @@ version = "1.1.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a2eb9349b6444b326872e140eb1cf5e7c522154d69e7a0ffb0fb81c06b37543f"
 
+[[package]]
+name = "statrs"
+version = "0.16.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2d08e5e1748192713cc281da8b16924fb46be7b0c2431854eadc785823e5696e"
+dependencies = [
+ "approx",
+ "lazy_static",
+ "nalgebra",
+ "num-traits",
+ "rand",
+]
+
+[[package]]
+name = "streaming-stats"
+version = "0.2.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "b0d670ce4e348a2081843569e0f79b21c99c91bb9028b3b3ecb0f050306de547"
+dependencies = [
+ "num-traits",
+]
+
 [[package]]
 name = "syn"
 version = "1.0.104"
@@ -1666,6 +1810,12 @@ dependencies = [
  "syn 2.0.46",
 ]
 
+[[package]]
+name = "typenum"
+version = "1.17.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "42ff0bf0c66b8238c6f3b578df37d0b7848e55df8577b3f74f92a69acceeb825"
+
 [[package]]
 name = "unarray"
 version = "0.1.4"
@@ -1837,6 +1987,16 @@ dependencies = [
  "wasm-bindgen",
 ]
 
+[[package]]
+name = "wide"
+version = "0.7.15"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "89beec544f246e679fc25490e3f8e08003bc4bf612068f325120dad4cea02c1c"
+dependencies = [
+ "bytemuck",
+ "safe_arch",
+]
+
 [[package]]
 name = "winapi"
 version = "0.3.9"
diff --git a/include/sourmash.h b/include/sourmash.h
index d647378da7..410a07eab0 100644
--- a/include/sourmash.h
+++ b/include/sourmash.h
@@ -37,6 +37,7 @@ enum SourmashErrorCode {
   SOURMASH_ERROR_CODE_READ_DATA = 1201,
   SOURMASH_ERROR_CODE_STORAGE = 1202,
   SOURMASH_ERROR_CODE_HLL_PRECISION_BOUNDS = 1301,
+  SOURMASH_ERROR_CODE_ANI_ESTIMATION_ERROR = 1401,
   SOURMASH_ERROR_CODE_IO = 100001,
   SOURMASH_ERROR_CODE_UTF8_ERROR = 100002,
   SOURMASH_ERROR_CODE_PARSE_INT = 100003,
diff --git a/src/core/Cargo.toml b/src/core/Cargo.toml
index 12d03ebc43..c177c0f807 100644
--- a/src/core/Cargo.toml
+++ b/src/core/Cargo.toml
@@ -38,6 +38,7 @@ fixedbitset = "0.4.0"
 getrandom = { version = "0.2", features = ["js"] }
 getset = "0.1.1"
 histogram = "0.9.1"
+itertools = "0.12.0"
 log = "0.4.20"
 md5 = "0.7.0"
 memmap2 = "0.9.4"
@@ -49,11 +50,14 @@ once_cell = "1.18.0"
 ouroboros = "0.18.3"
 piz = "0.5.0"
 primal-check = "0.3.1"
+rayon = { version = "1.8.1", optional = true }
 rkyv = { version = "0.7.44", optional = true }
 roaring = "0.10.3"
-rayon = { version = "1.8.1", optional = true }
+roots = "0.0.8"
 serde = { version = "1.0.196", features = ["derive"] }
 serde_json = "1.0.113"
+statrs = "0.16.0"
+streaming-stats = "0.2.3"
 thiserror = "1.0"
 twox-hash = "1.6.0"
 typed-builder = "0.18.0"
@@ -78,6 +82,10 @@ harness = false
 name = "minhash"
 harness = false
 
+[[bench]]
+name = "gather"
+harness = false
+
 [package.metadata.cargo-all-features]
 skip_optional_dependencies = true
 denylist = ["maturin"]
@@ -96,7 +104,7 @@ version = "0.3.68"
 features = ["console", "File"]
 
 [target.'cfg(all(target_arch = "wasm32"))'.dependencies.chrono]
-version = "0.4.34"
+version = "0.4.32"
 features = ["wasmbind"]
 
 [target.'cfg(all(target_arch = "wasm32", target_os="unknown"))'.dev-dependencies]
diff --git a/src/core/benches/gather.rs b/src/core/benches/gather.rs
new file mode 100644
index 0000000000..81fe2f87aa
--- /dev/null
+++ b/src/core/benches/gather.rs
@@ -0,0 +1,60 @@
+use std::fs::File;
+use std::io::BufReader;
+use std::path::PathBuf;
+
+use sourmash::collection::Collection;
+use sourmash::signature::Signature;
+use sourmash::sketch::Sketch;
+use sourmash::{index::calculate_gather_stats, storage::SigStore};
+
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+
+fn gather_stats_benchmarks(c: &mut Criterion) {
+    let mut filename = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
+    filename.push("../../tests/test-data/track_abund/47.fa.sig");
+    let file = File::open(filename).unwrap();
+    let reader = BufReader::new(file);
+    let sig = Signature::from_reader(reader)
+        .expect("Loading error")
+        .swap_remove(0);
+    let orig_query = sig.minhash().unwrap();
+    let query = orig_query.clone();
+    let total_weighted_hashes = orig_query.sum_abunds();
+
+    let mut filename = PathBuf::from(env!("CARGO_MANIFEST_DIR"));
+    filename.push("../../tests/test-data/track_abund/63.fa.sig");
+    // load collection to get sig in sigstore
+    let signatures = Signature::from_path(filename).expect("cant find file");
+    let collection = Collection::from_sigs(signatures).expect("cant make collection");
+    let match_sig: SigStore = collection.sig_for_dataset(0).expect("cant load sig");
+    let test_cases = vec![(false, false), (true, false), (false, true), (true, true)];
+
+    let mut group = c.benchmark_group("gather_stats");
+    for (calc_abund_stats, calc_ani_ci) in test_cases {
+        let test_name = format!(
+            "abund{}_ani_ci{}",
+            calc_abund_stats as u8, calc_ani_ci as u8
+        );
+        group.bench_function(&test_name, |b| {
+            b.iter(|| {
+                calculate_gather_stats(
+                    black_box(&orig_query),
+                    black_box(query.clone()),
+                    black_box(match_sig.clone()),
+                    black_box(42), // Example match_size
+                    black_box(1),  // Example gather_result_rank
+                    black_box(200),
+                    black_box(total_weighted_hashes.try_into().unwrap()),
+                    black_box(calc_abund_stats),
+                    black_box(calc_ani_ci),
+                    black_box(None), // don't set custom confidence intervals
+                )
+                .expect("error calculating gather stats");
+            });
+        });
+    }
+    group.finish();
+}
+
+criterion_group!(gather, gather_stats_benchmarks);
+criterion_main!(gather);
diff --git a/src/core/src/ani_utils.rs b/src/core/src/ani_utils.rs
new file mode 100644
index 0000000000..073cf6050b
--- /dev/null
+++ b/src/core/src/ani_utils.rs
@@ -0,0 +1,362 @@
+// Streamlined set of utils for containment --> ANI estimation
+// Equations based off of: https://github.com/KoslickiLab/mutation-rate-ci-calculator
+// Reference: https://doi.org/10.1101/2022.01.11.475870
+
+use roots::{find_root_brent, SimpleConvergency};
+use statrs::distribution::{ContinuousCDF, Normal};
+
+use crate::Error;
+
+fn exp_n_mutated(l: f64, k: f64, r1: f64) -> f64 {
+    let q = r1_to_q(k, r1);
+    l * q
+}
+
+fn var_n_mutated(l: f64, k: f64, r1: f64, q: Option<f64>) -> Result<f64, Error> {
+    if r1 == 0.0 {
+        return Ok(0.0);
+    }
+
+    let q = q.unwrap_or_else(|| r1_to_q(k, r1));
+
+    let var_n = l * (1.0 - q) * (q * (2.0 * k + (2.0 / r1) - 1.0) - 2.0 * k)
+        + k * (k - 1.0) * (1.0 - q).powi(2)
+        + (2.0 * (1.0 - q) / (r1.powi(2))) * ((1.0 + (k - 1.0) * (1.0 - q)) * r1 - q);
+
+    if var_n < 0.0 {
+        Err(Error::ANIEstimationError {
+            message: "varN is less than 0.0".into(),
+        })
+    } else {
+        Ok(var_n)
+    }
+}
+
+fn exp_n_mutated_squared(l: f64, k: f64, p: f64) -> Result<f64, Error> {
+    let var_n = var_n_mutated(l, k, p, None)?;
+    let exp_n_squared = exp_n_mutated(l, k, p).powi(2);
+    Ok(var_n + exp_n_squared)
+}
+
+fn probit(p: f64) -> f64 {
+    Normal::new(0.0, 1.0).unwrap().inverse_cdf(p)
+}
+
+fn r1_to_q(k: f64, r1: f64) -> f64 {
+    1.0 - (1.0 - r1).powi(k as i32)
+}
+
+// prior versions of ani estimation also returned this value by default.
+// BUT, it's not really something we need to calculate for prefetch/gather,
+// since there will always be something in common (or comparison would not be happening)
+// I'm not sure where the right place to put this back in is.. perhaps add minhash.contained_ani()
+// and return this there.
+// Usage:
+// let prob_nothing_in_common =
+//     get_exp_probability_nothing_common(n_unique_kmers, ksize, point_estimate, f_scaled)?;
+// let prob_threshold = prob_threshold.unwrap_or(1e-3);
+#[allow(dead_code)]
+fn get_exp_probability_nothing_common(
+    ani_estimate: f64,
+    ksize: f64,
+    f_scaled: f64,
+    n_unique_kmers: f64,
+) -> Result<f64, Error> {
+    if ani_estimate == 0.0 || ani_estimate == 1.0 {
+        Ok(1.0 - ani_estimate)
+    } else {
+        // Calculate the expected log probability.
+        let exp_nmut = exp_n_mutated(n_unique_kmers, ksize, 1.0 - ani_estimate);
+        let mut expected_log_probability = (n_unique_kmers - exp_nmut) * (1.0 - f_scaled).ln();
+
+        if expected_log_probability.is_infinite() {
+            expected_log_probability = f64::NEG_INFINITY;
+        }
+        // Return the exponential of the expected log probability.
+        Ok(expected_log_probability.exp())
+    }
+}
+
+/// Streamlined function for ANI from containment.
+/// todo: report ANI as % in 5.0?
+pub fn ani_from_containment(containment: f64, ksize: f64) -> f64 {
+    if containment == 0.0 {
+        0.0
+    } else if containment == 1.0 {
+        1.0
+    } else {
+        1.0 - (1.0 - containment.powf(1.0 / ksize))
+    }
+}
+
+// Calculate containment to ANI with confidence intervals
+pub fn ani_ci_from_containment(
+    containment: f64,
+    ksize: f64,
+    scaled: u64,
+    n_unique_kmers: u64,
+    confidence: Option<f64>,
+) -> Result<(f64, f64), Error> {
+    if containment == 0.0 {
+        return Ok((0.0, 0.0));
+    } else if containment == 1.0 {
+        return Ok((1.0, 1.0));
+    }
+    let confidence = confidence.unwrap_or(0.95);
+
+    // conversions needed throughout
+    let scaled_f64 = scaled as f64;
+    let f_scaled = 1.0 / scaled_f64;
+    let n_unique_kmers = n_unique_kmers as f64;
+    let alpha = 1.0 - confidence;
+
+    let z_alpha = probit(1.0 - alpha / 2.0);
+    let bias_factor = 1.0 - (1.0 - f_scaled).powi(n_unique_kmers as i32);
+    let term_1 = (1.0 - f_scaled) / (f_scaled * (n_unique_kmers).powi(3) * bias_factor.powi(2));
+    let term_2 = |pest: f64| {
+        n_unique_kmers * exp_n_mutated(n_unique_kmers, ksize, pest)
+            - exp_n_mutated_squared(n_unique_kmers, ksize, pest).unwrap_or(0.0)
+    };
+    let term_3 = |pest: f64| {
+        var_n_mutated(n_unique_kmers, ksize, pest, None).unwrap_or(0.0) / (n_unique_kmers).powi(2)
+    };
+
+    let var_direct = |pest: f64| term_1 * term_2(pest) + term_3(pest);
+
+    let f1 = |pest: f64| {
+        (1.0 - pest).powi(ksize as i32) + z_alpha * var_direct(pest).sqrt() - containment
+    };
+    let f2 = |pest: f64| {
+        (1.0 - pest).powi(ksize as i32) - z_alpha * var_direct(pest).sqrt() - containment
+    };
+
+    let mut convergency = SimpleConvergency {
+        eps: 1e-15,
+        max_iter: 1000,
+    };
+
+    let dist_sol1 =
+        find_root_brent(0.0000001, 0.9999999, &f1, &mut convergency).unwrap_or_default();
+    let dist_sol2 =
+        find_root_brent(0.0000001, 0.9999999, &f2, &mut convergency).unwrap_or_default();
+
+    Ok((1.0 - dist_sol1, 1.0 - dist_sol2))
+}
+
+#[cfg(test)]
+mod tests {
+
+    use super::*;
+    use std::f64::EPSILON;
+
+    #[test]
+    fn test_containment_to_ani_zero() {
+        let contain = 0.0;
+        let ksize = 21;
+        let scaled = 10;
+        let n_unique_kmers = 100;
+        let confidence = Some(0.95);
+        let res = ani_from_containment(contain, ksize as f64);
+        assert_eq!(res, 0.0);
+        let (ci_low, ci_high) =
+            ani_ci_from_containment(contain, ksize as f64, scaled, n_unique_kmers, confidence)
+                .unwrap();
+
+        eprintln!("{}", ci_low);
+        eprintln!("{}", ci_high);
+        assert_eq!(ci_low, 0.0);
+        assert_eq!(ci_high, 0.0);
+    }
+
+    #[test]
+    fn test_containment_to_ani_one() {
+        let contain = 1.0;
+        let ksize = 21;
+        let scaled = 10;
+        let n_unique_kmers = 100;
+        let confidence = None;
+        let res = ani_from_containment(contain, ksize as f64);
+        assert_eq!(res, 1.0);
+        let (ci_low, ci_high) =
+            ani_ci_from_containment(contain, ksize as f64, scaled, n_unique_kmers, confidence)
+                .unwrap();
+        assert_eq!(ci_low, 1.0);
+        assert_eq!(ci_high, 1.0);
+    }
+
+    #[test]
+    fn test_containment_to_ani_scaled1() {
+        let contain = 0.5;
+        let ksize = 21;
+        let scaled = 1;
+        let n_unique_kmers = 10000;
+        let confidence = None;
+        let ani = ani_from_containment(contain, ksize as f64);
+        assert!((ani - 0.9675317785238916) < EPSILON);
+        let (ci_low, ci_high) =
+            ani_ci_from_containment(contain, ksize as f64, scaled, n_unique_kmers, confidence)
+                .unwrap();
+        assert!((ci_low - 0.9635213980271021) < EPSILON);
+        assert!((ci_high - 0.9712900870335944) < EPSILON);
+    }
+
+    #[test]
+    fn test_containment_to_ani_scaled100() {
+        let contain = 0.1;
+        let ksize = 31;
+        let scaled = 100;
+        let n_unique_kmers = 10000;
+        let confidence = None;
+        let ani = ani_from_containment(contain, ksize as f64);
+        assert!((ani - 0.9284145445194744) < EPSILON);
+        let (ci_low, ci_high) =
+            ani_ci_from_containment(contain, ksize as f64, scaled, n_unique_kmers, confidence)
+                .unwrap();
+        assert!((ci_low - 0.9094445232754665) < EPSILON);
+        assert!((ci_high - 0.9467922076143345) < EPSILON);
+    }
+
+    #[test]
+    fn test_containment_to_ani_scaled100_2() {
+        let contain = 0.5;
+        let ksize = 21;
+        let scaled = 100;
+        let n_unique_kmers = 10000;
+        let confidence = None;
+        let ani = ani_from_containment(contain, ksize as f64);
+        assert!((ani - 0.9675317785238916) < EPSILON);
+        let (ci_low, ci_high) =
+            ani_ci_from_containment(contain, ksize as f64, scaled, n_unique_kmers, confidence)
+                .unwrap();
+        assert!((ci_low - 0.9569003945603415) < EPSILON);
+        assert!((ci_high - 0.9762879360833708) < EPSILON);
+    }
+
+    #[test]
+    fn test_var_n_mutated_zero() {
+        let r = 0.0;
+        let ksize = 31;
+        let nkmers = 200;
+        let var_n_mut = var_n_mutated(nkmers as f64, ksize as f64, r, None).unwrap(); // Assuming the function returns a Result
+        assert_eq!(var_n_mut, 0.0, "Expected variance to be 0 for r=0");
+    }
+
+    #[test]
+    fn test_var_n_mutated_value_error() {
+        let r = 10.0;
+        let ksize = 31;
+        let nkmers = 200;
+        match var_n_mutated(nkmers as f64, ksize as f64, r, None) {
+            Err(e) => assert_eq!(
+                e.to_string(),
+                "error while calculating ANI confidence intervals: varN is less than 0.0",
+                "Unexpected error message"
+            ),
+            Ok(_) => panic!("Expected an error, but got Ok"),
+        }
+    }
+
+    #[test]
+    fn test_var_n_mutated_success() {
+        let r = 0.4;
+        let ksize = 31;
+        let nkmers = 200_000;
+        let var_n_mut = var_n_mutated(nkmers as f64, ksize as f64, r, None).unwrap(); // Assuming the function returns a Result
+        let expected = 0.10611425440741508;
+        assert!(
+            (var_n_mut - expected).abs() < f64::EPSILON,
+            "Variance did not match expected value"
+        );
+    }
+
+    #[test]
+    fn test_r1_to_q() {
+        let k = 2.0;
+        let r1 = 0.5;
+        let result = r1_to_q(k, r1);
+        let expected = 0.75;
+
+        assert!(
+            (result - expected).abs() < EPSILON,
+            "The result of r1_to_q({}, {}) was {}, but {} was expected",
+            k,
+            r1,
+            result,
+            expected
+        );
+    }
+
+    #[test]
+    fn test_exp_n_mutated() {
+        let l = 100.0;
+        let k = 2.0;
+        let r1 = 0.5;
+
+        // Calculate the expected result based on the inputs
+        let expected_q = r1_to_q(k, r1);
+        let expected_result = l * expected_q;
+
+        let result = exp_n_mutated(l, k, r1);
+
+        assert!(
+            (result - expected_result).abs() < EPSILON,
+            "The result of exp_n_mutated({}, {}, {}) was {}, but {} was expected",
+            l,
+            k,
+            r1,
+            result,
+            expected_result
+        );
+    }
+
+    #[test]
+    fn test_get_exp_probability_nothing_common_ani_zero() {
+        let ani_estimate = 0.0;
+        let ksize = 31.0;
+        let f_scaled = 0.1;
+        let n_unique_kmers = 1000.0;
+
+        let result =
+            get_exp_probability_nothing_common(ani_estimate, ksize, f_scaled, n_unique_kmers)
+                .unwrap();
+        assert_eq!(
+            result, 1.0,
+            "Expected probability for ani_estimate of 0 to be 1.0"
+        );
+    }
+
+    #[test]
+    fn test_get_exp_probability_nothing_common_ani_one() {
+        let ani_estimate = 1.0;
+        let ksize = 31.0;
+        let f_scaled = 0.1;
+        let n_unique_kmers = 1000.0;
+
+        let result =
+            get_exp_probability_nothing_common(ani_estimate, ksize, f_scaled, n_unique_kmers)
+                .unwrap();
+        assert_eq!(
+            result, 0.0,
+            "Expected probability for ani_estimate of 1 to be 0.0"
+        );
+    }
+
+    #[test]
+    fn test_get_exp_probability_nothing_common() {
+        let contain = 0.1;
+        let ksize = 31 as f64;
+        let scaled = 10;
+        let f_scaled = 1.0 / scaled as f64;
+        let n_unique_kmers = 1000;
+
+        let ani = ani_from_containment(contain, ksize);
+        let result =
+            get_exp_probability_nothing_common(ani, ksize, f_scaled, n_unique_kmers as f64)
+                .unwrap();
+        assert!(
+            result >= 0.0 && result <= 1.0,
+            "The result should be a valid probability"
+        );
+        assert!((result - 0.000026561398887587855) < EPSILON);
+    }
+}
diff --git a/src/core/src/errors.rs b/src/core/src/errors.rs
index 28dea7a268..90c028eb38 100644
--- a/src/core/src/errors.rs
+++ b/src/core/src/errors.rs
@@ -46,6 +46,9 @@ pub enum SourmashError {
     #[error("Set error rate to a value smaller than 0.367696 and larger than 0.00203125")]
     HLLPrecisionBounds,
 
+    #[error("error while calculating ANI confidence intervals: {message}")]
+    ANIEstimationError { message: String },
+
     #[error(transparent)]
     ReadDataError(#[from] ReadDataError),
 
@@ -111,6 +114,8 @@ pub enum SourmashErrorCode {
     Storage = 12_02,
     // HLL errors
     HLLPrecisionBounds = 13_01,
+    // ANI errors
+    ANIEstimationError = 14_01,
     // external errors
     Io = 100_001,
     Utf8Error = 100_002,
@@ -144,6 +149,7 @@ impl SourmashErrorCode {
             SourmashError::ReadDataError { .. } => SourmashErrorCode::ReadData,
             SourmashError::StorageError { .. } => SourmashErrorCode::Storage,
             SourmashError::HLLPrecisionBounds { .. } => SourmashErrorCode::HLLPrecisionBounds,
+            SourmashError::ANIEstimationError { .. } => SourmashErrorCode::ANIEstimationError,
             SourmashError::SerdeError { .. } => SourmashErrorCode::SerdeError,
             SourmashError::IOError { .. } => SourmashErrorCode::Io,
             SourmashError::NifflerError { .. } => SourmashErrorCode::NifflerError,
diff --git a/src/core/src/ffi/minhash.rs b/src/core/src/ffi/minhash.rs
index 11863ba265..7d7c050fb0 100644
--- a/src/core/src/ffi/minhash.rs
+++ b/src/core/src/ffi/minhash.rs
@@ -199,7 +199,8 @@ pub unsafe extern "C" fn kmerminhash_remove_many(
     };
 
     // FIXME: proper exception here
-    mh.remove_many(hashes).expect("Hash removal error");
+    mh.remove_many(hashes.iter().copied())
+        .expect("Hash removal error");
 }
 
 ffi_fn! {
diff --git a/src/core/src/index/linear.rs b/src/core/src/index/linear.rs
index ff919b6f57..f58bbd1542 100644
--- a/src/core/src/index/linear.rs
+++ b/src/core/src/index/linear.rs
@@ -11,7 +11,7 @@ use crate::collection::CollectionSet;
 use crate::encodings::Idx;
 use crate::index::{GatherResult, Index, Selection, SigCounter};
 use crate::selection::Select;
-use crate::signature::{Signature, SigsTrait};
+use crate::signature::SigsTrait;
 use crate::sketch::minhash::KmerMinHash;
 use crate::sketch::Sketch;
 use crate::storage::SigStore;
@@ -150,13 +150,13 @@ impl LinearIndex {
             .internal_location()
             .into();
         let match_sig = self.collection.sig_for_dataset(dataset_id)?;
-        let result = self.stats_for_match(&match_sig, query, match_size, match_path, round)?;
+        let result = self.stats_for_match(match_sig, query, match_size, match_path, round)?;
         Ok(result)
     }
 
     fn stats_for_match(
         &self,
-        match_sig: &Signature,
+        match_sig: SigStore,
         query: &KmerMinHash,
         match_size: usize,
         match_path: PathBuf,
@@ -181,16 +181,27 @@ impl LinearIndex {
         let intersect_bp = (match_mh.scaled() * intersect_orig) as usize;
 
         let f_unique_to_query = intersect_orig as f64 / query.size() as f64;
-        let match_ = match_sig.clone();
+        let match_ = match_sig;
 
         // TODO: all of these
         let f_unique_weighted = 0.;
-        let average_abund = 0;
-        let median_abund = 0;
-        let std_abund = 0;
+        let average_abund = 0.;
+        let median_abund = 0.;
+        let std_abund = 0.;
         let md5 = "".into();
         let f_match_orig = 0.;
         let remaining_bp = 0;
+        let total_weighted_hashes = 0;
+        let n_unique_weighted_found = 0;
+        let query_containment_ani = 0.0;
+        let match_containment_ani = 0.0;
+        let max_containment_ani = 0.0;
+        let average_containment_ani = 0.0;
+        let query_containment_ani_ci_low = None;
+        let query_containment_ani_ci_high = None;
+        let match_containment_ani_ci_low = None;
+        let match_containment_ani_ci_high = None;
+        let sum_weighted_found = 0;
 
         Ok(GatherResult {
             intersect_bp,
@@ -209,6 +220,17 @@ impl LinearIndex {
             unique_intersect_bp,
             gather_result_rank,
             remaining_bp,
+            sum_weighted_found,
+            total_weighted_hashes,
+            n_unique_weighted_found,
+            query_containment_ani,
+            query_containment_ani_ci_low,
+            query_containment_ani_ci_high,
+            match_containment_ani,
+            match_containment_ani_ci_low,
+            match_containment_ani_ci_high,
+            max_containment_ani,
+            average_containment_ani,
         })
     }
 
diff --git a/src/core/src/index/mod.rs b/src/core/src/index/mod.rs
index ec55249b04..09891e4532 100644
--- a/src/core/src/index/mod.rs
+++ b/src/core/src/index/mod.rs
@@ -14,13 +14,18 @@ pub mod search;
 use std::path::Path;
 
 use getset::{CopyGetters, Getters, Setters};
-
 use serde::{Deserialize, Serialize};
+use stats::{median, stddev};
 use typed_builder::TypedBuilder;
 
+use crate::ani_utils::{ani_ci_from_containment, ani_from_containment};
 use crate::encodings::Idx;
 use crate::index::search::{search_minhashes, search_minhashes_containment};
 use crate::prelude::*;
+use crate::selection::Selection;
+use crate::signature::SigsTrait;
+use crate::sketch::minhash::KmerMinHash;
+use crate::storage::SigStore;
 use crate::Result;
 
 #[derive(TypedBuilder, CopyGetters, Getters, Setters, Serialize, Deserialize, Debug, PartialEq)]
@@ -34,11 +39,20 @@ pub struct GatherResult {
     #[getset(get_copy = "pub")]
     f_match: f64,
 
+    #[getset(get_copy = "pub")]
     f_unique_to_query: f64,
+
+    #[getset(get_copy = "pub")]
     f_unique_weighted: f64,
-    average_abund: usize,
-    median_abund: usize,
-    std_abund: usize,
+
+    #[getset(get_copy = "pub")]
+    average_abund: f64,
+
+    #[getset(get_copy = "pub")]
+    median_abund: f64,
+
+    #[getset(get_copy = "pub")]
+    std_abund: f64,
 
     #[getset(get = "pub")]
     filename: String,
@@ -50,17 +64,61 @@ pub struct GatherResult {
     md5: String,
 
     #[serde(skip)]
-    match_: Signature,
+    match_: SigStore,
 
+    #[getset(get_copy = "pub")]
     f_match_orig: f64,
+
+    #[getset(get_copy = "pub")]
     unique_intersect_bp: usize,
+
+    #[getset(get_copy = "pub")]
     gather_result_rank: usize,
+
+    #[getset(get_copy = "pub")]
     remaining_bp: usize,
+
+    #[getset(get_copy = "pub")]
+    n_unique_weighted_found: usize,
+
+    #[getset(get_copy = "pub")]
+    total_weighted_hashes: usize,
+
+    #[getset(get_copy = "pub")]
+    sum_weighted_found: usize,
+
+    #[getset(get_copy = "pub")]
+    query_containment_ani: f64,
+
+    #[getset(get_copy = "pub")]
+    #[serde(skip_serializing_if = "Option::is_none")]
+    query_containment_ani_ci_low: Option<f64>,
+
+    #[getset(get_copy = "pub")]
+    #[serde(skip_serializing_if = "Option::is_none")]
+    query_containment_ani_ci_high: Option<f64>,
+
+    #[getset(get_copy = "pub")]
+    match_containment_ani: f64,
+
+    #[getset(get_copy = "pub")]
+    #[serde(skip_serializing_if = "Option::is_none")]
+    match_containment_ani_ci_low: Option<f64>,
+
+    #[getset(get_copy = "pub")]
+    #[serde(skip_serializing_if = "Option::is_none")]
+    match_containment_ani_ci_high: Option<f64>,
+
+    #[getset(get_copy = "pub")]
+    average_containment_ani: f64,
+
+    #[getset(get_copy = "pub")]
+    max_containment_ani: f64,
 }
 
 impl GatherResult {
     pub fn get_match(&self) -> Signature {
-        self.match_.clone()
+        self.match_.clone().into()
     }
 }
 
@@ -145,3 +203,228 @@ where
         (*self).containment(other)
     }
 }
+
+// note all mh should be selected/downsampled prior to being passed in here.
+#[allow(clippy::too_many_arguments)]
+pub fn calculate_gather_stats(
+    orig_query: &KmerMinHash,
+    query: KmerMinHash,
+    match_sig: SigStore,
+    match_size: usize,
+    gather_result_rank: usize,
+    sum_weighted_found: usize,
+    total_weighted_hashes: usize,
+    calc_abund_stats: bool,
+    calc_ani_ci: bool,
+    confidence: Option<f64>,
+) -> Result<GatherResult> {
+    // get match_mh
+    let match_mh = match_sig.minhash().unwrap();
+    //bp remaining in subtracted query
+    let remaining_bp = (query.size() - match_size) * query.scaled() as usize;
+
+    // stats for this match vs original query
+    let (intersect_orig, _) = match_mh.intersection_size(orig_query).unwrap();
+    let intersect_bp = (match_mh.scaled() * intersect_orig) as usize;
+    let f_orig_query = intersect_orig as f64 / orig_query.size() as f64;
+    let f_match_orig = intersect_orig as f64 / match_mh.size() as f64;
+
+    // stats for this match vs current (subtracted) query
+    let f_match = match_size as f64 / match_mh.size() as f64;
+    let unique_intersect_bp = match_mh.scaled() as usize * match_size;
+    let f_unique_to_query = match_size as f64 / query.size() as f64;
+
+    // // get ANI values
+    let ksize = match_mh.ksize() as f64;
+    let query_containment_ani = ani_from_containment(f_unique_to_query, ksize);
+    let match_containment_ani = ani_from_containment(f_match, ksize);
+    let mut query_containment_ani_ci_low = None;
+    let mut query_containment_ani_ci_high = None;
+    let mut match_containment_ani_ci_low = None;
+    let mut match_containment_ani_ci_high = None;
+
+    if calc_ani_ci {
+        let n_unique_kmers = match_mh.n_unique_kmers();
+        let (qani_low, qani_high) = ani_ci_from_containment(
+            f_unique_to_query,
+            ksize,
+            match_mh.scaled(),
+            n_unique_kmers,
+            confidence,
+        )?;
+        query_containment_ani_ci_low = Some(qani_low);
+        query_containment_ani_ci_high = Some(qani_high);
+
+        let (mani_low, mani_high) = ani_ci_from_containment(
+            f_match,
+            ksize,
+            match_mh.scaled(),
+            n_unique_kmers,
+            confidence,
+        )?;
+        match_containment_ani_ci_low = Some(mani_low);
+        match_containment_ani_ci_high = Some(mani_high);
+    }
+
+    let average_containment_ani = (query_containment_ani + match_containment_ani) / 2.0;
+    let max_containment_ani = f64::max(query_containment_ani, match_containment_ani);
+
+    // set up non-abundance weighted values
+    let mut f_unique_weighted = f_unique_to_query;
+    let mut average_abund = 1.0;
+    let mut median_abund = 1.0;
+    let mut std_abund = 0.0;
+    // should these default to the unweighted numbers?
+    let mut n_unique_weighted_found = 0;
+    let mut sum_total_weighted_found = 0;
+
+    // If abundance, calculate abund-related metrics (vs current query)
+    if calc_abund_stats {
+        // need current downsampled query here to get f_unique_weighted
+        let (abunds, unique_weighted_found) = match match_mh.inflated_abundances(&query) {
+            Ok((abunds, unique_weighted_found)) => (abunds, unique_weighted_found),
+            Err(e) => {
+                return Err(e);
+            }
+        };
+        n_unique_weighted_found = unique_weighted_found as usize;
+        sum_total_weighted_found = sum_weighted_found + n_unique_weighted_found;
+        f_unique_weighted = n_unique_weighted_found as f64 / total_weighted_hashes as f64;
+
+        average_abund = n_unique_weighted_found as f64 / abunds.len() as f64;
+
+        // todo: try to avoid clone for these?
+        median_abund = median(abunds.iter().cloned()).unwrap();
+        std_abund = stddev(abunds.iter().cloned());
+    }
+
+    let result = GatherResult::builder()
+        .intersect_bp(intersect_bp)
+        .f_orig_query(f_orig_query)
+        .f_match(f_match)
+        .f_unique_to_query(f_unique_to_query)
+        .f_unique_weighted(f_unique_weighted)
+        .average_abund(average_abund)
+        .median_abund(median_abund)
+        .std_abund(std_abund)
+        .filename(match_sig.filename())
+        .name(match_sig.name())
+        .md5(match_sig.md5sum())
+        .match_(match_sig)
+        .f_match_orig(f_match_orig)
+        .unique_intersect_bp(unique_intersect_bp)
+        .gather_result_rank(gather_result_rank)
+        .remaining_bp(remaining_bp)
+        .n_unique_weighted_found(n_unique_weighted_found)
+        .query_containment_ani(query_containment_ani)
+        .query_containment_ani_ci_low(query_containment_ani_ci_low)
+        .query_containment_ani_ci_high(query_containment_ani_ci_high)
+        .match_containment_ani_ci_low(match_containment_ani_ci_low)
+        .match_containment_ani_ci_high(match_containment_ani_ci_high)
+        .match_containment_ani(match_containment_ani)
+        .average_containment_ani(average_containment_ani)
+        .max_containment_ani(max_containment_ani)
+        .sum_weighted_found(sum_total_weighted_found)
+        .total_weighted_hashes(total_weighted_hashes)
+        .build();
+    Ok(result)
+}
+
+#[cfg(test)]
+mod test_calculate_gather_stats {
+    use super::*;
+    use crate::cmd::ComputeParameters;
+    use crate::encodings::HashFunctions;
+    use crate::signature::Signature;
+    use crate::sketch::minhash::KmerMinHash;
+    use crate::sketch::Sketch;
+    // use std::f64::EPSILON;
+    // TODO: use f64::EPSILON when we bump MSRV
+    const EPSILON: f64 = 0.01;
+
+    #[test]
+    fn test_calculate_gather_stats() {
+        let scaled = 10;
+        let params = ComputeParameters::builder()
+            .ksizes(vec![31])
+            .scaled(scaled)
+            .build();
+
+        let mut match_sig = Signature::from_params(&params);
+        // create two minhash
+        let mut match_mh = KmerMinHash::new(scaled, 31, HashFunctions::Murmur64Dna, 42, true, 0);
+        match_mh.add_hash_with_abundance(1, 5);
+        match_mh.add_hash_with_abundance(3, 3);
+        match_mh.add_hash_with_abundance(5, 2);
+        match_mh.add_hash_with_abundance(8, 2);
+        match_mh.add_hash_with_abundance(11, 2); // Non-matching hash
+
+        match_sig.reset_sketches();
+        match_sig.push(Sketch::MinHash(match_mh.clone()));
+        match_sig.set_filename("match-filename");
+        match_sig.set_name("match-name");
+
+        eprintln!("num_sketches: {:?}", match_sig.size());
+        eprintln!("match_md5: {:?}", match_sig.md5sum());
+
+        // Setup orig_query minhash with abundances and non-matching hash
+        let mut orig_query = KmerMinHash::new(scaled, 31, HashFunctions::Murmur64Dna, 42, true, 0);
+        orig_query.add_hash_with_abundance(1, 3);
+        orig_query.add_hash_with_abundance(3, 2);
+        orig_query.add_hash_with_abundance(5, 1);
+        orig_query.add_hash_with_abundance(6, 1); // Non-matching hash
+        orig_query.add_hash_with_abundance(8, 1);
+        orig_query.add_hash_with_abundance(10, 1); // Non-matching hash
+
+        let query = orig_query.clone();
+        let total_weighted_hashes = orig_query.sum_abunds();
+
+        let match_size = 4;
+        let gather_result_rank = 0;
+        let calc_abund_stats = true;
+        let calc_ani_ci = false;
+        let result = calculate_gather_stats(
+            &orig_query,
+            query,
+            match_sig.into(),
+            match_size,
+            gather_result_rank,
+            0,
+            total_weighted_hashes.try_into().unwrap(),
+            calc_abund_stats,
+            calc_ani_ci,
+            None,
+        )
+        .unwrap();
+        // first, print all results
+        assert_eq!(result.filename(), "match-filename");
+        assert_eq!(result.name(), "match-name");
+        assert_eq!(result.md5(), "f54b271a62fb7e2856e7b8a33e741b6e");
+        assert_eq!(result.gather_result_rank, 0);
+        assert_eq!(result.remaining_bp, 20);
+
+        // results from match vs current query
+        assert_eq!(result.f_match, 0.8);
+        assert_eq!(result.unique_intersect_bp, 40);
+        assert_eq!(result.f_unique_to_query, 4.0 / 6.0);
+        eprintln!("{}", result.f_unique_weighted);
+        assert_eq!(result.f_unique_weighted, 7. / 9.);
+        assert_eq!(result.average_abund, 1.75);
+        assert_eq!(result.median_abund, 1.5);
+        assert_eq!(result.std_abund, 0.82915619758885);
+
+        // results from match vs orig_query
+        assert_eq!(result.intersect_bp, 40);
+        assert_eq!(result.f_orig_query, 4.0 / 6.0);
+        assert_eq!(result.f_match_orig, 4.0 / 5.0);
+
+        assert!((result.average_containment_ani - 0.98991665567826).abs() < EPSILON);
+        assert!((result.match_containment_ani - 0.9928276657672302).abs() < EPSILON);
+        assert!((result.query_containment_ani - 0.9870056455892898).abs() < EPSILON);
+        assert!((result.max_containment_ani - 0.9928276657672302).abs() < EPSILON);
+
+        assert_eq!(result.total_weighted_hashes, 9);
+        assert_eq!(result.n_unique_weighted_found, 7);
+        assert_eq!(result.sum_weighted_found, 7);
+    }
+}
diff --git a/src/core/src/index/revindex/disk_revindex.rs b/src/core/src/index/revindex/disk_revindex.rs
index f26c3c9593..afa74b8f80 100644
--- a/src/core/src/index/revindex/disk_revindex.rs
+++ b/src/core/src/index/revindex/disk_revindex.rs
@@ -11,14 +11,13 @@ use rocksdb::{ColumnFamilyDescriptor, MergeOperands, Options};
 use crate::collection::{Collection, CollectionSet};
 use crate::encodings::{Color, Idx};
 use crate::index::revindex::{
-    self as module, prepare_query, stats_for_cf, Datasets, DbStats, HashToColor, QueryColors,
-    RevIndexOps, DB, HASHES, MANIFEST, METADATA, STORAGE_SPEC, VERSION,
+    self as module, stats_for_cf, Datasets, DbStats, HashToColor, QueryColors, RevIndexOps, DB,
+    HASHES, MANIFEST, METADATA, STORAGE_SPEC, VERSION,
 };
-use crate::index::{GatherResult, SigCounter};
+use crate::index::{calculate_gather_stats, GatherResult, SigCounter};
 use crate::manifest::Manifest;
 use crate::prelude::*;
-use crate::signature::SigsTrait;
-use crate::sketch::minhash::KmerMinHash;
+use crate::sketch::minhash::{KmerMinHash, KmerMinHashBTree};
 use crate::sketch::Sketch;
 use crate::storage::{InnerStorage, Storage};
 use crate::Result;
@@ -167,7 +166,7 @@ impl RevIndex {
         // save DB version
         // TODO: probably should go together with a more general
         //       saving procedure used in create/update
-        self.db.put_cf(&cf_metadata, VERSION, &[DB_VERSION])?;
+        self.db.put_cf(&cf_metadata, VERSION, [DB_VERSION])?;
 
         // write manifest
         let mut wtr = vec![];
@@ -312,8 +311,18 @@ impl RevIndexOps for RevIndex {
     ) -> Result<Vec<GatherResult>> {
         let mut match_size = usize::max_value();
         let mut matches = vec![];
-        //let mut query: KmerMinHashBTree = orig_query.clone().into();
-        let selection = selection.unwrap_or_else(|| self.collection.selection());
+        let mut query = KmerMinHashBTree::from(orig_query.clone());
+        let mut sum_weighted_found = 0;
+        let _selection = selection.unwrap_or_else(|| self.collection.selection());
+        let mut orig_query_ds = orig_query.clone();
+        let total_weighted_hashes = orig_query.sum_abunds();
+
+        // or set this with user --track-abundance?
+        let calc_abund_stats = orig_query.track_abundance();
+
+        // todo: let user pass these options in
+        let calc_ani_ci = false;
+        let ani_confidence_interval_fraction = None;
 
         while match_size > threshold && !counter.is_empty() {
             trace!("counter len: {}", counter.len());
@@ -322,57 +331,39 @@ impl RevIndexOps for RevIndex {
             let (dataset_id, size) = counter.k_most_common_ordered(1)[0];
             match_size = if size >= threshold { size } else { break };
 
+            // this should downsample mh for us
             let match_sig = self.collection.sig_for_dataset(dataset_id)?;
 
-            // Calculate stats
-            let f_orig_query = match_size as f64 / orig_query.size() as f64;
-            let name = match_sig.name();
+            // get downsampled minhashes for comparison.
+            let match_mh = match_sig.minhash().unwrap().clone();
+            query = query.downsample_scaled(match_mh.scaled())?;
+            orig_query_ds = orig_query_ds.downsample_scaled(match_mh.scaled())?;
+
+            // just calculate essentials here
             let gather_result_rank = matches.len();
-            let match_ = match_sig.clone();
-            let md5 = match_sig.md5sum();
-
-            let match_mh = prepare_query(match_sig.into(), &selection)
-                .expect("Couldn't find a compatible MinHash");
-            let f_match = match_size as f64 / match_mh.size() as f64;
-            let unique_intersect_bp = match_mh.scaled() as usize * match_size;
-            let (intersect_orig, _) = match_mh.intersection_size(orig_query)?;
-            let intersect_bp = (match_mh.scaled() * intersect_orig) as usize;
-            let f_unique_to_query = intersect_orig as f64 / orig_query.size() as f64;
-
-            // TODO: all of these
-            let filename = "".into();
-            let f_unique_weighted = 0.;
-            let average_abund = 0;
-            let median_abund = 0;
-            let std_abund = 0;
-            let f_match_orig = 0.;
-            let remaining_bp = 0;
-
-            let result = GatherResult::builder()
-                .intersect_bp(intersect_bp)
-                .f_orig_query(f_orig_query)
-                .f_match(f_match)
-                .f_unique_to_query(f_unique_to_query)
-                .f_unique_weighted(f_unique_weighted)
-                .average_abund(average_abund)
-                .median_abund(median_abund)
-                .std_abund(std_abund)
-                .filename(filename)
-                .name(name)
-                .md5(md5)
-                .match_(match_.into())
-                .f_match_orig(f_match_orig)
-                .unique_intersect_bp(unique_intersect_bp)
-                .gather_result_rank(gather_result_rank)
-                .remaining_bp(remaining_bp)
-                .build();
-            matches.push(result);
+
+            // Calculate stats
+            let gather_result = calculate_gather_stats(
+                &orig_query_ds,
+                KmerMinHash::from(query.clone()),
+                match_sig,
+                match_size,
+                gather_result_rank,
+                sum_weighted_found,
+                total_weighted_hashes.try_into().unwrap(),
+                calc_abund_stats,
+                calc_ani_ci,
+                ani_confidence_interval_fraction,
+            )?;
+            // keep track of the sum weighted found
+            sum_weighted_found = gather_result.sum_weighted_found();
+            matches.push(gather_result);
 
             trace!("Preparing counter for next round");
             // Prepare counter for finding the next match by decrementing
             // all hashes found in the current match in other datasets
             // TODO: not used at the moment, so just skip.
-            //query.remove_many(match_mh.to_vec().as_slice())?;
+            query.remove_many(match_mh.iter_mins().copied())?; // is there a better way?
 
             // TODO: Use HashesToColors here instead. If not initialized,
             //       build it.
diff --git a/src/core/src/index/revindex/mem_revindex.rs b/src/core/src/index/revindex/mem_revindex.rs
index 5264c8550d..20f4c3aabd 100644
--- a/src/core/src/index/revindex/mem_revindex.rs
+++ b/src/core/src/index/revindex/mem_revindex.rs
@@ -92,7 +92,7 @@ impl RevIndex {
         // If threshold is zero, let's merge all queries and save time later
         let merged_query = queries.and_then(|qs| Self::merge_queries(qs, threshold));
 
-        let collection = Collection::from_paths(search_sigs)?.select(&selection)?;
+        let collection = Collection::from_paths(search_sigs)?.select(selection)?;
         let linear = LinearIndex::from_collection(collection.try_into()?);
 
         Ok(linear.index(threshold, merged_query, queries))
@@ -108,7 +108,7 @@ impl RevIndex {
         // If threshold is zero, let's merge all queries and save time later
         let merged_query = queries.and_then(|qs| Self::merge_queries(qs, threshold));
 
-        let collection = Collection::from_zipfile(zipfile)?.select(&selection)?;
+        let collection = Collection::from_zipfile(zipfile)?.select(selection)?;
         let linear = LinearIndex::from_collection(collection.try_into()?);
 
         Ok(linear.index(threshold, merged_query, queries))
diff --git a/src/core/src/index/revindex/mod.rs b/src/core/src/index/revindex/mod.rs
index bce95b93f6..a9b7091c8b 100644
--- a/src/core/src/index/revindex/mod.rs
+++ b/src/core/src/index/revindex/mod.rs
@@ -234,7 +234,7 @@ pub fn prepare_query(search_sig: Signature, selection: &Selection) -> Option<Kme
     })
 }
 
-#[derive(Debug, Default, PartialEq, Clone)]
+#[derive(Debug, Default, Clone)]
 pub enum Datasets {
     #[default]
     Empty,
diff --git a/src/core/src/lib.rs b/src/core/src/lib.rs
index da383372a0..8676f683dc 100644
--- a/src/core/src/lib.rs
+++ b/src/core/src/lib.rs
@@ -27,7 +27,9 @@ pub mod prelude;
 
 pub mod cmd;
 
+pub mod ani_utils;
 pub mod collection;
+pub mod encodings;
 pub mod index;
 pub mod manifest;
 pub mod selection;
@@ -35,8 +37,6 @@ pub mod signature;
 pub mod sketch;
 pub mod storage;
 
-pub mod encodings;
-
 #[cfg(feature = "from-finch")]
 pub mod from;
 
diff --git a/src/core/src/manifest.rs b/src/core/src/manifest.rs
index aefedeb05f..ca01e81acf 100644
--- a/src/core/src/manifest.rs
+++ b/src/core/src/manifest.rs
@@ -16,7 +16,7 @@ use crate::signature::{Signature, SigsTrait};
 use crate::sketch::Sketch;
 use crate::Result;
 
-#[derive(Debug, Serialize, Deserialize, Clone, CopyGetters, Getters, Setters, PartialEq)]
+#[derive(Debug, Serialize, Deserialize, Clone, CopyGetters, Getters, Setters, PartialEq, Eq)]
 pub struct Record {
     #[getset(get = "pub", set = "pub")]
     internal_location: PathBuf,
@@ -379,7 +379,7 @@ mod test {
         filename = PathBuf::from(filename);
 
         let mut pathfile = File::create(&filename).unwrap();
-        write!(pathfile, "Valid line\n").unwrap();
+        writeln!(pathfile, "Valid line").unwrap();
         pathfile.write_all(&[0xED, 0xA0, 0x80]).unwrap(); // invalid UTF-8
 
         // load into manifest
diff --git a/src/core/src/signature.rs b/src/core/src/signature.rs
index 381d45c643..2c795a0ed2 100644
--- a/src/core/src/signature.rs
+++ b/src/core/src/signature.rs
@@ -1033,9 +1033,8 @@ mod test {
             sig.add_sequence(&record.seq(), false).unwrap();
         }
 
-        let sketches = sig.sketches();
-        assert_eq!(sketches.len(), 12);
-        for sk in sketches {
+        assert_eq!(sig.size(), 12);
+        for sk in sig.iter() {
             assert_eq!(sk.size(), 500);
         }
     }
@@ -1144,7 +1143,7 @@ mod test {
         // iterate and check scaled
         for sig in &sigs {
             let modified_sig = sig.clone().select(&selection).unwrap();
-            for sketch in modified_sig.sketches() {
+            for sketch in modified_sig.iter() {
                 if let Sketch::MinHash(mh) = sketch {
                     dbg!("scaled: {:?}", mh.scaled());
                     assert_eq!(mh.scaled(), 2000);
diff --git a/src/core/src/sketch/minhash.rs b/src/core/src/sketch/minhash.rs
index 85690a271f..24cdc9539f 100644
--- a/src/core/src/sketch/minhash.rs
+++ b/src/core/src/sketch/minhash.rs
@@ -6,6 +6,7 @@ use std::iter::{Iterator, Peekable};
 use std::str;
 use std::sync::Mutex;
 
+use itertools::Itertools;
 use serde::de::Deserializer;
 use serde::ser::{SerializeStruct, Serializer};
 use serde::{Deserialize, Serialize};
@@ -421,9 +422,9 @@ impl KmerMinHash {
         Ok(())
     }
 
-    pub fn remove_many(&mut self, hashes: &[u64]) -> Result<(), Error> {
+    pub fn remove_many<T: IntoIterator<Item = u64>>(&mut self, hashes: T) -> Result<(), Error> {
         for min in hashes {
-            self.remove_hash(*min);
+            self.remove_hash(min);
         }
         Ok(())
     }
@@ -739,6 +740,14 @@ impl KmerMinHash {
         Ok(new_mh)
     }
 
+    pub fn sum_abunds(&self) -> u64 {
+        if let Some(abunds) = &self.abunds {
+            abunds.iter().sum()
+        } else {
+            self.size() as u64
+        }
+    }
+
     pub fn to_vec_abunds(&self) -> Vec<(u64, u64)> {
         if let Some(abunds) = &self.abunds {
             self.mins
@@ -765,11 +774,78 @@ impl KmerMinHash {
         hll
     }
 
+    // Approximate total number of kmers
+    // this could be improved by generating an HLL estimate while sketching instead
+    // (for scaled minhashes)
+    pub fn n_unique_kmers(&self) -> u64 {
+        self.size() as u64 * self.scaled() // + (self.ksize - 1) for bp estimation
+    }
+
     // create a downsampled copy of self
     pub fn downsample_scaled(&self, scaled: u64) -> Result<KmerMinHash, Error> {
         let max_hash = max_hash_for_scaled(scaled);
         self.downsample_max_hash(max_hash)
     }
+
+    pub fn inflate(&mut self, abunds_from: &KmerMinHash) -> Result<(), Error> {
+        self.check_compatible(abunds_from)?;
+
+        // check that abunds_from has abundances
+        if abunds_from.abunds.is_none() {
+            return Err(Error::NeedsAbundanceTracking);
+        }
+
+        let self_iter = self.mins.iter();
+        let abunds_iter = abunds_from.abunds.as_ref().unwrap().iter();
+        let abunds_from_iter = abunds_from.mins.iter().zip(abunds_iter);
+
+        let (mins, abunds): (Vec<u64>, Vec<u64>) = self_iter
+            .merge_join_by(abunds_from_iter, |&self_val, &(other_val, _)| {
+                self_val.cmp(other_val)
+            })
+            .filter_map(|either| match either {
+                itertools::EitherOrBoth::Both(self_val, (_other_val, &other_abund)) => {
+                    Some((self_val, other_abund))
+                }
+                _ => None,
+            })
+            .unzip();
+
+        self.mins = mins;
+        self.abunds = Some(abunds);
+
+        self.reset_md5sum();
+        Ok(())
+    }
+
+    pub fn inflated_abundances(&self, abunds_from: &KmerMinHash) -> Result<(Vec<u64>, u64), Error> {
+        self.check_compatible(abunds_from)?;
+        // check that abunds_from has abundances
+        if abunds_from.abunds.is_none() {
+            return Err(Error::NeedsAbundanceTracking);
+        }
+
+        let self_iter = self.mins.iter();
+        let abunds_iter = abunds_from.abunds.as_ref().unwrap().iter();
+        let abunds_from_iter = abunds_from.mins.iter().zip(abunds_iter);
+
+        let (abundances, total_abundance): (Vec<u64>, u64) = self_iter
+            .merge_join_by(abunds_from_iter, |&self_val, &(other_val, _)| {
+                self_val.cmp(other_val)
+            })
+            .filter_map(|either| match either {
+                itertools::EitherOrBoth::Both(_self_val, (_other_val, other_abund)) => {
+                    Some(*other_abund)
+                }
+                _ => None,
+            })
+            .fold((Vec::new(), 0u64), |(mut acc_vec, acc_sum), abund| {
+                acc_vec.push(abund);
+                (acc_vec, acc_sum + abund)
+            });
+
+        Ok((abundances, total_abundance))
+    }
 }
 
 impl SigsTrait for KmerMinHash {
@@ -1271,9 +1347,9 @@ impl KmerMinHashBTree {
         }
     }
 
-    pub fn remove_many(&mut self, hashes: &[u64]) -> Result<(), Error> {
+    pub fn remove_many<T: IntoIterator<Item = u64>>(&mut self, hashes: T) -> Result<(), Error> {
         for min in hashes {
-            self.remove_hash(*min);
+            self.remove_hash(min);
         }
         Ok(())
     }
@@ -1539,6 +1615,14 @@ impl KmerMinHashBTree {
                 .collect()
         }
     }
+
+    pub fn sum_abunds(&self) -> u64 {
+        if let Some(abunds) = &self.abunds {
+            abunds.values().sum()
+        } else {
+            self.size() as u64
+        }
+    }
 }
 
 impl SigsTrait for KmerMinHashBTree {
diff --git a/src/core/src/storage.rs b/src/core/src/storage.rs
index 17cbb7701c..26ca5ba6f7 100644
--- a/src/core/src/storage.rs
+++ b/src/core/src/storage.rs
@@ -77,6 +77,15 @@ pub struct SigStore {
     data: OnceCell<Signature>,
 }
 
+impl PartialEq for SigStore {
+    fn eq(&self, other: &Self) -> bool {
+        self.filename == other.filename
+            && self.name == other.name
+            && self.metadata == other.metadata
+            && self.data == other.data
+    }
+}
+
 #[derive(Serialize, Deserialize)]
 pub(crate) struct StorageInfo {
     pub backend: String,
@@ -343,7 +352,7 @@ impl Storage for ZipStorage {
     }
 
     fn spec(&self) -> String {
-        format!("zip://{}", self.path().unwrap_or("".into()))
+        format!("zip://{}", self.path().unwrap_or_else(|| "".into()))
     }
 }
 


Test patch:
diff --git a/src/core/tests/minhash.rs b/src/core/tests/minhash.rs
index 12477ed0d2..195091056f 100644
--- a/src/core/tests/minhash.rs
+++ b/src/core/tests/minhash.rs
@@ -204,10 +204,10 @@ fn oracle_mins(hashes in vec(u64::ANY, 1..10000)) {
     }
 
     c.add_from(&a).unwrap();
-    c.remove_many(&to_remove).unwrap();
+    c.remove_many(to_remove.iter().copied()).unwrap();
 
     d.add_from(&b).unwrap();
-    d.remove_many(&to_remove).unwrap();
+    d.remove_many(to_remove.iter().copied()).unwrap();
 
     assert_eq!(a.mins(), b.mins());
     assert_eq!(c.mins(), d.mins());
@@ -245,8 +245,8 @@ fn oracle_mins_scaled(hashes in vec(u64::ANY, 1..10000)) {
     c.add_many(&hashes).unwrap();
     d.add_many(&hashes).unwrap();
 
-    c.remove_many(&to_remove).unwrap();
-    d.remove_many(&to_remove).unwrap();
+    c.remove_many(to_remove.iter().copied()).unwrap();
+    d.remove_many(to_remove.iter().copied()).unwrap();
 
     a.remove_hash(hashes[0]);
     b.remove_hash(hashes[0]);
@@ -789,3 +789,103 @@ fn seq_to_hashes_2(seq in "QRMTHINK") {
 }
 
 }
+
+#[test]
+fn test_inflate() {
+    // Setup minhash_a with some mins but no abundances
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Hp, 42, true, 0);
+    a.add_hash(10);
+    a.add_hash(20);
+    a.add_hash(30);
+
+    // Setup minhash_b with mins, some of which match a, and with abundances
+    let mut b = KmerMinHash::new(5, 3, HashFunctions::Murmur64Hp, 42, true, 0);
+    b.add_hash_with_abundance(10, 2);
+    b.add_hash_with_abundance(20, 4);
+    b.add_hash_with_abundance(40, 6); // Non-matching hash
+
+    // Attempt to inflate minhash_a using minhash_b's abundances
+    assert!(a.inflate(&b).is_ok());
+
+    a.inflate(&b).unwrap();
+    eprintln!("{:?}", a.to_vec_abunds());
+    assert_eq!(a.to_vec_abunds(), vec![(10, 2), (20, 4)]);
+}
+
+#[test]
+fn test_inflated_abundances() {
+    // Setup minhash_a with some mins but no abundances
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Hp, 42, false, 0);
+    a.add_hash(10);
+    a.add_hash(20);
+    a.add_hash(30);
+
+    // Setup minhash_b with mins, some of which match a, and with abundances
+    let mut b = KmerMinHash::new(5, 3, HashFunctions::Murmur64Hp, 42, true, 0);
+    b.add_hash_with_abundance(10, 2);
+    b.add_hash_with_abundance(20, 4);
+    b.add_hash_with_abundance(40, 9); // Non-matching hash
+
+    // Attempt to inflate minhash_a using minhash_b's abundances
+    assert!(a.inflate(&b).is_ok());
+
+    let (abunds, total_abund) = a.inflated_abundances(&b).unwrap();
+    assert_eq!(abunds, vec![2, 4]);
+    assert_eq!(total_abund, 6);
+}
+
+#[test]
+fn test_inflate_noabund() {
+    // Setup minhash a with some mins but no abundances
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Dna, 42, false, 0);
+    a.add_hash(10);
+    a.add_hash(20);
+    a.add_hash(30);
+    let b = a.clone();
+    let result = a.inflate(&b);
+    assert!(matches!(
+        result,
+        Err(sourmash::Error::NeedsAbundanceTracking)
+    ));
+}
+
+#[test]
+fn test_inflated_abunds_noabund() {
+    // Setup minhash a with some mins but no abundances
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Dna, 42, false, 0);
+    a.add_hash(10);
+    a.add_hash(20);
+    a.add_hash(30);
+    let result = a.inflated_abundances(&a);
+    assert!(matches!(
+        result,
+        Err(sourmash::Error::NeedsAbundanceTracking)
+    ));
+}
+
+#[test]
+fn test_sum_abunds() {
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Dna, 42, true, 0);
+    a.add_hash_with_abundance(10, 2);
+    a.add_hash_with_abundance(20, 4);
+    a.add_hash_with_abundance(40, 9);
+    assert_eq!(a.sum_abunds(), 15);
+}
+
+#[test]
+fn test_sum_abunds_noabund() {
+    let mut a = KmerMinHash::new(5, 3, HashFunctions::Murmur64Dna, 42, false, 0);
+    a.add_hash(10);
+    a.add_hash(20);
+    a.add_hash(30);
+    assert_eq!(a.sum_abunds(), 3);
+}
+
+#[test]
+fn test_n_unique_kmers() {
+    let mut mh = KmerMinHash::new(10, 21, HashFunctions::Murmur64Dna, 42, true, 0);
+    mh.add_hash(10);
+    mh.add_hash(20);
+    mh.add_hash(30);
+    assert_eq!(mh.n_unique_kmers(), 30)
+}