pointsource_algs: comparison src/sliding

-:aec67cdd6b14
+:efa60bc4f743
 use numeric_literals::replace_float_literals;
 use serde::{Serialize, Deserialize};
 //use colored::Colorize;
 //use nalgebra::{DVector, DMatrix};
 use itertools::izip;
-use std::iter::{Map, Flatten};
+use std::iter::Iterator;
 use alg_tools::iterate::{
 AlgIteratorFactory,
 AlgIteratorState
 };
-use alg_tools::euclidean::{
+use alg_tools::euclidean::Euclidean;
-Euclidean,
-Dot
-};
 use alg_tools::sets::Cube;
 use alg_tools::loc::Loc;
 use alg_tools::mapping::{Apply, Differentiable};
+use alg_tools::norms::{Norm, L2};
 use alg_tools::bisection_tree::{
 BTFN,
 PreBTFN,
 Bounds,
 BTNodeLookup,
 };
 use alg_tools::mapping::RealMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use crate::types::*;
-use crate::measures::{
+use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon};
-DiscreteMeasure,
-DeltaMeasure,
-};
 use crate::measures::merging::{
 //SpikeMergingMethod,
 SpikeMerging,
 };
 use crate::forward_model::ForwardModel;
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct SlidingFBConfig<F : Float> {
 /// Step length scaling
 pub τ0 : F,
-/// Transport smoothness assumption
+/// Transport step length $θ$ normalised to $(0, 1)$.
-pub ℓ0 : F,
+pub θ0 : F,
-/// Inverse of the scaling factor $θ$ of the 2-norm-squared transport cost.
+/// Maximum transport mass scaling.
-/// This means that $τθ$ is the step length for the transport step.
+// /// The maximum transported mass is this factor times $\norm{b}^2/(2α)$.
-pub inverse_transport_scaling : F,
+// pub max_transport_scale : F,
-/// Factor for deciding transport reduction based on smoothness assumption violation
+/// Transport tolerance wrt. ω
-pub minimum_goodness_factor : F,
+pub transport_tolerance_ω : F,
-/// Maximum rays to retain in transports from each source.
+/// Transport tolerance wrt. ∇v
-pub maximum_rays : usize,
+pub transport_tolerance_dv : F,
 /// Generic parameters
 pub insertion : FBGenericConfig<F>,
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> Default for SlidingFBConfig<F> {
 fn default() -> Self {
 SlidingFBConfig {
 τ0 : 0.99,
-ℓ0 : 1.5,
+θ0 : 0.99,
-inverse_transport_scaling : 1.0,
+//max_transport_scale : 10.0,
-minimum_goodness_factor : 1.0, // TODO: totally arbitrary choice,
+transport_tolerance_ω : 1.0, // TODO: no idea what this should be
-// should be scaled by problem data?
+transport_tolerance_dv : 1.0, // TODO: no idea what this should be
-maximum_rays : 10,
 insertion : Default::default()
 }
 }
 }
-/// A transport ray (including various additional computational information).
+/// Scale each |γ|_i ≠ 0 by q_i=q̄/g(γ_i)
-#[derive(Clone, Debug)]
-pub struct Ray<Domain, F : Num> {
-/// The destination of the ray, and the mass. The source is indicated in a [`RaySet`].
-δ : DeltaMeasure<Domain, F>,
-/// Goodness of the data term for the aray: $v(z)-v(y)-⟨∇v(x), z-y⟩ + ℓ‖z-y‖^2$.
-goodness : F,
-/// Goodness of the regularisation term for the ray: $w(z)-w(y)$.
-/// Initially zero until $w$ can be constructed.
-reg_goodness : F,
-/// Indicates that this ray also forms a component in γ^{k+1} with the mass `to_return`.
-to_return : F,
-}
-/// A set of transport rays with the same source point.
-#[derive(Clone, Debug)]
-pub struct RaySet<Domain, F : Num> {
-/// Source of every ray in thset
-source : Domain,
-/// Mass of the diagonal ray, with destination the same as the source.
-diagonal: F,
-/// Goodness of the data term for the diagonal ray with $z=x$:
-/// $v(x)-v(y)-⟨∇v(x), x-y⟩ + ℓ‖x-y‖^2$.
-diagonal_goodness : F,
-/// Goodness of the data term for the diagonal ray with $z=x$: $w(x)-w(y)$.
-diagonal_reg_goodness : F,
-/// The non-diagonal rays.
-rays : Vec<Ray<Domain, F>>,
-}
 #[replace_float_literals(F::cast_from(literal))]
-impl<Domain, F : Float> RaySet<Domain, F> {
+fn scale_down<'a, I, F, G, const N : usize>(
-fn non_diagonal_mass(&self) -> F {
+iter : I,
-self.rays
+q̄ : F,
-.iter()
+mut g : G
-.map(|Ray{ δ : DeltaMeasure{ α, .. }, .. }| *α)
+) where F : Float,
-.sum()
+I : Iterator<Item = &'a mut DeltaMeasure<Loc<F,N>, F>>,
-}
+G : FnMut(&DeltaMeasure<Loc<F,N>, F>) -> F {
+iter.for_each(|δ| {
-fn total_mass(&self) -> F {
+if δ.α != 0.0 {
-self.non_diagonal_mass() + self.diagonal
+let b = g(δ);
-}
+if b * δ.α > 0.0 {
+δ.α *= q̄/b;
-fn targets<'a>(&'a self)
+}
--> Map<
+}
-std::slice::Iter<'a, Ray<Domain, F>>,
+});
-fn(&'a Ray<Domain, F>) -> &'a DeltaMeasure<Domain, F>
-> {
-fn get_δ<'b, Domain, F : Float>(Ray{ δ, .. }: &'b Ray<Domain, F>)
--> &'b DeltaMeasure<Domain, F> {
-δ
-}
-self.rays
-.iter()
-.map(get_δ)
-}
-// fn non_diagonal_goodness(&self) -> F {
-//     self.rays
-//         .iter()
-//         .map(|&Ray{ δ : DeltaMeasure{ α, .. }, goodness, reg_goodness, .. }| {
-//             α * (goodness + reg_goodness)
-//         })
-//         .sum()
-// }
-// fn total_goodness(&self) -> F {
-//     self.non_diagonal_goodness() + (self.diagonal_goodness + self.diagonal_reg_goodness)
-// }
-fn non_diagonal_badness(&self) -> F {
-self.rays
-.iter()
-.map(|&Ray{ δ : DeltaMeasure{ α, .. }, goodness, reg_goodness, .. }| {
-0.0.max(- α * (goodness + reg_goodness))
-})
-.sum()
-}
-fn total_badness(&self) -> F {
-self.non_diagonal_badness()
-+ 0.0.max(- self.diagonal * (self.diagonal_goodness + self.diagonal_reg_goodness))
-}
-fn total_return(&self) -> F {
-self.rays
-.iter()
-.map(|&Ray{ to_return, .. }| to_return)
-.sum()
-}
-}
-#[replace_float_literals(F::cast_from(literal))]
-impl<Domain : Clone, F : Num> RaySet<Domain, F> {
-fn return_targets<'a>(&'a self)
--> Flatten<Map<
-std::slice::Iter<'a, Ray<Domain, F>>,
-fn(&'a Ray<Domain, F>) -> Option<DeltaMeasure<Domain, F>>
->> {
-fn get_return<'b, Domain : Clone, F : Num>(ray: &'b Ray<Domain, F>)
--> Option<DeltaMeasure<Domain, F>> {
-(ray.to_return != 0.0).then_some(
-DeltaMeasure{x : ray.δ.x.clone(), α : ray.to_return}
-)
-}
-let tmp : Map<
-std::slice::Iter<'a, Ray<Domain, F>>,
-fn(&'a Ray<Domain, F>) -> Option<DeltaMeasure<Domain, F>>
-> = self.rays
-.iter()
-.map(get_return);
-tmp.flatten()
-}
 }
 /// Iteratively solve the pointsource localisation problem using sliding forward-backward
 /// splitting
 ///
 /// The parametrisatio is as for [`pointsource_fb_reg`].
 /// Inertia is currently not supported.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_sliding_fb_reg<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Reg, const N : usize>(
+pub fn pointsource_sliding_fb_reg<'a, F, I, A, GA, 𝒟, BTA, BT𝒟, G𝒟, S, K, Reg, const N : usize>(
 opA : &'a A,
 b : &A::Observable,
 reg : Reg,
 op𝒟 : &'a 𝒟,
 sfbconfig : &SlidingFBConfig<F>,
 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
 A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
 + Lipschitz<&'a 𝒟, FloatType=F> + TransportLipschitz<L2Squared, FloatType=F>,
 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
-𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
+𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>,
-𝒟::Codomain : RealMapping<F, N>,
+Codomain = BTFN<F, G𝒟, BT𝒟, N>>,
+BT𝒟 : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>
 + Differentiable<Loc<F, N>, Output=Loc<F,N>>,
 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
 //+ Differentiable<Loc<F, N>, Output=Loc<F,N>>,
 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
 PlotLookup : Plotting<N>,
 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
 Reg : SlidingRegTerm<F, N> {
 assert!(sfbconfig.τ0 > 0.0 &&
-sfbconfig.inverse_transport_scaling > 0.0 &&
+sfbconfig.θ0 > 0.0);
-sfbconfig.ℓ0 > 0.0);
 // Set up parameters
 let config = &sfbconfig.insertion;
 let op𝒟norm = op𝒟.opnorm_bound();
-let θ = sfbconfig.inverse_transport_scaling;
+//let max_transport = sfbconfig.max_transport_scale
-let τ = sfbconfig.τ0/opA.lipschitz_factor(&op𝒟).unwrap()
+//                    * reg.radon_norm_bound(b.norm2_squared() / 2.0);
-.max(opA.transport_lipschitz_factor(L2Squared) * θ);
+//let tlip = opA.transport_lipschitz_factor(L2Squared) * max_transport;
-let ℓ = sfbconfig.ℓ0; // TODO: v scaling?
+//let ℓ = 0.0;
+let θ = sfbconfig.θ0; // (ℓ + tlip);
+let τ = sfbconfig.τ0/opA.lipschitz_factor(&op𝒟).unwrap();
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 // by τ compared to the conditional gradient approach.
 let tolerance = config.tolerance * τ * reg.tolerance_scaling();
 let mut ε = tolerance.initial();
 // Initialise iterates
-let mut μ : DiscreteMeasure<Loc<F, N>, F> = DiscreteMeasure::new();
+let mut μ = DiscreteMeasure::new();
-let mut μ_transported_base = DiscreteMeasure::new();
+let mut γ1 = DiscreteMeasure::new();
-let mut γ_hat : Vec<RaySet<Loc<F, N>, F>> = Vec::new();   // γ̂_k and extra info
 let mut residual = -b;
 let mut stats = IterInfo::new();
 // Run the algorithm
 iterator.iterate(|state| {
 // Calculate smooth part of surrogate model.
 // Using `std::mem::replace` here is not ideal, and expects that `empty_observable`
 // has no significant overhead. For some reosn Rust doesn't allow us simply moving
 // the residual and replacing it below before the end of this closure.
-residual *= -τ;
 let r = std::mem::replace(&mut residual, opA.empty_observable());
-let minus_τv = opA.preadjoint().apply(r);
+let v = opA.preadjoint().apply(r);
-// Save current base point and shift μ to new positions.
+// Save current base point and shift μ to new positions. Idea is that
-let μ_base = μ.clone();
+//  μ_base(_masses) = μ^k (vector of masses)
-for δ in μ.iter_spikes_mut() {
+//  μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1}
-δ.x += minus_τv.differential(&δ.x) * θ;
+//  γ1 = π_♯^1γ^{k+1}
-}
+//  μ = μ^{k+1}
-let mut μ_transported = μ.clone();
+let μ_base_masses : Vec<F> = μ.iter_masses().collect();
+let mut μ_base_minus_γ0 = μ.clone(); // Weights will be set in the loop below.
-assert_eq!(μ.len(), γ_hat.len());
+// Construct μ^{k+1} and π_♯^1γ^{k+1} initial candidates
+let mut sum_norm_dv_times_γinit = 0.0;
-// Calculate the goodness λ formed from γ_hat (≈ γ̂_k) and γ^{k+1}, where the latter
+let mut sum_abs_γinit = 0.0;
-// transports points x from μ_base to points y in μ as shifted above, or “returns”
+//let mut sum_norm_dv = 0.0;
-// them “home” to z given by the rays in γ_hat. Returning is necessary if the rays
+let γ_prev_len = γ1.len();
-// are not “good” for the smoothness assumptions, or if γ_hat has more mass than
+assert!(μ.len() >= γ_prev_len);
-// μ_base.
+γ1.extend(μ[γ_prev_len..].iter().cloned());
-let mut total_goodness = 0.0;     // data term goodness
+for (δ, ρ) in izip!(μ.iter_spikes_mut(), γ1.iter_spikes_mut()) {
-let mut total_reg_goodness = 0.0; // regulariser goodness
+let d_v_x = v.differential(&δ.x);
-let minimum_goodness = - ε * sfbconfig.minimum_goodness_factor;
+// If old transport has opposing sign, the new transport will be none.
+ρ.α = if (ρ.α > 0.0 && δ.α < 0.0) || (ρ.α < 0.0 && δ.α > 0.0) {
-for (δ, r) in izip!(μ.iter_spikes(), γ_hat.iter_mut()) {
+0.0
-// Calculate data term goodness for all rays.
+} else {
-let &DeltaMeasure{ x : ref y, α : δ_mass } = δ;
+δ.α
-let x = &r.source;
-let mvy = minus_τv.apply(y);
-let mdvx = minus_τv.differential(x);
-let mut r_total_mass = 0.0; // Total mass of all rays with source r.source.
-let mut bad_mass = 0.0;
-let mut calc_goodness = |goodness : &mut F, reg_goodness : &mut F, α, z : &Loc<F, N>| {
-*reg_goodness = 0.0; // Initial guess
-*goodness = mvy - minus_τv.apply(z) + mdvx.dot(&(z-y))
-+ ℓ * z.dist2_squared(&y);
-total_goodness += *goodness * α;
-r_total_mass += α; // TODO: should this include to_return from staging? (Probably not)
-if *goodness < 0.0 {
-bad_mass += α;
-}
 };
-for ray in r.rays.iter_mut() {
+δ.x -= d_v_x * (θ * δ.α.signum()); // This is δ.α.signum() when δ.α ≠ 0.
-calc_goodness(&mut ray.goodness, &mut ray.reg_goodness, ray.δ.α, &ray.δ.x);
+ρ.x = δ.x;
+let nrm = d_v_x.norm(L2);
+let a = ρ.α.abs();
+let v = nrm * a;
+if v > 0.0 {
+sum_norm_dv_times_γinit += v;
+sum_abs_γinit += a;
 }
-calc_goodness(&mut r.diagonal_goodness, &mut r.diagonal_reg_goodness, r.diagonal, x);
+}
-// If the total mass of the ray set is less than that of μ at the same source,
+// A priori transport adaptation based on bounding ∫ ⟨∇v(x), z-y⟩ dλ(x, y, z).
-// a diagonal component needs to be added to be able to (attempt to) transport
+// This is just one option, there are many.
-// all mass of μ. In the opposite case, we need to construct γ_{k+1} to ‘return’
+let t = ε * sfbconfig.transport_tolerance_dv;
-// the the extra mass of γ̂_k to the target z. We return mass from the oldest “bad”
+if sum_norm_dv_times_γinit > t {
-// rays in the set.
+// Scale each |γ|_i by q_i=q̄/‖vx‖_i such that ∑_i |γ|_i q_i ‖vx‖_i = t
-if δ_mass >= r_total_mass {
+// TODO: store the closure values above?
-r.diagonal += δ_mass - r_total_mass;
+scale_down(γ1.iter_spikes_mut(),
-} else {
+t / sum_abs_γinit,
-let mut reduce_transport = r_total_mass - δ_mass;
+|δ| v.differential(&δ.x).norm(L2));
-let mut good_needed = (bad_mass - reduce_transport).max(0.0);
+}
-// NOTE: reg_goodness is zero at this point, so it is not used in this code.
+//println!("|γ| = {}, |μ| = {}", γ1.norm(crate::measures::Radon), μ.norm(crate::measures::Radon));
-let mut reduce_ray = |goodness, to_return : Option<&mut F>, α : &mut F| {
-if reduce_transport > 0.0 {
-let return_amount = if goodness < 0.0 {
-α.min(reduce_transport)
-} else {
-let amount = α.min(good_needed);
-good_needed -= amount;
-amount
-};
-if return_amount > 0.0 {
-reduce_transport -= return_amount;
-// Adjust total goodness by returned amount
-total_goodness -= goodness * return_amount;
-to_return.map(|tr| *tr += return_amount);
-*α -= return_amount;
-*α > 0.0
-} else {
-true
-}
-} else {
-true
-}
-};
-r.rays.retain_mut(|ray| {
-reduce_ray(ray.goodness, Some(&mut ray.to_return), &mut ray.δ.α)
-});
-// A bad diagonal is simply reduced without any 'return'.
-// It was, after all, just added to match μ, but there is no need to match it.
-// It's just a heuristic.
-// TODO: Maybe a bad diagonal should be the first to go.
-reduce_ray(r.diagonal_goodness, None, &mut r.diagonal);
-}
-}
 // Solve finite-dimensional subproblem several times until the dual variable for the
 // regularisation term conforms to the assumptions made for the transport above.
 let (d, within_tolerances) = 'adapt_transport: loop {
-// If transport violates goodness requirements, shift it to ‘return’ mass to z,
+// Update weights for μ_base_minus_γ0 = μ^k - π_♯^0γ^{k+1}
-// forcing y = z. Based on the badness of each ray set (sum of bad rays' goodness),
+for (δ_γ1, δ_μ_base_minus_γ0, &α_μ_base) in izip!(γ1.iter_spikes(),
-// we proportionally distribute the reductions to each ray set, and within each ray
+μ_base_minus_γ0.iter_spikes_mut(),
-// set, prioritise reducing the oldest bad rays' weight.
+μ_base_masses.iter()) {
-let tg = total_goodness + total_reg_goodness;
+δ_μ_base_minus_γ0.set_mass(α_μ_base - δ_γ1.get_mass());
-let adaptation_needed = minimum_goodness - tg;
+}
-if adaptation_needed > 0.0 {
-let total_badness = γ_hat.iter().map(|r| r.total_badness()).sum();
+// Calculate transported_minus_τv = -τA_*(A[μ_transported + μ_transported_base]-b)
+let residual_μ̆ = calculate_residual2(&γ1, &μ_base_minus_γ0, opA, b);
-let mut return_ray = |goodness : F,
+let transported_minus_τv̆ = opA.preadjoint().apply(residual_μ̆ * (-τ));
-reg_goodness : F,
-to_return : Option<&mut F>,
+// Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
-α : &mut F,
+let (d, within_tolerances) = insert_and_reweigh(
-left_to_return : &mut F| {
+&mut μ, &transported_minus_τv̆, &γ1, Some(&μ_base_minus_γ0),
-let g = goodness + reg_goodness;
+op𝒟, op𝒟norm,
-assert!(*α >= 0.0 && *left_to_return >= 0.0);
+τ, ε,
-if *left_to_return > 0.0 && g < 0.0 {
+config,
-let return_amount = (*left_to_return / (-g)).min(*α);
+&reg, state, &mut stats,
-*left_to_return -= (-g) * return_amount;
+);
-total_goodness -= goodness * return_amount;
-total_reg_goodness -= reg_goodness * return_amount;
+// A posteriori transport adaptation based on bounding (1/τ)∫ ω(z) - ω(y) dλ(x, y, z).
-to_return.map(|tr| *tr += return_amount);
+let all_ok = if false { // Basic check
-*α -= return_amount;
+// If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not,
-*α > 0.0
+// then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1
-} else {
+// at that point to zero, and retry.
-true
+let mut all_ok = true;
-}
+for (α_μ, α_γ1) in izip!(μ.iter_masses(), γ1.iter_masses_mut()) {
-};
+if α_μ == 0.0 && *α_γ1 != 0.0 {
+all_ok = false;
-for r in γ_hat.iter_mut() {
+*α_γ1 = 0.0;
-let mut left_to_return = adaptation_needed * r.total_badness() / total_badness;
-if left_to_return > 0.0 {
-for ray in r.rays.iter_mut() {
-return_ray(ray.goodness, ray.reg_goodness,
-Some(&mut ray.to_return), &mut ray.δ.α, &mut left_to_return);
-}
-return_ray(r.diagonal_goodness, r.diagonal_reg_goodness,
-None, &mut r.diagonal, &mut left_to_return);
 }
 }
-}
+all_ok
+} else {
-// Construct μ_k + (π_#^1-π_#^0)γ_{k+1}.
+// TODO: Could maybe optimise, as this is also formed in insert_and_reweigh above.
-// This can be broken down into
+let mut minus_ω = op𝒟.apply(γ1.sub_matching(&μ) + &μ_base_minus_γ0);
-//
-// μ_transported_base = [μ - π_#^0 (γ_shift + γ_return)] + π_#^1 γ_return, and
+// let vpos = γ1.iter_spikes()
-// μ_transported = π_#^1 γ_shift
+//              .filter(|δ| δ.α > 0.0)
-//
+//              .map(|δ| minus_ω.apply(&δ.x))
-// where γ_shift is our “true” γ_{k+1}, and γ_return is the return compoennt.
+//              .reduce(F::max)
-// The former can be constructed from δ.x and δ_new.x for δ in μ_base and δ_new in μ
+//              .and_then(|threshold| {
-// (which has already been shifted), and the mass stored in a γ_hat ray's δ measure
+//                 minus_ω.minimise_below(threshold,
-// The latter can be constructed from γ_hat rays' source and destination with the
+//                                         ε * config.refinement.tolerance_mult,
-// to_return mass.
+//                                         config.refinement.max_steps)
-//
+//                        .map(|(_z, minus_ω_z)| minus_ω_z)
-// Note that μ_transported is constructed to have the same spike locations as μ, but
+//              });
-// to have same length as μ_base. This loop does not iterate over the spikes of μ
-// (and corresponding transports of γ_hat) that have been newly     added in the current
+// let vneg = γ1.iter_spikes()
-// 'adapt_transport loop.
+//              .filter(|δ| δ.α < 0.0)
-for (δ, δ_transported, r) in izip!(μ_base.iter_spikes(),
+//              .map(|δ| minus_ω.apply(&δ.x))
-μ_transported.iter_spikes_mut(),
+//              .reduce(F::min)
-γ_hat.iter()) {
+//              .and_then(|threshold| {
-let &DeltaMeasure{ref x, α} = δ;
+//                 minus_ω.maximise_above(threshold,
-debug_assert_eq!(*x, r.source);
+//                                         ε * config.refinement.tolerance_mult,
-let shifted_mass = r.total_mass();
+//                                         config.refinement.max_steps)
-let ret_mass = r.total_return();
+//                        .map(|(_z, minus_ω_z)| minus_ω_z)
-// μ - π_#^0 (γ_shift + γ_return)
+//              });
-μ_transported_base += DeltaMeasure { x : *x, α : α - shifted_mass - ret_mass };
+let (_, vpos) = minus_ω.minimise(ε * config.refinement.tolerance_mult,
-// π_#^1 γ_return
+config.refinement.max_steps);
-μ_transported_base.extend(r.return_targets());
+let (_, vneg) = minus_ω.maximise(ε * config.refinement.tolerance_mult,
-// π_#^1 γ_shift
+config.refinement.max_steps);
-δ_transported.set_mass(shifted_mass);
-}
+let t = τ * ε * sfbconfig.transport_tolerance_ω;
-// Calculate transported_minus_τv = -τA_*(A[μ_transported + μ_transported_base]-b)
+let val = |δ : &DeltaMeasure<Loc<F, N>, F>| {
-let transported_residual = calculate_residual2(&μ_transported,
+δ.α * (minus_ω.apply(&δ.x) - if δ.α >= 0.0 { vpos } else { vneg })
-&μ_transported_base,
+// match if δ.α >= 0.0 { vpos } else { vneg } {
-opA, b);
+//     None => 0.0,
-let transported_minus_τv = opA.preadjoint()
+//     Some(v) => δ.α * (minus_ω.apply(&δ.x) - v)
-.apply(transported_residual);
+// }
+};
-// Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
+// Calculate positive/bad (rp) values under the integral.
-let (mut d, within_tolerances) = insert_and_reweigh(
+// Also store sum of masses for the positive entries.
-&mut μ, &transported_minus_τv, &μ_transported, Some(&μ_transported_base),
+let (rp, w) = γ1.iter_spikes().fold((0.0, 0.0), |(p, w), δ| {
-op𝒟, op𝒟norm,
+let v = val(δ);
-τ, ε,
+if v <= 0.0 { (p, w) } else { (p + v, w + δ.α.abs()) }
-config, &reg, state, &mut stats
+});
-);
+if rp > t {
-// We have  d = ω0 - τv - 𝒟μ = -𝒟(μ - μ^k) - τv; more precisely
+// TODO: store v above?
-//          d = minus_τv + op𝒟.preapply(μ_diff(μ, μ_transported, config));
+scale_down(γ1.iter_spikes_mut(), t / w, val);
-// We “essentially” assume that the subdifferential w of the regularisation term
+false
-// satisfies w'(y)=0, so for a “goodness” estimate τ[w(y)-w(z)-w'(y)(z-y)]
+} else {
-// that incorporates the assumption, we need to calculate τ[w(z) - w(y)] for
+true
-// some w in the subdifferential of the regularisation term, such that
-// -ε ≤ τw - d ≤ ε. This is done by [`RegTerm::goodness`].
-for r in γ_hat.iter_mut() {
-for ray in r.rays.iter_mut() {
-ray.reg_goodness = reg.goodness(&mut d, &μ, &r.source, &ray.δ.x, τ, ε, config);
-total_reg_goodness += ray.reg_goodness * ray.δ.α;
 }
-}
+};
-// If update of regularisation term goodness didn't invalidate minimum goodness
+if all_ok {
-// requirements, we have found our step. Otherwise we need to keep reducing
-// transport by repeating the loop.
-if total_goodness + total_reg_goodness >= minimum_goodness {
 break 'adapt_transport (d, within_tolerances)
 }
 };
-// Update γ_hat to new location
+stats.untransported_fraction = Some({
-for (δ_new, r) in izip!(μ.iter_spikes(), γ_hat.iter_mut()) {
+assert_eq!(μ_base_masses.len(), γ1.len());
-// Prune rays that only had a return component, as the return component becomes
+let (a, b) = stats.untransported_fraction.unwrap_or((0.0, 0.0));
-// a diagonal in γ̂^{k+1}.
+let source = μ_base_masses.iter().map(|v| v.abs()).sum();
-r.rays.retain(|ray| ray.δ.α != 0.0);
+(a + μ_base_minus_γ0.norm(Radon), b + source)
-// Otherwise zero out the return component, or stage rays for pruning
+});
-// to keep memory and computational demands reasonable.
+stats.transport_error = Some({
-let n_rays = r.rays.len();
+assert_eq!(μ_base_masses.len(), γ1.len());
-for (ray, ir) in izip!(r.rays.iter_mut(), (0..n_rays).rev()) {
+let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0));
-if ir >= sfbconfig.maximum_rays {
+let err = izip!(μ.iter_masses(), γ1.iter_masses()).map(|(v,w)| (v-w).abs()).sum();
-// Only keep sfbconfig.maximum_rays - 1 previous rays, staging others for
+(a + err, b + γ1.norm(Radon))
-// pruning in next step.
+});
-ray.to_return = ray.δ.α;
-ray.δ.α = 0.0;
+// Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the
-} else {
+// latter needs to be pruned when μ is.
-ray.to_return = 0.0;
+// TODO: This could do with a two-vector Vec::retain to avoid copies.
-}
+let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned());
-ray.goodness = 0.0; // TODO: probably not needed
+if μ_new.len() != μ.len() {
-ray.reg_goodness = 0.0;
+let mut μ_iter = μ.iter_spikes();
-}
+γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO);
-// Add a new ray for the currently diagonal component
+μ = μ_new;
-if r.diagonal > 0.0 {
-r.rays.push(Ray{
-δ : DeltaMeasure{x : r.source, α : r.diagonal},
-goodness : 0.0,
-reg_goodness : 0.0,
-to_return : 0.0,
-});
-// TODO: Maybe this does not need to be done here, and is sufficent to to do where
-// the goodness is calculated.
-r.diagonal = 0.0;
-}
-r.diagonal_goodness = 0.0;
-// Shift source
-r.source = δ_new.x;
-}
-// Extend to new spikes
-γ_hat.extend(μ[γ_hat.len()..].iter().map(|δ_new| {
-RaySet{
-source : δ_new.x,
-rays : [].into(),
-diagonal : 0.0,
-diagonal_goodness : 0.0,
-diagonal_reg_goodness : 0.0
-}
-}));
-// Prune spikes with zero weight. This also moves the marginal differences of corresponding
-// transports from γ_hat to γ_pruned_marginal_diff.
-// TODO: optimise standard prune with swap_remove.
-μ_transported_base.clear();
-let mut i = 0;
-assert_eq!(μ.len(), γ_hat.len());
-while i < μ.len() {
-if μ[i].α == F::ZERO {
-μ.swap_remove(i);
-let r = γ_hat.swap_remove(i);
-μ_transported_base.extend(r.targets().cloned());
-μ_transported_base -= DeltaMeasure{ α : r.non_diagonal_mass(), x : r.source };
-} else {
-i += 1;
-}
 }
 // TODO: how to merge?
 // Update residual
 // Give function value if needed
 state.if_verbose(|| {
 // Plot if so requested
 plotter.plot_spikes(
 format!("iter {} end; {}", state.iteration(), within_tolerances), &d,
-"start".to_string(), Some(&minus_τv),
+"start".to_string(), None::<&A::PreadjointCodomain>, // TODO: Should be Some(&((-τ) * v)), but not implemented
 reg.target_bounds(τ, ε_prev), &μ,
 );
 // Calculate mean inner iterations and reset relevant counters.
 // Return the statistics
 let res = IterInfo {

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs