Mercurial > repos > measures / changeset

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/Cargo.toml	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,38 @@
+[package]
+name = "measures"
+version = "0.1.0"
+edition = "2021"
+rust-version = "1.85"
+authors = ["Tuomo Valkonen <tuomov@iki.fi>"]
+description = "Representation of discrete measures"
+homepage = "https://tuomov.iki.fi/software/pointsource_algs/"
+repository = "https://tuomov.iki.fi/repos/pointsource_algs/"
+license-file = "LICENSE"
+keywords = [
+    "measure",
+]
+categories = ["mathematics", "science"]
+
+[dependencies.alg_tools]
+version = "~0.4.0-dev"
+path = "../alg_tools"
+default-features = false
+#features = ["nightly"]
+
+[dependencies]
+serde = { version = "1.0", features = ["derive"] }
+nalgebra = { version = "~0.33.0", features = ["rand-no-std"] }
+numeric_literals = "~0.2.0"
+# All of these dependencies are only required with the `pyo3` feature
+pyo3 = { version = "~0.27.0", optional = true }
+numpy = { version = "~0.27.0", features = ["nalgebra"], optional = true }
+
+[features]
+default = []
+pyo3 = ["dep:pyo3", "dep:numpy", "alg_tools/pyo3"]
+
+[build-dependencies]
+regex = "~1.11.0"
+
+[profile.release]
+debug = true

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README.md	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,4 @@
+# measures
+
+This crates implements measures, in particular, sums of Dirac measures.
+It was separated from `pointsource_algs`.

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/base.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,20 @@
+//! Basic definitions for measures
+
+use alg_tools::norms::{Norm, NormExponent};
+use alg_tools::types::Num;
+use serde::Serialize;
+
+/// This is used with [`Norm::norm`] to indicate that a Radon norm is to be computed.
+#[derive(Copy, Clone, Serialize, Debug)]
+pub struct Radon;
+impl NormExponent for Radon {}
+
+/// A trait for (Radon) measures.
+///
+/// Currently has no methods, just the requirement that the Radon norm be implemented.
+pub trait Measure<F: Num>: Norm<Radon, F> {
+    type Domain;
+}
+
+/// Decomposition of measures
+pub struct MeasureDecomp;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/delta.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,320 @@
+/*!
+This module implementes delta measures, i.e., single spikes $\alpha \delta_x$ for some
+location $x$ and mass $\alpha$.
+*/
+
+use super::base::*;
+use alg_tools::instance::{ClosedSpace, Instance, Space};
+use alg_tools::linops::{Linear, Mapping};
+use alg_tools::loc::Loc;
+use alg_tools::norms::Norm;
+use alg_tools::self_ownable;
+use alg_tools::types::*;
+use serde::ser::{Serialize, SerializeStruct, Serializer};
+use std::ops::{Div, DivAssign, Mul, MulAssign, Neg};
+
+/// Representation of a delta measure.
+///
+/// This is a single spike $\alpha \delta\_x$ for some location $x$ in `Domain` and
+/// a mass $\alpha$ in `F`.
+#[derive(Clone, Copy, Debug)]
+pub struct DeltaMeasure<Domain, F: Num> {
+    // This causes [`csv`] to crash.
+    //#[serde(flatten)]
+    /// Location of the spike
+    pub x: Domain,
+    /// Mass of the spike
+    pub α: F,
+}
+
+const COORDINATE_NAMES: &'static [&'static str] = &["x0", "x1", "x2", "x3", "x4", "x5", "x6", "x7"];
+
+// Need to manually implement serialisation as [`csv`] writer fails on
+// structs with nested arrays as well as with #[serde(flatten)].
+impl<F: Num, const N: usize> Serialize for DeltaMeasure<Loc<N, F>, F>
+where
+    F: Serialize,
+{
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        assert!(N <= COORDINATE_NAMES.len());
+
+        let mut s = serializer.serialize_struct("DeltaMeasure", N + 1)?;
+        for (i, e) in (0..).zip(self.x.iter()) {
+            s.serialize_field(COORDINATE_NAMES[i], e)?;
+        }
+        s.serialize_field("weight", &self.α)?;
+        s.end()
+    }
+}
+
+impl<Domain, F: Float> Measure<F> for DeltaMeasure<Domain, F> {
+    type Domain = Domain;
+}
+
+impl<Domain, F: Float> Norm<Radon, F> for DeltaMeasure<Domain, F> {
+    #[inline]
+    fn norm(&self, _: Radon) -> F {
+        self.α.abs()
+    }
+}
+
+// impl<Domain : PartialEq, F : Float> Dist<Radon, F> for DeltaMeasure<Domain, F> {
+//     #[inline]
+//     fn dist(&self, other : &Self, _ : Radon) -> F {
+//         if self.x == other. x {
+//             (self.α - other.α).abs()
+//         } else {
+//             self.α.abs() + other.α.abs()
+//         }
+//     }
+// }
+
+impl<Domain, G, F: Num> Mapping<G> for DeltaMeasure<Domain, F>
+where
+    Domain: Space,
+    G::Codomain: Mul<F, Output = G::Codomain>,
+    G: Mapping<Domain> + Clone + ClosedSpace,
+    for<'b> &'b Domain: Instance<Domain>,
+{
+    type Codomain = G::Codomain;
+
+    #[inline]
+    fn apply<I: Instance<G>>(&self, g: I) -> Self::Codomain {
+        g.eval(|g̃| g̃.apply(&self.x) * self.α)
+    }
+}
+
+impl<Domain, G, F: Num> Linear<G> for DeltaMeasure<Domain, F>
+where
+    Domain: Space,
+    G::Codomain: Mul<F, Output = G::Codomain>,
+    G: Mapping<Domain> + Clone + ClosedSpace,
+    for<'b> &'b Domain: Instance<Domain>,
+{
+}
+
+// /// Partial blanket implementation of [`DeltaMeasure`] as a linear functional of [`Mapping`]s.
+// /// A full blanket implementation is not possible due to annoying Rust limitations: only [`Apply`]
+// /// on a reference is implemented, but a consuming [`Apply`] has to be implemented on a case-by-case
+// /// basis, not because an implementation could not be written, but because the Rust trait system
+// /// chokes up.
+// impl<Domain, G, F : Num, V> Linear<G> for DeltaMeasure<Domain, F>
+// where G: for<'a> Apply<&'a Domain, Output = V>,
+//       V : Mul<F>,
+//       Self: Apply<G, Output =  <V as Mul<F>>::Output> {
+//     type Codomain = <V as Mul<F>>::Output;
+// }
+
+// impl<'b, Domain, G, F : Num, V> Apply<&'b G> for DeltaMeasure<Domain, F>
+// where G: for<'a> Apply<&'a Domain, Output = V>,
+//       V : Mul<F> {
+//     type Output = <V as Mul<F>>::Output;
+
+//     #[inline]
+//     fn apply(&self, g : &'b G) -> Self::Output {
+//         g.apply(&self.x) * self.α
+//     }
+// }
+
+// /// Implementation of the necessary apply for BTFNs
+// mod btfn_apply {
+//     use super::*;
+//     use alg_tools::bisection_tree::{BTFN, BTImpl, SupportGenerator, LocalAnalysis};
+
+//     impl<F : Float, BT, G, V, const N : usize> Apply<BTFN<F, G, BT, N>>
+//     for DeltaMeasure<Loc<N, F>, F>
+//     where BT : BTImpl< N, F>,
+//         G : SupportGenerator< N, F, Id=BT::Data>,
+//         G::SupportType : LocalAnalysis<F, BT::Agg, N> + for<'a> Apply<&'a Loc<N, F>, Output = V>,
+//         V : std::iter::Sum + Mul<F> {
+
+//         type Output = <V as Mul<F>>::Output;
+
+//         #[inline]
+//         fn apply(&self, g : BTFN<F, G, BT, N>) -> Self::Output {
+//             g.apply(&self.x) * self.α
+//         }
+//     }
+// }
+
+impl<D, Domain, F: Num> From<(D, F)> for DeltaMeasure<Domain, F>
+where
+    D: Into<Domain>,
+{
+    #[inline]
+    fn from((x, α): (D, F)) -> Self {
+        DeltaMeasure { x: x.into(), α: α }
+    }
+}
+
+impl<'a, Domain: Clone, F: Num> From<&'a DeltaMeasure<Domain, F>> for DeltaMeasure<Domain, F> {
+    #[inline]
+    fn from(d: &'a DeltaMeasure<Domain, F>) -> Self {
+        d.clone()
+    }
+}
+
+impl<Domain, F: Num> DeltaMeasure<Domain, F> {
+    /// Set the mass of the spike.
+    #[inline]
+    pub fn set_mass(&mut self, α: F) {
+        self.α = α
+    }
+
+    /// Set the location of the spike.
+    #[inline]
+    pub fn set_location(&mut self, x: Domain) {
+        self.x = x
+    }
+
+    /// Get the mass of the spike.
+    #[inline]
+    pub fn get_mass(&self) -> F {
+        self.α
+    }
+
+    /// Get a mutable reference to the mass of the spike.
+    #[inline]
+    pub fn get_mass_mut(&mut self) -> &mut F {
+        &mut self.α
+    }
+
+    /// Get a reference to the location of the spike.
+    #[inline]
+    pub fn get_location(&self) -> &Domain {
+        &self.x
+    }
+
+    /// Get a mutable reference to the location of the spike.
+    #[inline]
+    pub fn get_location_mut(&mut self) -> &mut Domain {
+        &mut self.x
+    }
+}
+
+impl<Domain, F: Num> IntoIterator for DeltaMeasure<Domain, F> {
+    type Item = Self;
+    type IntoIter = std::iter::Once<Self>;
+
+    #[inline]
+    fn into_iter(self) -> Self::IntoIter {
+        std::iter::once(self)
+    }
+}
+
+impl<'a, Domain, F: Num> IntoIterator for &'a DeltaMeasure<Domain, F> {
+    type Item = Self;
+    type IntoIter = std::iter::Once<Self>;
+
+    #[inline]
+    fn into_iter(self) -> Self::IntoIter {
+        std::iter::once(self)
+    }
+}
+
+macro_rules! make_delta_scalarop_rhs {
+    ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
+        impl<F: Num, Domain> $trait<F> for DeltaMeasure<Domain, F> {
+            type Output = Self;
+            fn $fn(mut self, b: F) -> Self {
+                self.α.$fn_assign(b);
+                self
+            }
+        }
+
+        impl<'a, F: Num, Domain> $trait<&'a F> for DeltaMeasure<Domain, F> {
+            type Output = Self;
+            fn $fn(mut self, b: &'a F) -> Self {
+                self.α.$fn_assign(*b);
+                self
+            }
+        }
+
+        impl<'b, F: Num, Domain: Clone> $trait<F> for &'b DeltaMeasure<Domain, F> {
+            type Output = DeltaMeasure<Domain, F>;
+            fn $fn(self, b: F) -> Self::Output {
+                DeltaMeasure { α: self.α.$fn(b), x: self.x.clone() }
+            }
+        }
+
+        impl<'a, 'b, F: Num, Domain: Clone> $trait<&'a F> for &'b DeltaMeasure<Domain, F> {
+            type Output = DeltaMeasure<Domain, F>;
+            fn $fn(self, b: &'a F) -> Self::Output {
+                DeltaMeasure { α: self.α.$fn(*b), x: self.x.clone() }
+            }
+        }
+
+        impl<F: Num, Domain> $trait_assign<F> for DeltaMeasure<Domain, F> {
+            fn $fn_assign(&mut self, b: F) {
+                self.α.$fn_assign(b)
+            }
+        }
+
+        impl<'a, F: Num, Domain> $trait_assign<&'a F> for DeltaMeasure<Domain, F> {
+            fn $fn_assign(&mut self, b: &'a F) {
+                self.α.$fn_assign(*b)
+            }
+        }
+    };
+}
+
+make_delta_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
+make_delta_scalarop_rhs!(Div, div, DivAssign, div_assign);
+
+macro_rules! make_delta_scalarop_lhs {
+    ($trait:ident, $fn:ident; $($f:ident)+) => { $(
+        impl<Domain> $trait<DeltaMeasure<Domain, $f>> for $f {
+            type Output = DeltaMeasure<Domain, $f>;
+            fn $fn(self, mut δ : DeltaMeasure<Domain, $f>) -> Self::Output {
+                δ.α = self.$fn(δ.α);
+                δ
+            }
+        }
+
+        impl<'a, Domain : Clone> $trait<&'a DeltaMeasure<Domain, $f>> for $f {
+            type Output = DeltaMeasure<Domain, $f>;
+            fn $fn(self, δ : &'a DeltaMeasure<Domain, $f>) -> Self::Output {
+                DeltaMeasure{ x : δ.x.clone(), α : self.$fn(δ.α) }
+            }
+        }
+
+        impl<'b, Domain> $trait<DeltaMeasure<Domain, $f>> for &'b $f {
+            type Output = DeltaMeasure<Domain, $f>;
+            fn $fn(self, mut δ : DeltaMeasure<Domain, $f>) -> Self::Output {
+                δ.α = self.$fn(δ.α);
+                δ
+            }
+        }
+
+        impl<'a, 'b, Domain : Clone> $trait<&'a DeltaMeasure<Domain, $f>> for &'b $f {
+            type Output = DeltaMeasure<Domain, $f>;
+            fn $fn(self, δ : &'a DeltaMeasure<Domain, $f>) -> Self::Output {
+                DeltaMeasure{ x : δ.x.clone(), α : self.$fn(δ.α) }
+            }
+        }
+    )+ }
+}
+
+make_delta_scalarop_lhs!(Mul, mul; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);
+make_delta_scalarop_lhs!(Div, div; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);
+
+macro_rules! make_delta_unary {
+    ($trait:ident, $fn:ident, $type:ty) => {
+        impl<'a, F: Num + Neg<Output = F>, Domain: Clone> Neg for $type {
+            type Output = DeltaMeasure<Domain, F>;
+            fn $fn(self) -> Self::Output {
+                let mut tmp = self.clone();
+                tmp.α = tmp.α.$fn();
+                tmp
+            }
+        }
+    };
+}
+
+make_delta_unary!(Neg, neg, DeltaMeasure<Domain, F>);
+make_delta_unary!(Neg, neg, &'a DeltaMeasure<Domain, F>);
+
+self_ownable!(DeltaMeasure<Domain, F> where Domain: Clone, F: Num);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/discrete.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,1143 @@
+//! This module implementes discrete measures.
+
+use super::base::*;
+use super::delta::*;
+use alg_tools::collection::Collection;
+use alg_tools::instance::{ClosedSpace, Decomposition, EitherDecomp, Instance, MyCow, Space};
+use alg_tools::iter::{MapF, Mappable};
+use alg_tools::linops::{Linear, Mapping};
+use alg_tools::loc::Loc;
+use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::norms::{Norm, Normed};
+use alg_tools::self_ownable;
+use alg_tools::tabledump::TableDump;
+use alg_tools::types::*;
+use nalgebra::DVector;
+use serde::ser::{Serialize, SerializeSeq, Serializer};
+use std::iter::Sum;
+use std::ops::{
+    Add, AddAssign, Div, DivAssign, Index, IndexMut, Mul, MulAssign, Neg, Sub, SubAssign,
+};
+
+/// Representation of a discrete measure.
+///
+/// This is the measure $μ = ∑_{k=1}^n α_k δ_{x_k}$, consisting of several
+/// [`DeltaMeasure`], i.e., “spikes” $α_k δ_{x_k}$ with weights $\alpha_k$ in `F` at locations
+/// $x_k$ in `Domain`.
+#[derive(Clone, Debug)]
+pub struct DiscreteMeasure<Domain, F: Num> {
+    pub(super) spikes: Vec<DeltaMeasure<Domain, F>>,
+}
+
+pub type RNDM<const N: usize, F = f64> = DiscreteMeasure<Loc<N, F>, F>;
+
+/// Iterator over the [`DeltaMeasure`] spikes of a [`DiscreteMeasure`].
+pub type SpikeIter<'a, Domain, F> = std::slice::Iter<'a, DeltaMeasure<Domain, F>>;
+
+/// Iterator over mutable [`DeltaMeasure`] spikes of a [`DiscreteMeasure`].
+pub type SpikeIterMut<'a, Domain, F> = std::slice::IterMut<'a, DeltaMeasure<Domain, F>>;
+
+/// Iterator over the locations of the spikes of a [`DiscreteMeasure`].
+pub type LocationIter<'a, Domain, F> =
+    std::iter::Map<SpikeIter<'a, Domain, F>, fn(&'a DeltaMeasure<Domain, F>) -> &'a Domain>;
+
+/// Iterator over the masses of the spikes of a [`DiscreteMeasure`].
+pub type MassIter<'a, Domain, F> =
+    std::iter::Map<SpikeIter<'a, Domain, F>, fn(&'a DeltaMeasure<Domain, F>) -> F>;
+
+/// Iterator over the mutable locations of the spikes of a [`DiscreteMeasure`].
+pub type MassIterMut<'a, Domain, F> = std::iter::Map<
+    SpikeIterMut<'a, Domain, F>,
+    for<'r> fn(&'r mut DeltaMeasure<Domain, F>) -> &'r mut F,
+>;
+
+impl<Domain, F: Num> DiscreteMeasure<Domain, F> {
+    /// Create a new zero measure (empty spike set).
+    pub fn new() -> Self {
+        DiscreteMeasure { spikes: Vec::new() }
+    }
+
+    /// Number of [`DeltaMeasure`] spikes in the measure
+    #[inline]
+    pub fn len(&self) -> usize {
+        self.spikes.len()
+    }
+
+    /// Replace with the zero measure.
+    #[inline]
+    pub fn clear(&mut self) {
+        self.spikes.clear()
+    }
+
+    /// Remove `i`:th spike, not maintaining order.
+    ///
+    /// Panics if indiex is out of bounds.
+    #[inline]
+    pub fn swap_remove(&mut self, i: usize) -> DeltaMeasure<Domain, F> {
+        self.spikes.swap_remove(i)
+    }
+
+    /// Iterate over (references to) the [`DeltaMeasure`] spikes in this measure
+    #[inline]
+    pub fn iter_spikes(&self) -> SpikeIter<'_, Domain, F> {
+        self.spikes.iter()
+    }
+
+    /// Iterate over mutable references to the [`DeltaMeasure`] spikes in this measure
+    #[inline]
+    pub fn iter_spikes_mut(&mut self) -> SpikeIterMut<'_, Domain, F> {
+        self.spikes.iter_mut()
+    }
+
+    /// Iterate over the location of the spikes in this measure
+    #[inline]
+    pub fn iter_locations(&self) -> LocationIter<'_, Domain, F> {
+        self.iter_spikes().map(DeltaMeasure::get_location)
+    }
+
+    /// Iterate over the masses of the spikes in this measure
+    #[inline]
+    pub fn iter_masses(&self) -> MassIter<'_, Domain, F> {
+        self.iter_spikes().map(DeltaMeasure::get_mass)
+    }
+
+    /// Iterate over the masses of the spikes in this measure
+    #[inline]
+    pub fn iter_masses_mut(&mut self) -> MassIterMut<'_, Domain, F> {
+        self.iter_spikes_mut().map(DeltaMeasure::get_mass_mut)
+    }
+
+    /// Update the masses of all the spikes to those produced by an iterator.
+    #[inline]
+    pub fn set_masses<I: Iterator<Item = F>>(&mut self, iter: I) {
+        self.spikes
+            .iter_mut()
+            .zip(iter)
+            .for_each(|(δ, α)| δ.set_mass(α));
+    }
+
+    /// Update the locations of all the spikes to those produced by an iterator.
+    #[inline]
+    pub fn set_locations<'a, I: Iterator<Item = &'a Domain>>(&mut self, iter: I)
+    where
+        Domain: 'static + Clone,
+    {
+        self.spikes
+            .iter_mut()
+            .zip(iter.cloned())
+            .for_each(|(δ, α)| δ.set_location(α));
+    }
+
+    // /// Map the masses of all the spikes using a function and an iterator
+    // #[inline]
+    // pub fn zipmap_masses<
+    //     I : Iterator<Item=F>,
+    //     G : Fn(F, I::Item) -> F
+    // > (&mut self, iter : I, g : G) {
+    //     self.spikes.iter_mut().zip(iter).for_each(|(δ, v)| δ.set_mass(g(δ.get_mass(), v)));
+    // }
+
+    /// Prune all spikes with zero mass.
+    #[inline]
+    pub fn prune(&mut self) {
+        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.
+    #[inline]
+    pub fn extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(&mut self, iter: I) {
+        self.spikes.extend(iter);
+    }
+
+    /// Add a spike to the measure
+    #[inline]
+    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> {
+    type Item = DeltaMeasure<Domain, F>;
+    type IntoIter = std::vec::IntoIter<DeltaMeasure<Domain, F>>;
+
+    #[inline]
+    fn into_iter(self) -> Self::IntoIter {
+        self.spikes.into_iter()
+    }
+}
+
+impl<'a, Domain, F: Num> IntoIterator for &'a DiscreteMeasure<Domain, F> {
+    type Item = &'a DeltaMeasure<Domain, F>;
+    type IntoIter = SpikeIter<'a, Domain, F>;
+
+    #[inline]
+    fn into_iter(self) -> Self::IntoIter {
+        self.spikes.iter()
+    }
+}
+
+impl<Domain, F: Num> Sum<DeltaMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    // Required method
+    fn sum<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = DeltaMeasure<Domain, F>>,
+    {
+        Self::from_iter(iter)
+    }
+}
+
+impl<'a, Domain: Clone, F: Num> Sum<&'a DeltaMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    // Required method
+    fn sum<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = &'a DeltaMeasure<Domain, F>>,
+    {
+        Self::from_iter(iter.cloned())
+    }
+}
+
+impl<Domain, F: Num> Sum<DiscreteMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    // Required method
+    fn sum<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = DiscreteMeasure<Domain, F>>,
+    {
+        Self::from_iter(iter.map(|μ| μ.into_iter()).flatten())
+    }
+}
+
+impl<'a, Domain: Clone, F: Num> Sum<&'a DiscreteMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    // Required method
+    fn sum<I>(iter: I) -> Self
+    where
+        I: Iterator<Item = &'a DiscreteMeasure<Domain, F>>,
+    {
+        Self::from_iter(iter.map(|μ| μ.iter_spikes()).flatten().cloned())
+    }
+}
+
+impl<Domain: Clone, F: Float> DiscreteMeasure<Domain, F> {
+    /// Computes `μ1 ← θ * μ1 - ζ * μ2`, pruning entries where both `μ1` (`self`) and `μ2` have
+    // zero weight. `μ2` will contain a pruned copy of pruned original `μ1` without arithmetic
+    /// performed. **This expects `self` and `μ2` to have matching coordinates in each index**.
+    // `μ2` can be than `self`, but not longer.
+    pub fn pruning_sub(&mut self, θ: F, ζ: F, μ2: &mut Self) {
+        for δ in &self[μ2.len()..] {
+            μ2.push(DeltaMeasure { x: δ.x.clone(), α: F::ZERO });
+        }
+        debug_assert_eq!(self.len(), μ2.len());
+        let mut dest = 0;
+        for i in 0..self.len() {
+            let α = self[i].α;
+            let α_new = θ * α - ζ * μ2[i].α;
+            if dest < i {
+                μ2[dest] = DeltaMeasure { x: self[i].x.clone(), α };
+                self[dest] = DeltaMeasure { x: self[i].x.clone(), α: α_new };
+            } else {
+                μ2[i].α = α;
+                self[i].α = α_new;
+            }
+            dest += 1;
+        }
+        self.spikes.truncate(dest);
+        μ2.spikes.truncate(dest);
+    }
+}
+
+impl<Domain, F: Float> DiscreteMeasure<Domain, F> {
+    /// Prune all spikes with mass absolute value less than the given `tolerance`.
+    #[inline]
+    pub fn prune_approx(&mut self, tolerance: F) {
+        self.spikes.retain(|δ| δ.α.abs() > tolerance);
+    }
+}
+
+impl<Domain, F: Float + ToNalgebraRealField> DiscreteMeasure<Domain, F> {
+    /// Extracts the masses of the spikes as a [`DVector`].
+    pub fn masses_dvector(&self) -> DVector<F::MixedType> {
+        DVector::from_iterator(
+            self.len(),
+            self.iter_masses().map(|α| α.to_nalgebra_mixed()),
+        )
+    }
+
+    /// Sets the masses of the spikes from the values of a [`DVector`].
+    pub fn set_masses_dvector(&mut self, x: &DVector<F::MixedType>) {
+        self.set_masses(x.iter().map(|&α| F::from_nalgebra_mixed(α)));
+    }
+
+    // /// Extracts the masses of the spikes as a [`Vec`].
+    // pub fn masses_vec(&self) -> Vec<F::MixedType> {
+    //     self.iter_masses()
+    //         .map(|α| α.to_nalgebra_mixed())
+    //         .collect()
+    // }
+
+    // /// Sets the masses of the spikes from the values of a [`Vec`].
+    // pub fn set_masses_vec(&mut self, x : &Vec<F::MixedType>) {
+    //     self.set_masses(x.iter().map(|&α| F::from_nalgebra_mixed(α)));
+    // }
+}
+
+// impl<Domain, F :Num> Index<usize> for DiscreteMeasure<Domain, F> {
+//     type Output = DeltaMeasure<Domain, F>;
+//     #[inline]
+//     fn index(&self, i : usize) -> &Self::Output {
+//         self.spikes.index(i)
+//     }
+// }
+
+// impl<Domain, F :Num> IndexMut<usize> for DiscreteMeasure<Domain, F> {
+//     #[inline]
+//     fn index_mut(&mut self, i : usize) -> &mut Self::Output {
+//         self.spikes.index_mut(i)
+//     }
+// }
+
+impl<Domain, F: Num, I: std::slice::SliceIndex<[DeltaMeasure<Domain, F>]>> Index<I>
+    for DiscreteMeasure<Domain, F>
+{
+    type Output = <I as std::slice::SliceIndex<[DeltaMeasure<Domain, F>]>>::Output;
+    #[inline]
+    fn index(&self, i: I) -> &Self::Output {
+        self.spikes.index(i)
+    }
+}
+
+impl<Domain, F: Num, I: std::slice::SliceIndex<[DeltaMeasure<Domain, F>]>> IndexMut<I>
+    for DiscreteMeasure<Domain, F>
+{
+    #[inline]
+    fn index_mut(&mut self, i: I) -> &mut Self::Output {
+        self.spikes.index_mut(i)
+    }
+}
+
+impl<Domain, F: Num, D: Into<DeltaMeasure<Domain, F>>, const K: usize> From<[D; K]>
+    for DiscreteMeasure<Domain, F>
+{
+    #[inline]
+    fn from(list: [D; K]) -> Self {
+        list.into_iter().collect()
+    }
+}
+
+impl<Domain, F: Num> From<Vec<DeltaMeasure<Domain, F>>> for DiscreteMeasure<Domain, F> {
+    #[inline]
+    fn from(spikes: Vec<DeltaMeasure<Domain, F>>) -> Self {
+        DiscreteMeasure { spikes }
+    }
+}
+
+impl<'a, Domain, F: Num, D> From<&'a [D]> for DiscreteMeasure<Domain, F>
+where
+    &'a D: Into<DeltaMeasure<Domain, F>>,
+{
+    #[inline]
+    fn from(list: &'a [D]) -> Self {
+        list.into_iter().map(|d| d.into()).collect()
+    }
+}
+
+impl<Domain, F: Num> From<DeltaMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    #[inline]
+    fn from(δ: DeltaMeasure<Domain, F>) -> Self {
+        DiscreteMeasure { spikes: vec![δ] }
+    }
+}
+
+impl<'a, Domain: Clone, F: Num> From<&'a DeltaMeasure<Domain, F>> for DiscreteMeasure<Domain, F> {
+    #[inline]
+    fn from(δ: &'a DeltaMeasure<Domain, F>) -> Self {
+        DiscreteMeasure { spikes: vec![δ.clone()] }
+    }
+}
+
+impl<Domain, F: Num, D: Into<DeltaMeasure<Domain, F>>> FromIterator<D>
+    for DiscreteMeasure<Domain, F>
+{
+    #[inline]
+    fn from_iter<T>(iter: T) -> Self
+    where
+        T: IntoIterator<Item = D>,
+    {
+        DiscreteMeasure { spikes: iter.into_iter().map(|m| m.into()).collect() }
+    }
+}
+
+impl<'a, F: Num, const N: usize> TableDump<'a> for DiscreteMeasure<Loc<N, F>, F>
+where
+    DeltaMeasure<Loc<N, F>, F>: Serialize + 'a,
+{
+    type Iter = std::slice::Iter<'a, DeltaMeasure<Loc<N, F>, F>>;
+
+    // fn tabledump_headers(&'a self) -> Vec<String> {
+    //     let mut v : Vec<String> = (0..N).map(|i| format!("x{}", i)).collect();
+    //     v.push("weight".into());
+    //     v
+    // }
+
+    fn tabledump_entries(&'a self) -> Self::Iter {
+        // Ensure order matching the headers above
+        self.spikes.iter()
+    }
+}
+
+// Need to manually implement serialisation for DeltaMeasure<Loc<N, F>, F> [`csv`] writer fails on
+// structs with nested arrays as well as with #[serde(flatten)].
+// Then derive no longer works for DiscreteMeasure
+impl<F: Num, const N: usize> Serialize for DiscreteMeasure<Loc<N, F>, F>
+where
+    F: Serialize,
+{
+    fn serialize<S>(&self, serializer: S) -> Result<S::Ok, S::Error>
+    where
+        S: Serializer,
+    {
+        let mut s = serializer.serialize_seq(Some(self.spikes.len()))?;
+        for δ in self.spikes.iter() {
+            s.serialize_element(δ)?;
+        }
+        s.end()
+    }
+}
+
+impl<Domain: PartialEq, F: Float> Measure<F> for DiscreteMeasure<Domain, F> {
+    type Domain = Domain;
+}
+
+impl<Domain: PartialEq, F: Float> Norm<Radon, F> for DiscreteMeasure<Domain, F>
+where
+    DeltaMeasure<Domain, F>: Norm<Radon, F>,
+{
+    #[inline]
+    fn norm(&self, _: Radon) -> F {
+        self.spikes.iter().map(|m| m.norm(Radon)).sum()
+    }
+}
+
+impl<Domain, G, F: Num> Mapping<G> for DiscreteMeasure<Domain, F>
+where
+    Domain: Space,
+    G::Codomain: Sum + Mul<F, Output = G::Codomain>,
+    G: Mapping<Domain, Codomain = F> + Clone + ClosedSpace,
+    for<'b> &'b Domain: Instance<Domain>,
+{
+    type Codomain = G::Codomain;
+
+    #[inline]
+    fn apply<I: Instance<G>>(&self, g: I) -> Self::Codomain {
+        g.eval(|g| self.spikes.iter().map(|m| g.apply(&m.x) * m.α).sum())
+    }
+}
+
+impl<Domain, G, F: Num> Linear<G> for DiscreteMeasure<Domain, F>
+where
+    Domain: Space,
+    G::Codomain: Sum + Mul<F, Output = G::Codomain>,
+    G: Mapping<Domain, Codomain = F> + Clone + ClosedSpace,
+    for<'b> &'b Domain: Instance<Domain>,
+{
+}
+
+/// Helper trait for constructing arithmetic operations for combinations
+/// of [`DiscreteMeasure`] and [`DeltaMeasure`], and their references.
+trait Lift<F: Num, Domain> {
+    type Producer: Iterator<Item = DeltaMeasure<Domain, F>>;
+
+    #[allow(dead_code)]
+    /// Lifts `self` into a [`DiscreteMeasure`].
+    fn lift(self) -> DiscreteMeasure<Domain, F>;
+
+    /// Lifts `self` into a [`DiscreteMeasure`], apply either `f` or `f_mut` whether the type
+    /// this method is implemented for is a reference or or not.
+    fn lift_with(
+        self,
+        f: impl Fn(&DeltaMeasure<Domain, F>) -> DeltaMeasure<Domain, F>,
+        f_mut: impl FnMut(&mut DeltaMeasure<Domain, F>),
+    ) -> DiscreteMeasure<Domain, F>;
+
+    /// Extend `self` into a [`DiscreteMeasure`] with the spikes produced by `iter`.
+    fn lift_extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(
+        self,
+        iter: I,
+    ) -> DiscreteMeasure<Domain, F>;
+
+    /// Returns an iterator for producing copies of the spikes of `self`.
+    fn produce(self) -> Self::Producer;
+}
+
+impl<F: Num, Domain> Lift<F, Domain> for DiscreteMeasure<Domain, F> {
+    type Producer = std::vec::IntoIter<DeltaMeasure<Domain, F>>;
+
+    #[inline]
+    fn lift(self) -> DiscreteMeasure<Domain, F> {
+        self
+    }
+
+    fn lift_with(
+        mut self,
+        _f: impl Fn(&DeltaMeasure<Domain, F>) -> DeltaMeasure<Domain, F>,
+        f_mut: impl FnMut(&mut DeltaMeasure<Domain, F>),
+    ) -> DiscreteMeasure<Domain, F> {
+        self.spikes.iter_mut().for_each(f_mut);
+        self
+    }
+
+    #[inline]
+    fn lift_extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(
+        mut self,
+        iter: I,
+    ) -> DiscreteMeasure<Domain, F> {
+        self.spikes.extend(iter);
+        self
+    }
+
+    #[inline]
+    fn produce(self) -> Self::Producer {
+        self.spikes.into_iter()
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Lift<F, Domain> for &'a DiscreteMeasure<Domain, F> {
+    type Producer = MapF<std::slice::Iter<'a, DeltaMeasure<Domain, F>>, DeltaMeasure<Domain, F>>;
+
+    #[inline]
+    fn lift(self) -> DiscreteMeasure<Domain, F> {
+        self.clone()
+    }
+
+    fn lift_with(
+        self,
+        f: impl Fn(&DeltaMeasure<Domain, F>) -> DeltaMeasure<Domain, F>,
+        _f_mut: impl FnMut(&mut DeltaMeasure<Domain, F>),
+    ) -> DiscreteMeasure<Domain, F> {
+        DiscreteMeasure { spikes: self.spikes.iter().map(f).collect() }
+    }
+
+    #[inline]
+    fn lift_extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(
+        self,
+        iter: I,
+    ) -> DiscreteMeasure<Domain, F> {
+        let mut res = self.clone();
+        res.spikes.extend(iter);
+        res
+    }
+
+    #[inline]
+    fn produce(self) -> Self::Producer {
+        // TODO: maybe not optimal to clone here and would benefit from
+        // a reference version of lift_extend.
+        self.spikes.iter().mapF(Clone::clone)
+    }
+}
+
+impl<F: Num, Domain> Lift<F, Domain> for DeltaMeasure<Domain, F> {
+    type Producer = std::iter::Once<DeltaMeasure<Domain, F>>;
+
+    #[inline]
+    fn lift(self) -> DiscreteMeasure<Domain, F> {
+        DiscreteMeasure { spikes: vec![self] }
+    }
+
+    #[inline]
+    fn lift_with(
+        mut self,
+        _f: impl Fn(&DeltaMeasure<Domain, F>) -> DeltaMeasure<Domain, F>,
+        mut f_mut: impl FnMut(&mut DeltaMeasure<Domain, F>),
+    ) -> DiscreteMeasure<Domain, F> {
+        f_mut(&mut self);
+        DiscreteMeasure { spikes: vec![self] }
+    }
+
+    #[inline]
+    fn lift_extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(
+        self,
+        iter: I,
+    ) -> DiscreteMeasure<Domain, F> {
+        let mut spikes = vec![self];
+        spikes.extend(iter);
+        DiscreteMeasure { spikes: spikes }
+    }
+
+    #[inline]
+    fn produce(self) -> Self::Producer {
+        std::iter::once(self)
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Lift<F, Domain> for &'a DeltaMeasure<Domain, F> {
+    type Producer = std::iter::Once<DeltaMeasure<Domain, F>>;
+
+    #[inline]
+    fn lift(self) -> DiscreteMeasure<Domain, F> {
+        DiscreteMeasure { spikes: vec![self.clone()] }
+    }
+
+    #[inline]
+    fn lift_with(
+        self,
+        f: impl Fn(&DeltaMeasure<Domain, F>) -> DeltaMeasure<Domain, F>,
+        _f_mut: impl FnMut(&mut DeltaMeasure<Domain, F>),
+    ) -> DiscreteMeasure<Domain, F> {
+        DiscreteMeasure { spikes: vec![f(self)] }
+    }
+
+    #[inline]
+    fn lift_extend<I: Iterator<Item = DeltaMeasure<Domain, F>>>(
+        self,
+        iter: I,
+    ) -> DiscreteMeasure<Domain, F> {
+        let mut spikes = vec![self.clone()];
+        spikes.extend(iter);
+        DiscreteMeasure { spikes: spikes }
+    }
+
+    #[inline]
+    fn produce(self) -> Self::Producer {
+        std::iter::once(self.clone())
+    }
+}
+
+macro_rules! make_discrete_addsub_assign {
+    ($rhs:ty) => {
+        // Discrete += (&)Discrete
+        impl<'a, F: Num, Domain: Clone> AddAssign<$rhs> for DiscreteMeasure<Domain, F> {
+            fn add_assign(&mut self, other: $rhs) {
+                self.spikes.extend(other.produce());
+            }
+        }
+
+        impl<'a, F: Num + Neg<Output = F>, Domain: Clone> SubAssign<$rhs>
+            for DiscreteMeasure<Domain, F>
+        {
+            fn sub_assign(&mut self, other: $rhs) {
+                self.spikes.extend(other.produce().map(|δ| -δ));
+            }
+        }
+    };
+}
+
+make_discrete_addsub_assign!(DiscreteMeasure<Domain, F>);
+make_discrete_addsub_assign!(&'a DiscreteMeasure<Domain, F>);
+make_discrete_addsub_assign!(DeltaMeasure<Domain, F>);
+make_discrete_addsub_assign!(&'a DeltaMeasure<Domain, F>);
+
+macro_rules! make_discrete_addsub {
+    ($lhs:ty, $rhs:ty, $alt_order:expr) => {
+        impl<'a, 'b, F: Num, Domain: Clone> Add<$rhs> for $lhs {
+            type Output = DiscreteMeasure<Domain, F>;
+            fn add(self, other: $rhs) -> DiscreteMeasure<Domain, F> {
+                if !$alt_order {
+                    self.lift_extend(other.produce())
+                } else {
+                    other.lift_extend(self.produce())
+                }
+            }
+        }
+
+        impl<'a, 'b, F: Num + Neg<Output = F>, Domain: Clone> Sub<$rhs> for $lhs {
+            type Output = DiscreteMeasure<Domain, F>;
+            fn sub(self, other: $rhs) -> DiscreteMeasure<Domain, F> {
+                self.lift_extend(other.produce().map(|δ| -δ))
+            }
+        }
+    };
+}
+
+make_discrete_addsub!(DiscreteMeasure<Domain, F>,     DiscreteMeasure<Domain, F>,     false);
+make_discrete_addsub!(DiscreteMeasure<Domain, F>,     &'b DiscreteMeasure<Domain, F>, false);
+make_discrete_addsub!(&'a DiscreteMeasure<Domain, F>, DiscreteMeasure<Domain, F>,     true);
+make_discrete_addsub!(
+    &'a DiscreteMeasure<Domain, F>,
+    &'b DiscreteMeasure<Domain, F>,
+    false
+);
+make_discrete_addsub!(DeltaMeasure<Domain, F>,        DiscreteMeasure<Domain, F>,     false);
+make_discrete_addsub!(DeltaMeasure<Domain, F>,        &'b DiscreteMeasure<Domain, F>, false);
+make_discrete_addsub!(&'a DeltaMeasure<Domain, F>,    DiscreteMeasure<Domain, F>,     true);
+make_discrete_addsub!(
+    &'a DeltaMeasure<Domain, F>,
+    &'b DiscreteMeasure<Domain, F>,
+    false
+);
+make_discrete_addsub!(DiscreteMeasure<Domain, F>,     DeltaMeasure<Domain, F>,        false);
+make_discrete_addsub!(DiscreteMeasure<Domain, F>,     &'b DeltaMeasure<Domain, F>,    false);
+make_discrete_addsub!(&'a DiscreteMeasure<Domain, F>, DeltaMeasure<Domain, F>,        false);
+make_discrete_addsub!(
+    &'a DiscreteMeasure<Domain, F>,
+    &'b DeltaMeasure<Domain, F>,
+    false
+);
+make_discrete_addsub!(DeltaMeasure<Domain, F>,        DeltaMeasure<Domain, F>,        false);
+make_discrete_addsub!(DeltaMeasure<Domain, F>,        &'b DeltaMeasure<Domain, F>,    false);
+make_discrete_addsub!(&'a DeltaMeasure<Domain, F>,    DeltaMeasure<Domain, F>,        false);
+make_discrete_addsub!(
+    &'a DeltaMeasure<Domain, F>,
+    &'b DeltaMeasure<Domain, F>,
+    false
+);
+
+macro_rules! make_discrete_scalarop_rhs {
+    ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
+        make_discrete_scalarop_rhs!(@assign DiscreteMeasure<Domain, F>, F, $trait_assign, $fn_assign);
+        make_discrete_scalarop_rhs!(@assign DiscreteMeasure<Domain, F>, &'a F, $trait_assign, $fn_assign);
+        make_discrete_scalarop_rhs!(@new DiscreteMeasure<Domain, F>, F, $trait, $fn, $fn_assign);
+        make_discrete_scalarop_rhs!(@new DiscreteMeasure<Domain, F>, &'a F, $trait, $fn, $fn_assign);
+        make_discrete_scalarop_rhs!(@new &'b DiscreteMeasure<Domain, F>, F, $trait, $fn, $fn_assign);
+        make_discrete_scalarop_rhs!(@new &'b DiscreteMeasure<Domain, F>, &'a F, $trait, $fn, $fn_assign);
+    };
+
+    (@assign $lhs:ty, $rhs:ty, $trait_assign:ident, $fn_assign:ident) => {
+        impl<'a, 'b, F : Num, Domain> $trait_assign<$rhs> for $lhs {
+            fn $fn_assign(&mut self, b : $rhs) {
+                self.spikes.iter_mut().for_each(|δ| δ.$fn_assign(b));
+            }
+        }
+    };
+    (@new $lhs:ty, $rhs:ty, $trait:ident, $fn:ident, $fn_assign:ident) => {
+        impl<'a, 'b, F : Num, Domain : Clone> $trait<$rhs> for $lhs {
+            type Output = DiscreteMeasure<Domain, F>;
+            fn $fn(self, b : $rhs) -> Self::Output {
+                self.lift_with(|δ| δ.$fn(b), |δ| δ.$fn_assign(b))
+            }
+        }
+    };
+}
+
+make_discrete_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
+make_discrete_scalarop_rhs!(Div, div, DivAssign, div_assign);
+
+macro_rules! make_discrete_unary {
+    ($trait:ident, $fn:ident, $type:ty) => {
+        impl<'a, F: Num + Neg<Output = F>, Domain: Clone> Neg for $type {
+            type Output = DiscreteMeasure<Domain, F>;
+            fn $fn(self) -> Self::Output {
+                self.lift_with(|δ| δ.$fn(), |δ| δ.α = δ.α.$fn())
+            }
+        }
+    };
+}
+
+make_discrete_unary!(Neg, neg, DiscreteMeasure<Domain, F>);
+make_discrete_unary!(Neg, neg, &'a DiscreteMeasure<Domain, F>);
+
+// impl<F : Num, Domain> Neg for DiscreteMeasure<Domain, F> {
+//     type Output = Self;
+//     fn $fn(mut self, b : F) -> Self {
+//         self.lift().spikes.iter_mut().for_each(|δ| δ.neg(b));
+//         self
+//     }
+// }
+
+macro_rules! make_discrete_scalarop_lhs {
+    ($trait:ident, $fn:ident; $($f:ident)+) => { $(
+        impl<Domain> $trait<DiscreteMeasure<Domain, $f>> for $f {
+            type Output = DiscreteMeasure<Domain, $f>;
+            fn $fn(self, mut v : DiscreteMeasure<Domain, $f>) -> Self::Output {
+                v.spikes.iter_mut().for_each(|δ| δ.α = self.$fn(δ.α));
+                v
+            }
+        }
+
+        impl<'a, Domain : Copy> $trait<&'a DiscreteMeasure<Domain, $f>> for $f {
+            type Output = DiscreteMeasure<Domain, $f>;
+            fn $fn(self, v : &'a DiscreteMeasure<Domain, $f>) -> Self::Output {
+                DiscreteMeasure{
+                    spikes : v.spikes.iter().map(|δ| self.$fn(δ)).collect()
+                }
+            }
+        }
+
+        impl<'b, Domain> $trait<DiscreteMeasure<Domain, $f>> for &'b $f {
+            type Output = DiscreteMeasure<Domain, $f>;
+            fn $fn(self, mut v : DiscreteMeasure<Domain, $f>) -> Self::Output {
+                v.spikes.iter_mut().for_each(|δ| δ.α = self.$fn(δ.α));
+                v
+            }
+        }
+
+        impl<'a, 'b, Domain : Copy> $trait<&'a DiscreteMeasure<Domain, $f>> for &'b $f {
+            type Output = DiscreteMeasure<Domain, $f>;
+            fn $fn(self, v : &'a DiscreteMeasure<Domain, $f>) -> Self::Output {
+                DiscreteMeasure{
+                    spikes : v.spikes.iter().map(|δ| self.$fn(δ)).collect()
+                }
+            }
+        }
+    )+ }
+}
+
+make_discrete_scalarop_lhs!(Mul, mul; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);
+make_discrete_scalarop_lhs!(Div, div; f32 f64 i8 i16 i32 i64 isize u8 u16 u32 u64 usize);
+
+impl<F: Num, Domain> Collection for DiscreteMeasure<Domain, F> {
+    type Element = DeltaMeasure<Domain, F>;
+    type RefsIter<'a>
+        = std::slice::Iter<'a, Self::Element>
+    where
+        Self: 'a;
+
+    #[inline]
+    fn iter_refs(&self) -> Self::RefsIter<'_> {
+        self.iter_spikes()
+    }
+}
+
+impl<Domain: Clone, F: Num> Space for DiscreteMeasure<Domain, F> {
+    type Principal = Self;
+
+    type Decomp = MeasureDecomp;
+}
+
+pub type SpikeSlice<'b, Domain, F> = &'b [DeltaMeasure<Domain, F>];
+
+pub type EitherSlice<'b, Domain, F> =
+    EitherDecomp<Vec<DeltaMeasure<Domain, F>>, SpikeSlice<'b, Domain, F>>;
+
+impl<F: Num, Domain: Clone> Decomposition<DiscreteMeasure<Domain, F>> for MeasureDecomp {
+    type Decomposition<'b>
+        = EitherSlice<'b, Domain, F>
+    where
+        DiscreteMeasure<Domain, F>: 'b;
+    type Reference<'b>
+        = SpikeSlice<'b, Domain, F>
+    where
+        DiscreteMeasure<Domain, F>: 'b;
+
+    /// Left the lightweight reference type into a full decomposition type.
+    fn lift<'b>(r: Self::Reference<'b>) -> Self::Decomposition<'b> {
+        EitherDecomp::Borrowed(r)
+    }
+}
+
+impl<F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for DiscreteMeasure<Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(self.spikes.as_slice())
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self)
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        self
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Owned(self.spikes)
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for &'a DiscreteMeasure<Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(self.spikes.as_slice())
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Borrowed(self)
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        self.clone()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Borrowed(self.spikes.as_slice())
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for EitherSlice<'a, Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        match self {
+            EitherDecomp::Owned(v) => f(v.as_slice()),
+            EitherDecomp::Borrowed(s) => f(s),
+        }
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        match self {
+            EitherDecomp::Owned(v) => v.into(),
+            EitherDecomp::Borrowed(s) => s.into(),
+        }
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        self
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for &'a EitherSlice<'a, Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        match self {
+            EitherDecomp::Owned(v) => f(v.as_slice()),
+            EitherDecomp::Borrowed(s) => f(s),
+        }
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        match self {
+            EitherDecomp::Owned(v) => v.as_slice(),
+            EitherDecomp::Borrowed(s) => s,
+        }
+        .into()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        match self {
+            EitherDecomp::Owned(v) => EitherDecomp::Borrowed(v.as_slice()),
+            EitherDecomp::Borrowed(s) => EitherDecomp::Borrowed(s),
+        }
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for SpikeSlice<'a, Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(self)
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        self.into()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Borrowed(self)
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for &'a SpikeSlice<'a, Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(*self)
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        (*self).into()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Borrowed(*self)
+    }
+}
+
+impl<F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for DeltaMeasure<Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(std::slice::from_ref(self))
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        self.into()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Owned(vec![self])
+    }
+}
+
+impl<'a, F: Num, Domain: Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
+    for &'a DeltaMeasure<Domain, F>
+{
+    fn eval_ref<'b, R>(&'b self, f: impl FnOnce(SpikeSlice<'b, Domain, F>) -> R) -> R
+    where
+        DiscreteMeasure<Domain, F>: 'b,
+        Self: 'b,
+    {
+        f(std::slice::from_ref(*self))
+    }
+
+    fn cow<'b>(self) -> MyCow<'b, DiscreteMeasure<Domain, F>>
+    where
+        Self: 'b,
+    {
+        MyCow::Owned(self.own())
+    }
+
+    fn own(self) -> DiscreteMeasure<Domain, F> {
+        self.into()
+    }
+
+    fn decompose<'b>(self) -> EitherSlice<'b, Domain, F>
+    where
+        Self: 'b,
+    {
+        EitherDecomp::Borrowed(std::slice::from_ref(self))
+    }
+}
+
+impl<F, Domain> Normed<F> for DiscreteMeasure<Domain, F>
+where
+    F: Float,
+    Domain: Clone + PartialEq,
+    DeltaMeasure<Domain, F>: Norm<Radon, F>,
+{
+    type NormExp = Radon;
+
+    fn norm_exponent(&self) -> Radon {
+        Radon
+    }
+}
+
+self_ownable!(DiscreteMeasure<Domain, F> where Domain: Clone, F: Num);

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/lib.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,18 @@
+// The main documentation is in the README.
+#![doc = include_str!("../README.md")]
+// We use unicode. We would like to use much more of it than Rust allows.
+// Live with it. Embrace it.
+#![allow(uncommon_codepoints)]
+#![allow(mixed_script_confusables)]
+#![allow(confusable_idents)]
+
+mod base;
+pub use base::*;
+mod delta;
+pub use delta::*;
+mod discrete;
+pub use discrete::*;
+pub mod merging;
+
+#[cfg(feature = "pyo3")]
+pub mod python;

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/merging.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,316 @@
+/*!
+Spike merging heuristics for [`DiscreteMeasure`]s.
+
+This module primarily provides the [`SpikeMerging`] trait, and within it,
+the [`SpikeMerging::merge_spikes`] method. The trait is implemented on
+[`DiscreteMeasure<Loc<N, F>, F>`]s in dimensions `N=1` and `N=2`.
+*/
+
+use numeric_literals::replace_float_literals;
+use serde::{Deserialize, Serialize};
+use std::cmp::Ordering;
+//use clap::builder::{PossibleValuesParser, PossibleValue};
+use alg_tools::nanleast::NaNLeast;
+
+use super::delta::*;
+use super::discrete::*;
+use alg_tools::loc::Loc;
+use alg_tools::types::*;
+
+/// Spike merging heuristic selection
+#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
+#[allow(dead_code)]
+pub struct SpikeMergingMethod<F> {
+    // Merging radius
+    pub radius: F,
+    // Enabled
+    pub enabled: bool,
+    // Interpolate merged points
+    pub interp: bool,
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float> Default for SpikeMergingMethod<F> {
+    fn default() -> Self {
+        SpikeMergingMethod { radius: 0.01, enabled: false, interp: true }
+    }
+}
+
+/// Trait for dimension-dependent implementation of heuristic peak merging strategies.
+pub trait SpikeMerging<F> {
+    /// Attempt spike merging according to [`SpikeMerging`] method.
+    ///
+    /// Returns the last [`Some`] returned by the merging candidate acceptance decision closure
+    /// `accept` if any merging is performed. The closure should accept as its only parameter a
+    /// new candidate measure (it will generally be internally mutated `self`, although this is
+    /// not guaranteed), and return [`None`] if the merge is accepted, and otherwise a [`Some`] of
+    /// an arbitrary value. This method will return that value for the *last* accepted merge, or
+    /// [`None`] if no merge was accepted.
+    ///
+    /// This method is stable with respect to spike locations:  on merge, the weights of existing
+    /// removed spikes is set to zero, new ones inserted at the end of the spike vector.
+    /// They merge may also be performed by increasing the weights of the existing spikes,
+    /// without inserting new spikes.
+    fn merge_spikes<G>(&mut self, method: SpikeMergingMethod<F>, accept: G) -> usize
+    where
+        G: FnMut(&'_ Self) -> bool,
+    {
+        if method.enabled {
+            self.do_merge_spikes_radius(method.radius, method.interp, accept)
+        } else {
+            0
+        }
+    }
+
+    /// Attempt to merge spikes based on a value and a fitness function.
+    ///
+    /// Calls [`SpikeMerging::merge_spikes`] with `accept` constructed from the composition of
+    /// `value` and `fitness`, compared to initial fitness. Returns the last return value of `value`
+    // for a merge  accepted by `fitness`. If no merge was accepted, `value` applied to the initial
+    /// `self` is returned. also the number of merges is returned;
+    fn merge_spikes_fitness<G, H, V, O>(
+        &mut self,
+        method: SpikeMergingMethod<F>,
+        value: G,
+        fitness: H,
+    ) -> (V, usize)
+    where
+        G: Fn(&'_ Self) -> V,
+        H: Fn(&'_ V) -> O,
+        O: PartialOrd,
+    {
+        let mut res = value(self);
+        let initial_fitness = fitness(&res);
+        let count = self.merge_spikes(method, |μ| {
+            res = value(μ);
+            fitness(&res) <= initial_fitness
+        });
+        (res, count)
+    }
+
+    /// Attempt to merge spikes that are within radius $ρ$ of each other (unspecified norm).
+    ///
+    /// This method implements [`SpikeMerging::merge_spikes`].
+    fn do_merge_spikes_radius<G>(&mut self, ρ: F, interp: bool, accept: G) -> usize
+    where
+        G: FnMut(&'_ Self) -> bool;
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float, const N: usize> DiscreteMeasure<Loc<N, F>, F> {
+    /// Attempts to merge spikes with indices `i` and `j`.
+    ///
+    /// This assumes that the weights of the two spikes have already been checked not to be zero.
+    ///
+    /// The parameter `res` points to the current “result” for [`SpikeMerging::merge_spikes`].
+    /// If the merge is accepted by `accept` returning a [`Some`], `res` will be replaced by its
+    /// return value.
+    ///
+    /// Returns the index of `self.spikes` storing the new spike.
+    fn attempt_merge<G>(
+        &mut self,
+        i: usize,
+        j: usize,
+        interp: bool,
+        accept: &mut G,
+    ) -> Option<usize>
+    where
+        G: FnMut(&'_ Self) -> bool,
+    {
+        let &DeltaMeasure { x: xi, α: αi } = &self.spikes[i];
+        let &DeltaMeasure { x: xj, α: αj } = &self.spikes[j];
+
+        if interp {
+            // Merge inplace
+            self.spikes[i].α = 0.0;
+            self.spikes[j].α = 0.0;
+            let αia = αi.abs();
+            let αja = αj.abs();
+            self.spikes
+                .push(DeltaMeasure { α: αi + αj, x: (xi * αia + xj * αja) / (αia + αja) });
+            if accept(self) {
+                Some(self.spikes.len() - 1)
+            } else {
+                // Merge not accepted, restore modification
+                self.spikes[i].α = αi;
+                self.spikes[j].α = αj;
+                self.spikes.pop();
+                None
+            }
+        } else {
+            // Attempt merge inplace, first combination
+            self.spikes[i].α = αi + αj;
+            self.spikes[j].α = 0.0;
+            if accept(self) {
+                // Merge accepted
+                Some(i)
+            } else {
+                // Attempt merge inplace, second combination
+                self.spikes[i].α = 0.0;
+                self.spikes[j].α = αi + αj;
+                if accept(self) {
+                    // Merge accepted
+                    Some(j)
+                } else {
+                    // Merge not accepted, restore modification
+                    self.spikes[i].α = αi;
+                    self.spikes[j].α = αj;
+                    None
+                }
+            }
+        }
+    }
+}
+
+/// Sorts a vector of indices into `slice` by `compare`.
+///
+/// The closure `compare` operators on references to elements of `slice`.
+/// Returns the sorted vector of indices into `slice`.
+pub fn sort_indices_by<V, F>(slice: &[V], mut compare: F) -> Vec<usize>
+where
+    F: FnMut(&V, &V) -> Ordering,
+{
+    let mut indices = Vec::from_iter(0..slice.len());
+    indices.sort_by(|&i, &j| compare(&slice[i], &slice[j]));
+    indices
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float> SpikeMerging<F> for DiscreteMeasure<Loc<1, F>, F> {
+    fn do_merge_spikes_radius<G>(&mut self, ρ: F, interp: bool, mut accept: G) -> usize
+    where
+        G: FnMut(&'_ Self) -> bool,
+    {
+        // Sort by coordinate into an indexing array.
+        let mut indices = sort_indices_by(&self.spikes, |&δ1, &δ2| {
+            let &Loc([x1]) = &δ1.x;
+            let &Loc([x2]) = &δ2.x;
+            // nan-ignoring ordering of floats
+            NaNLeast(x1).cmp(&NaNLeast(x2))
+        });
+
+        // Initialise result
+        let mut count = 0;
+
+        // Scan consecutive pairs and merge if close enough and accepted by `accept`.
+        if indices.len() == 0 {
+            return count;
+        }
+        for k in 0..(indices.len() - 1) {
+            let i = indices[k];
+            let j = indices[k + 1];
+            let &DeltaMeasure { x: Loc([xi]), α: αi } = &self.spikes[i];
+            let &DeltaMeasure { x: Loc([xj]), α: αj } = &self.spikes[j];
+            debug_assert!(xi <= xj);
+            // If close enough, attempt merging
+            if αi != 0.0 && αj != 0.0 && xj <= xi + ρ {
+                if let Some(l) = self.attempt_merge(i, j, interp, &mut accept) {
+                    // For this to work (the debug_assert! to not trigger above), the new
+                    // coordinate produced by attempt_merge has to be at most xj.
+                    indices[k + 1] = l;
+                    count += 1
+                }
+            }
+        }
+
+        count
+    }
+}
+
+/// Orders `δ1` and `δ1` according to the first coordinate.
+fn compare_first_coordinate<F: Float>(
+    δ1: &DeltaMeasure<Loc<2, F>, F>,
+    δ2: &DeltaMeasure<Loc<2, F>, F>,
+) -> Ordering {
+    let &Loc([x11, ..]) = &δ1.x;
+    let &Loc([x21, ..]) = &δ2.x;
+    // nan-ignoring ordering of floats
+    NaNLeast(x11).cmp(&NaNLeast(x21))
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<F: Float> SpikeMerging<F> for DiscreteMeasure<Loc<2, F>, F> {
+    fn do_merge_spikes_radius<G>(&mut self, ρ: F, interp: bool, mut accept: G) -> usize
+    where
+        G: FnMut(&'_ Self) -> bool,
+    {
+        // Sort by first coordinate into an indexing array.
+        let mut indices = sort_indices_by(&self.spikes, compare_first_coordinate);
+
+        // Initialise result
+        let mut count = 0;
+        let mut start_scan_2nd = 0;
+
+        // Scan in order
+        if indices.len() == 0 {
+            return count;
+        }
+        for k in 0..indices.len() - 1 {
+            let i = indices[k];
+            let &DeltaMeasure { x: Loc([xi1, xi2]), α: αi } = &self[i];
+
+            if αi == 0.0 {
+                // Nothin to be done if the weight is already zero
+                continue;
+            }
+
+            let mut closest = None;
+
+            // Scan for second spike. We start from `start_scan_2nd + 1` with `start_scan_2nd`
+            // the smallest invalid merging index on the previous loop iteration, because a
+            // the _closest_ mergeable spike might have index less than `k` in `indices`, and a
+            // merge with it might have not been attempted with this spike if a different closer
+            // spike was discovered based on the second coordinate.
+            'scan_2nd: for l in (start_scan_2nd + 1)..indices.len() {
+                if l == k {
+                    // Do not attempt to merge a spike with itself
+                    continue;
+                }
+                let j = indices[l];
+                let &DeltaMeasure { x: Loc([xj1, xj2]), α: αj } = &self[j];
+
+                if xj1 < xi1 - ρ {
+                    // Spike `j = indices[l]` has too low first coordinate. Update starting index
+                    // for next iteration, and continue scanning.
+                    start_scan_2nd = l;
+                    continue 'scan_2nd;
+                } else if xj1 > xi1 + ρ {
+                    // Break out: spike `j = indices[l]` has already too high first coordinate, no
+                    // more close enough spikes can be found due to the sorting of `indices`.
+                    break 'scan_2nd;
+                }
+
+                // If also second coordinate is close enough, attempt merging if closer than
+                // previously discovered mergeable spikes.
+                let d2 = (xi2 - xj2).abs();
+                if αj != 0.0 && d2 <= ρ {
+                    let r1 = xi1 - xj1;
+                    let d = (d2 * d2 + r1 * r1).sqrt();
+                    match closest {
+                        None => closest = Some((l, j, d)),
+                        Some((_, _, r)) if r > d => closest = Some((l, j, d)),
+                        _ => {}
+                    }
+                }
+            }
+
+            // Attempt merging closest close-enough spike
+            if let Some((l, j, _)) = closest {
+                if let Some(n) = self.attempt_merge(i, j, interp, &mut accept) {
+                    // If merge was succesfull, make new spike candidate for merging.
+                    indices[l] = n;
+                    count += 1;
+                    let compare = |i, j| compare_first_coordinate(&self.spikes[i], &self.spikes[j]);
+                    // Re-sort relevant range of indices
+                    if l < k {
+                        indices[l..k].sort_by(|&i, &j| compare(i, j));
+                    } else {
+                        indices[k + 1..=l].sort_by(|&i, &j| compare(i, j));
+                    }
+                }
+            }
+        }
+
+        count
+    }
+}

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/python.rs	Fri Nov 28 12:48:17 2025 -0500
@@ -0,0 +1,142 @@
+/*!
+Python wrapper to measures. Only enabled with crate feature `pyo3`.
+
+These should really be in the `pointsource_pde` crate, but Rust doesn't allow that.
+*/
+
+use super::{DeltaMeasure, DiscreteMeasure};
+use alg_tools::loc::Loc;
+use numpy::{Ix1, PyArray};
+use pyo3::prelude::*;
+use pyo3::types::{PyList, PyModule, PyModuleMethods};
+use pyo3::{pymodule, Borrowed, Bound, PyResult};
+
+macro_rules! create_interface {
+    ($name: ident, $iter:ident, $N: literal, $F:ident) => {
+        #[allow(non_camel_case_types)]
+        #[pyclass(module = "pointsource_algs")]
+        /// Wrapper to [`DiscreteMeasure<Loc<$N, $F>, $F>,`]
+        ///
+        /// This is mainly needed because pyo3 does not support generics, not even instantiating
+        /// them to specific values.
+        pub struct $name {
+            inner: DiscreteMeasure<Loc<$N, $F>, $F>,
+        }
+
+        impl $name {
+            pub fn wrap(inner: DiscreteMeasure<Loc<$N, $F>, $F>) -> Self {
+                Self { inner }
+            }
+        }
+
+        #[pymethods]
+        impl $name {
+            #[new]
+            pub fn new(contents: &Bound<'_, PyList>) -> PyResult<Self> {
+                // let vec: Vec<(PyArray<$F, Ix1>, $F)> = contents.extract()?;
+                // Ok(Self { inner: vec.into() })
+                let mut res = DiscreteMeasure::new();
+                for v in contents.iter() {
+                    let (x, α): ([$F; $N], $F) = v.extract()?;
+                    res.push(DeltaMeasure { x: x.into(), α });
+                }
+                Ok(Self { inner: res })
+            }
+
+            /*pub fn from(vec: Vec<([$F; $N], $F)>) -> Self {
+                Self {
+                    inner: DiscreteMeasure::from_iter(vec.into_iter()),
+                }
+            }*/
+
+            fn __iter__(slf: Bound<'_, Self>) -> PyResult<Py<$iter>> {
+                Py::new(slf.py(), $iter { measure_ref: slf.unbind(), next: 0 })
+            }
+        }
+
+        #[allow(non_camel_case_types)]
+        #[pyclass(module = "pointsource_algs")]
+        /// Python-side iterator for [`DiscreteMeasure<Loc<$N, $F>, $F>,`]
+        /// Returns tuples (weight, coords)
+        pub struct $iter {
+            measure_ref: Py<$name>,
+            next: usize,
+        }
+
+        #[pymethods]
+        impl $iter {
+            fn __iter__(slf: PyRef<'_, Self>) -> PyRef<'_, Self> {
+                slf
+            }
+
+            // We cannot use lifetimes in types exported to Python, so have to use a
+            // very primitive iterator.
+            fn __next__(
+                mut slf: PyRefMut<'_, Self>,
+            ) -> PyResult<Option<(Bound<'_, PyArray<$F, Ix1>>, $F)>> {
+                let py = slf.py();
+                let meas_: PyRef<'_, $name> = slf.measure_ref.extract(py)?;
+                let meas = &(meas_.inner);
+                let next = &mut slf.next;
+                Ok((*next < meas.len()).then(|| {
+                    let δ = meas[*next];
+                    *next += 1;
+                    (PyArray::from_slice(py, &δ.x.0), δ.α)
+                }))
+            }
+        }
+
+        // Direct access without passing through the wrappers
+        impl<'a, 'py> FromPyObject<'a, 'py> for DiscreteMeasure<Loc<$N, $F>, $F> {
+            type Error = PyErr;
+            fn extract(obj: Borrowed<'a, 'py, PyAny>) -> PyResult<Self> {
+                let wrapper: PyRef<'_, $name> = obj.extract()?;
+                Ok(wrapper.inner.clone())
+            }
+        }
+
+        impl<'a, 'py> IntoPyObject<'py> for &'a mut DiscreteMeasure<Loc<$N, $F>, $F> {
+            type Target = $name;
+            type Error = PyErr;
+            type Output = Bound<'py, $name>;
+
+            fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
+                Bound::new(py, $name { inner: self.clone() })
+            }
+        }
+
+        impl<'py> IntoPyObject<'py> for DiscreteMeasure<Loc<$N, $F>, $F> {
+            type Target = $name;
+            type Error = PyErr;
+            type Output = Bound<'py, $name>;
+
+            fn into_pyobject(self, py: Python<'py>) -> Result<Self::Output, Self::Error> {
+                Bound::new(py, $name { inner: self })
+            }
+        }
+    };
+}
+
+create_interface!(DiscreteMeasure_1_f64, DiscreteMeasureIter_1_f64, 1, f64);
+create_interface!(DiscreteMeasure_2_f64, DiscreteMeasureIter_2_f64, 2, f64);
+create_interface!(DiscreteMeasure_3_f64, DiscreteMeasureIter_3_f64, 3, f64);
+
+/// Populates the Python module.
+///
+/// Needs to be called with [`pyo3::append_to_inittab`]:
+/// ```
+/// append_to_inittab!(pymod_pointsource_algs);
+/// ```
+/// before initialising the intepreter.
+#[pymodule]
+#[pyo3(name = "measures")]
+pub fn pymod(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    //m.add_class::<crate::run::DefaultAlgorithm>()?;
+    m.add_class::<crate::python::DiscreteMeasure_1_f64>()?;
+    m.add_class::<crate::python::DiscreteMeasure_2_f64>()?;
+    m.add_class::<crate::python::DiscreteMeasure_3_f64>()?;
+    m.add_class::<crate::python::DiscreteMeasureIter_1_f64>()?;
+    m.add_class::<crate::python::DiscreteMeasureIter_2_f64>()?;
+    m.add_class::<crate::python::DiscreteMeasureIter_3_f64>()?;
+    Ok(())
+}