diff -r 9738b51d90d7 -r 4f468d35fa29 src/measures/discrete.rs
--- a/src/measures/discrete.rs	Sun Apr 27 15:03:51 2025 -0500
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,1014 +0,0 @@
-//! This module implementes discrete measures.
-
-use std::ops::{
-    Div,Mul,DivAssign,MulAssign,Neg,
-    Add,Sub,AddAssign,SubAssign,
-    Index,IndexMut,
-};
-use std::iter::Sum;
-use serde::ser::{Serializer, Serialize, SerializeSeq};
-use nalgebra::DVector;
-
-use alg_tools::norms::Norm;
-use alg_tools::tabledump::TableDump;
-use alg_tools::linops::{Mapping, Linear};
-use alg_tools::iter::{MapF,Mappable};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::collection::Collection;
-use alg_tools::instance::{Instance, Decomposition, MyCow, EitherDecomp, Space};
-
-use crate::types::*;
-use super::base::*;
-use super::delta::*;
-
-/// 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<F, const N : usize> = DiscreteMeasure<Loc<F, N>, 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<F, N>,F>
-where DeltaMeasure<Loc<F, N>, F> : Serialize + 'a {
-    type Iter = std::slice::Iter<'a, DeltaMeasure<Loc<F, N>, 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<F, N>, 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<F, N>, 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<F, Radon> for DiscreteMeasure<Domain, F>
-where DeltaMeasure<Domain, F> : Norm<F, Radon> {
-    #[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 + Space,
-    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 + Space,
-    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 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 decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Owned(self.spikes)
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        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
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for &'a DiscreteMeasure<Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Borrowed(self.spikes.as_slice())
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        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()
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for EitherSlice<'a, Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        self
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        match self {
-            EitherDecomp::Owned(v) => v.as_slice(),
-            EitherDecomp::Borrowed(s) => s,
-        }
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        match self {
-            EitherDecomp::Owned(v) => v.into(),
-            EitherDecomp::Borrowed(s) => s.into(),
-        }
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for &'a EitherSlice<'a, Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        match self {
-            EitherDecomp::Owned(v) => EitherDecomp::Borrowed(v.as_slice()),
-            EitherDecomp::Borrowed(s) => EitherDecomp::Borrowed(s),
-        }
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        match self {
-            EitherDecomp::Owned(v) => v.as_slice(),
-            EitherDecomp::Borrowed(s) => s,
-        }
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        match self {
-            EitherDecomp::Owned(v) => v.as_slice(),
-            EitherDecomp::Borrowed(s) => s
-        }.into()
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for SpikeSlice<'a, Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Borrowed(self)
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        self
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        self.into()
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for &'a SpikeSlice<'a, Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Borrowed(*self)
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        *self
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        (*self).into()
-    }
-}
-
-impl<F : Num, Domain : Clone > Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for DeltaMeasure<Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Owned(vec![self])
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        std::slice::from_ref(self)
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        self.into()
-    }
-}
-
-impl<'a, F : Num, Domain : Clone> Instance<DiscreteMeasure<Domain, F>, MeasureDecomp>
-for &'a DeltaMeasure<Domain, F>
-{
-    fn decompose<'b>(self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Decomposition<'b>
-    where Self : 'b, DiscreteMeasure<Domain, F> : 'b {
-        EitherDecomp::Borrowed(std::slice::from_ref(self))
-    }
-  
-    fn ref_instance(&self)
-        -> <MeasureDecomp as Decomposition<DiscreteMeasure<Domain, F>>>::Reference<'_>
-    {
-        std::slice::from_ref(*self)
-    }
-
-    fn own(self) -> DiscreteMeasure<Domain, F> {
-        self.into()
-    }
-}