pointsource_algs: comparison src/measures/merging.rs

-:aec67cdd6b14
+:efa60bc4f743
 /// Spike merging heuristic selection
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[allow(dead_code)]
 pub enum SpikeMergingMethod<F> {
-/// Try to merge spikes within a given radius of eachother
+/// Try to merge spikes within a given radius of each other, averaging the location
 HeuristicRadius(F),
+/// Try to merge spikes within a given radius of each other, attempting original locations
+HeuristicRadiusNoInterp(F),
 /// No merging
 None,
 }
 // impl<F : Float> SpikeMergingMethod<F> {
 impl<F : ClapFloat> std::fmt::Display for SpikeMergingMethod<F> {
 fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> Result<(), std::fmt::Error> {
 match self {
 Self::None => write!(f, "none"),
-Self::HeuristicRadius(r) => std::fmt::Display::fmt(r, f),
+Self::HeuristicRadius(r) => write!(f, "i:{}", r),
+Self::HeuristicRadiusNoInterp(r) => write!(f, "n:{}", r),
 }
 }
 }
 impl<F : ClapFloat> std::str::FromStr for SpikeMergingMethod<F> {
 fn from_str(s: &str) -> Result<Self, Self::Err> {
 if s == "none" {
 Ok(Self::None)
 } else {
-Ok(Self::HeuristicRadius(F::from_str(s)?))
+let mut subs = s.split(':');
+match subs.next() {
+None =>  Ok(Self::HeuristicRadius(F::from_str(s)?)),
+Some(t) if t == "n" => match subs.next() {
+None => Err(core::num::dec2flt::pfe_invalid()),
+Some(v) => Ok(Self::HeuristicRadiusNoInterp(F::from_str(v)?))
+},
+Some(t) if t == "i" => match subs.next() {
+None => Err(core::num::dec2flt::pfe_invalid()),
+Some(v) => Ok(Self::HeuristicRadius(F::from_str(v)?))
+},
+Some(v) => Ok(Self::HeuristicRadius(F::from_str(v)?))
+}
 }
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 /// 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 weight of existing
 /// spikes is set to zero, and a new one inserted at the end of the spike vector.
-fn merge_spikes<G, V>(&mut self, method : SpikeMergingMethod<F>, accept : G) -> Option<V>
+fn merge_spikes<G>(&mut self, method : SpikeMergingMethod<F>, accept : G) -> usize
-where G : Fn(&'_ Self) -> Option<V> {
+where G : FnMut(&'_ Self) -> bool {
 match method {
-SpikeMergingMethod::HeuristicRadius(ρ) => self.do_merge_spikes_radius(ρ, accept),
+SpikeMergingMethod::HeuristicRadius(ρ) =>
-SpikeMergingMethod::None => None,
+self.do_merge_spikes_radius(ρ, true, accept),
+SpikeMergingMethod::HeuristicRadiusNoInterp(ρ) =>
+self.do_merge_spikes_radius(ρ, false, accept),
+SpikeMergingMethod::None => 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.
+/// `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
+) -> (V, usize)
 where G : Fn(&'_ Self) -> V,
 H : Fn(&'_ V) -> O,
 O : PartialOrd {
-let initial_res = value(self);
+let mut res = value(self);
-let initial_fitness = fitness(&initial_res);
+let initial_fitness = fitness(&res);
-self.merge_spikes(method, |μ| {
+let count = self.merge_spikes(method, |μ| {
-let res = value(μ);
+res = value(μ);
-(fitness(&res) <= initial_fitness).then_some(res)
+fitness(&res) <= initial_fitness
-}).unwrap_or(initial_res)
+});
+(res, count)
 }
 /// Attempt to merge spikes that are within radius $ρ$ of each other (unspecified norm).
 ///
 /// This method implements [`SpikeMerging::merge_spikes`] for
 /// [`SpikeMergingMethod::HeuristicRadius`]. The closure `accept` and the return value are
 /// as for that method.
-fn do_merge_spikes_radius<G, V>(&mut self, ρ : F, accept : G) -> Option<V>
+fn do_merge_spikes_radius<G>(&mut self, ρ : F, interp : bool, accept : G) -> usize
-where G : Fn(&'_ Self) -> Option<V>;
+where G : FnMut(&'_ Self) -> bool;
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float, const N : usize>  DiscreteMeasure<Loc<F, N>, 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.
-fn attempt_merge<G, V>(
+///
+/// Returns the index of `self.spikes` storing the new spike.
+fn attempt_merge<G>(
 &mut self,
-res : &mut Option<V>,
 i : usize,
 j : usize,
-accept : &G
+interp : bool,
-) -> bool
+accept : &mut G
-where G : Fn(&'_ Self) -> Option<V> {
+) -> Option<usize>
+where G : FnMut(&'_ Self) -> bool {
 let &DeltaMeasure{ x : xi, α : αi } = &self.spikes[i];
 let &DeltaMeasure{ x : xj, α : αj } = &self.spikes[j];
-// Merge inplace
+if interp {
-self.spikes[i].α = 0.0;
+// Merge inplace
-self.spikes[j].α = 0.0;
+self.spikes[i].α = 0.0;
-//self.spikes.push(DeltaMeasure{ α : αi + αj, x : (xi + xj)/2.0 });
+self.spikes[j].α = 0.0;
-self.spikes.push(DeltaMeasure{ α : αi + αj, x : (xi * αi + xj * αj) / (αi + αj) });
+let αia = αi.abs();
-match accept(self) {
+let αja = αj.abs();
-some@Some(..) => {
+self.spikes.push(DeltaMeasure{ α : αi + αj, x : (xi * αia + xj * αja) / (αia + αja) });
-// Merge accepted, update our return value
+if accept(self) {
-*res = some;
+Some(self.spikes.len()-1)
-// On next iteration process the newly merged spike.
+} else {
-//indices[k+1] = self.spikes.len() - 1;
-true
-},
-None => {
 // Merge not accepted, restore modification
 self.spikes[i].α = αi;
 self.spikes[j].α = αj;
 self.spikes.pop();
-false
+None
 }
-}
+} else {
-}
+// Attempt merge inplace, first combination
+self.spikes[i].α = αi + αj;
-/*
+self.spikes[j].α = 0.0;
-/// Attempts to merge spikes with indices i and j, acceptance through a delta.
+if accept(self) {
-fn attempt_merge_change<G, V>(
+// Merge accepted
-&mut self,
+Some(i)
-res : &mut Option<V>,
+} else {
-i : usize,
+// Attempt merge inplace, second combination
-j : usize,
-accept_change : &G
-) -> bool
-where G : Fn(&'_ Self) -> Option<V> {
-let &DeltaMeasure{ x : xi, α : αi } = &self.spikes[i];
-let &DeltaMeasure{ x : xj, α : αj } = &self.spikes[j];
-let δ = DeltaMeasure{ α : αi + αj, x : (xi + xj)/2.0 };
-let λ = [-self.spikes[i], -self.spikes[j], δ.clone()].into();
-match accept_change(&λ) {
-some@Some(..) => {
-// Merge accepted, update our return value
-*res = some;
 self.spikes[i].α = 0.0;
-self.spikes[j].α = 0.0;
+self.spikes[j].α = αi + αj;
-self.spikes.push(δ);
+if accept(self) {
-true
+// Merge accepted
-},
+Some(j)
-None => {
+} else {
-false
+// 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`.
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> SpikeMerging<F> for DiscreteMeasure<Loc<F, 1>, F> {
-fn do_merge_spikes_radius<G, V>(
+fn do_merge_spikes_radius<G>(
 &mut self,
 ρ : F,
-accept : G
+interp : bool,
-) -> Option<V>
+mut accept : G
-where G : Fn(&'_ Self) -> Option<V> {
+) -> 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 res = None;
+let mut count = 0;
 // Scan consecutive pairs and merge if close enough and accepted by `accept`.
 if indices.len() == 0 {
-return res
+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 self.attempt_merge(&mut res, i, j, &accept) {
+if let Some(l) = self.attempt_merge(i, j, interp, &mut accept) {
-indices[k+1] = self.spikes.len() - 1;
+// 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
+}
-res
+}
+}
+count
 }
 }
 /// Orders `δ1` and `δ1` according to the first coordinate.
 fn compare_first_coordinate<F : Float>(
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> SpikeMerging<F> for DiscreteMeasure<Loc<F, 2>, F> {
-fn do_merge_spikes_radius<G, V>(&mut self, ρ : F, accept : G) -> Option<V>
+fn do_merge_spikes_radius<G>(&mut self, ρ : F, interp : bool, mut accept : G) -> usize
-where G : Fn(&'_ Self) -> Option<V> {
+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 res = None;
+let mut count = 0;
 let mut start_scan_2nd = 0;
 // Scan in order
 if indices.len() == 0 {
-return res
+return count
 }
 for k in 0..indices.len()-1 {
 let i = indices[k];
 let &DeltaMeasure{ x : Loc([xi1, xi2]), α : αi } = &self[i];
 }
 }
 // Attempt merging closest close-enough spike
 if let Some((l, j, _)) = closest {
-if self.attempt_merge(&mut res, i, j, &accept) {
+if let Some(n) = self.attempt_merge(i, j, interp, &mut accept) {
 // If merge was succesfull, make new spike candidate for merging.
-indices[l] = self.spikes.len() - 1;
+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));
 }
 }
 }
 }
-res
+count
 }
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/measures/merging.rs

src/measures/merging.rs