diff -r aec67cdd6b14 -r efa60bc4f743 src/measures/discrete.rs
--- a/src/measures/discrete.rs	Tue Aug 01 10:32:12 2023 +0300
+++ b/src/measures/discrete.rs	Thu Aug 29 00:00:00 2024 -0500
@@ -121,7 +121,13 @@
     /// Prune all spikes with zero mass.
     #[inline]
     pub fn prune(&mut self) {
-        self.spikes.retain(|δ| δ.α != F::ZERO);
+        self.prune_by(|δ| δ.α != F::ZERO);
+    }
+
+    /// Prune spikes by the predicate `g`.
+    #[inline]
+    pub fn prune_by<G : FnMut(&DeltaMeasure<Domain, F>) -> bool>(&mut self, g : G) {
+        self.spikes.retain(g);
     }
 
     /// Add the spikes produced by `iter` to this measure.
@@ -138,6 +144,48 @@
     pub fn push(&mut self, δ : DeltaMeasure<Domain, F>) {
         self.spikes.push(δ);
     }
+
+    /// Iterate over triples of masses and locations of two discrete measures, which are assumed
+    /// to have equal locations of same spike indices.
+    pub fn both_matching<'a>(&'a self, other : &'a DiscreteMeasure<Domain, F>) ->
+      impl Iterator<Item=(F, F, &'a Domain)> {
+        let m = self.len().max(other.len());
+        self.iter_spikes().map(Some).chain(std::iter::repeat(None))
+            .zip(other.iter_spikes().map(Some).chain(std::iter::repeat(None)))
+            .take(m)
+            .map(|(oδ, orδ)| {
+                match (oδ, orδ) {
+                    (Some(δ), Some(rδ)) => (δ.α, rδ.α, &δ.x), // Assumed δ.x=rδ.x
+                    (Some(δ), None)     => (δ.α, F::ZERO,  &δ.x),
+                    (None, Some(rδ))    => (F::ZERO, rδ.α, &rδ.x),
+                    (None, None)        => panic!("This cannot happen!"),
+                }
+            })
+    }
+
+    /// Subtract `other` from `self`, assuming equal locations of same spike indices
+    pub fn sub_matching(&self, other : &DiscreteMeasure<Domain, F>) -> DiscreteMeasure<Domain, F>
+    where Domain : Clone {
+        self.both_matching(other)
+            .map(|(α, β, x)| (x.clone(), α - β))
+            .collect()
+    }
+
+    /// Add `other` to `self`, assuming equal locations of same spike indices
+    pub fn add_matching(&self, other : &DiscreteMeasure<Domain, F>) -> DiscreteMeasure<Domain, F>
+    where Domain : Clone {
+        self.both_matching(other)
+            .map(|(α, β, x)| (x.clone(), α + β))
+            .collect()
+    }
+
+    /// Calculate the Radon-norm distance of `self` to `other`,
+    /// assuming equal locations of same spike indices.
+    pub fn dist_matching(&self, other : &DiscreteMeasure<Domain, F>) -> F where F : Float {
+        self.both_matching(other)
+            .map(|(α, β, _)| (α-β).abs())
+            .sum()
+    }
 }
 
 impl<Domain, F : Num> IntoIterator for DiscreteMeasure<Domain, F> {