feat: add Mash distance metric and rolling entropy support

Implement the Mash distance metric across the CLI, index, and compact vector traits. This includes adding a `Mash` variant to the `DistanceMetric` enum and `MetricArg` CLI argument, implementing the conversion from Jaccard distances using the standard mutation-rate estimator formula, and updating documentation with supported metrics and algorithmic references. Additionally, add an `entropy` method to rolling statistics for computing order-specific entropy.
This commit is contained in:
Eric Coissac
2026-07-09 11:40:48 +02:00
parent dc3392865f
commit 2610a4af79
7 changed files with 61 additions and 8 deletions
+23
View File
@@ -1,5 +1,16 @@
use ndarray::{Array1, Array2};
/// Convert a Jaccard distance matrix (`1 - J`) into a Mash distance matrix, per
/// https://mash.readthedocs.io/en/latest/distances.html:
/// `D = -1/k * ln(2J / (1+J))`.
fn jaccard_to_mash(d_jaccard: &Array2<f64>, k: usize) -> Array2<f64> {
d_jaccard.mapv(|d| {
let j = 1.0 - d;
if j <= 0.0 { 1.0 }
else { -1.0 / k as f64 * (2.0 * j / (1.0 + j)).ln() }
})
}
// ── Column-level weight statistic — total count or presence count per column.
/// Additive across layers and partitions; used as denominator in normalised distances.
///
@@ -74,6 +85,12 @@ pub trait CountPartials: ColumnWeights {
m
}
/// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived
/// from the presence-threshold Jaccard distance.
fn threshold_mash_dist_matrix(&self, k: usize, threshold: u32) -> Array2<f64> {
jaccard_to_mash(&self.threshold_jaccard_dist_matrix(threshold), k)
}
fn relfreq_bray_dist_matrix(&self) -> Array2<f64> {
let global = self.col_weights();
let mut m = self.partial_relfreq_bray(&global).mapv(|v| 1.0 - v);
@@ -126,6 +143,12 @@ pub trait BitPartials: ColumnWeights {
m
}
/// Mash distance (https://mash.readthedocs.io/en/latest/distances.html), derived
/// from the Jaccard distance.
fn mash_dist_matrix(&self, k: usize) -> Array2<f64> {
jaccard_to_mash(&self.jaccard_dist_matrix(), k)
}
fn hamming_dist_matrix(&self) -> Array2<u64> {
self.partial_hamming()
}
+4
View File
@@ -14,6 +14,8 @@ pub enum DistanceMetric {
Jaccard,
/// Hamming distance (number of differing kmer positions) on presence/absence data.
Hamming,
/// Mash distance on presence/absence data (Jaccard-derived mutation-rate estimate).
Mash,
/// Bray-Curtis dissimilarity on raw counts.
BrayCurtis,
/// Bray-Curtis dissimilarity normalised by per-genome total counts.
@@ -84,6 +86,7 @@ impl KmerIndex {
DistanceMetric::Hellinger => CountPartials::hellinger_dist_matrix(&global),
DistanceMetric::HellingerEuclidean => CountPartials::hellinger_euclidean_dist_matrix(&global),
DistanceMetric::Jaccard => CountPartials::threshold_jaccard_dist_matrix(&global, presence_threshold),
DistanceMetric::Mash => CountPartials::threshold_mash_dist_matrix(&global, self.kmer_size(), presence_threshold),
DistanceMetric::Hamming => {
return Err(OKIError::InvalidInput(
"Hamming is only available for presence/absence indexes".into(),
@@ -108,6 +111,7 @@ impl KmerIndex {
let matrix = match metric {
DistanceMetric::Jaccard => BitPartials::jaccard_dist_matrix(&global),
DistanceMetric::Mash => BitPartials::mash_dist_matrix(&global, self.kmer_size()),
DistanceMetric::Hamming => {
BitPartials::hamming_dist_matrix(&global).mapv(|v| v as f64)
}
+2
View File
@@ -10,6 +10,7 @@ use tracing::info;
#[derive(clap::ValueEnum, Clone, Copy, Debug)]
pub enum MetricArg {
Jaccard,
Mash,
Hamming,
BrayCurtis,
#[value(name = "relfreq-bray-curtis")]
@@ -26,6 +27,7 @@ impl From<MetricArg> for DistanceMetric {
fn from(m: MetricArg) -> Self {
match m {
MetricArg::Jaccard => DistanceMetric::Jaccard,
MetricArg::Mash => DistanceMetric::Mash,
MetricArg::Hamming => DistanceMetric::Hamming,
MetricArg::BrayCurtis => DistanceMetric::BrayCurtis,
MetricArg::RelfreqBrayCurtis => DistanceMetric::RelfreqBrayCurtis,
-7
View File
@@ -196,13 +196,6 @@ impl RollingStat {
.map(|raw| Minimizer::from_raw_unchecked(raw << (64 - self.m * 2)))
}
pub fn entropy(&self, order: usize) -> Option<f64> {
if !self.ready() {
return None;
}
Some(self.entropy.entropy(order))
}
pub fn normalized_entropy(&self) -> Option<f64> {
if !self.ready() {
return None;