Fri, 09 Dec 2022 14:10:48 +0200
Added command line option for (power) tolerance
//! 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::{Apply, Linear}; use alg_tools::iter::{MapF,Mappable}; use alg_tools::nalgebra_support::ToNalgebraRealField; 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>>, } /// 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() } /// 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(α)); } // /// 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.spikes.retain(|δ| δ.α != F::ZERO); } } 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 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) { let mut μ2_get = 0; let mut μ2_insert = 0; self.spikes.drain_filter(|&mut DeltaMeasure{ α : ref mut α_ref, ref x }| { // Get weight of spike in μ2, zero if out of bounds. let β = μ2.spikes.get(μ2_get).map_or(F::ZERO, DeltaMeasure::get_mass); μ2_get += 1; if *α_ref == F::ZERO && β == F::ZERO { // Prune true } else { // Save self weight let α = *α_ref; // Modify self *α_ref = θ * α - ζ * β; // Make copy of old self weight in μ2 let δ = DeltaMeasure{ α, x : x.clone() }; match μ2.spikes.get_mut(μ2_insert) { Some(replace) => { *replace = δ; }, None => { debug_assert_eq!(μ2.len(), μ2_insert); μ2.spikes.push(δ); }, } μ2_insert += 1; // Keep false } }); // Truncate μ2 to same length as self. μ2.spikes.truncate(μ2_insert); debug_assert_eq!(μ2.len(), self.len()); } } 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(α))); } } 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, 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, 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, Y : Sum + Mul<F, Output=Y>> Apply<G> for DiscreteMeasure<Domain, F> where G: for<'a> Apply<&'a Domain, Output = Y> { type Output = Y; #[inline] fn apply(&self, g : G) -> Y { self.spikes.iter().map(|m| g.apply(&m.x) * m.α).sum() } } impl<Domain, G, F : Num, Y : Sum + Mul<F, Output=Y>> Linear<G> for DiscreteMeasure<Domain, F> where G : for<'a> Apply<&'a Domain, Output = Y> { type Codomain = Y; } /// 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>>; /// 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);