pointsource_algs: comparison src/sliding

-:c5d8bd1a7728
+:6316d68b58af
 };
 use crate::fb::*;
 use crate::regularisation::SlidingRegTerm;
 use crate::dataterm::{
 L2Squared,
-//DataTerm,
+DataTerm,
 calculate_residual,
 calculate_residual2,
 };
 //use crate::transport::TransportLipschitz;
 pub struct TransportConfig<F : Float> {
 /// Transport step length $θ$ normalised to $(0, 1)$.
 pub θ0 : F,
 /// Factor in $(0, 1)$ for decreasing transport to adapt to tolerance.
 pub adaptation : F,
-/// Transport tolerance wrt. ω
+/// A priori transport tolerance multiplier (C_pri)
-pub tolerance_ω : F,
+pub tolerance_mult_pri : F,
-/// Transport tolerance wrt. ∇v
+/// A posteriori transport tolerance multiplier (C_pos)
-pub tolerance_dv : F,
+pub tolerance_mult_pos : F,
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl <F : Float> TransportConfig<F> {
 /// Check that the parameters are ok. Panics if not.
 pub fn check(&self) {
 assert!(self.θ0 > 0.0);
 assert!(0.0 < self.adaptation && self.adaptation < 1.0);
-assert!(self.tolerance_dv > 0.0);
+assert!(self.tolerance_mult_pri > 0.0);
-assert!(self.tolerance_ω > 0.0);
+assert!(self.tolerance_mult_pos > 0.0);
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
 impl<F : Float> Default for TransportConfig<F> {
 fn default() -> Self {
 TransportConfig {
-θ0 : 0.01,
+θ0 : 0.4,
 adaptation : 0.9,
-tolerance_ω : 1000.0, // TODO: no idea what this should be
+tolerance_mult_pos : 100.0,
-tolerance_dv : 1000.0, // TODO: no idea what this should be
+tolerance_mult_pri : 1000.0,
 }
 }
 }
 /// Settings for [`pointsource_sliding_fb_reg`].
 }
 },
 FullyAdaptive{ l : ref mut adaptive_ℓ_v, ref mut max_transport, g : ref calculate_θ } => {
 *max_transport = max_transport.max(γ1.norm(Radon));
 let mut θ = calculate_θ(*adaptive_ℓ_v, *max_transport);
-loop {
+// Do two runs through the spikes to update θ, breaking if first run did not cause
-let θτ = τ * θ;
+// a change.
+for _i in 0..=1 {
+let mut changes = false;
 for (δ, ρ) in izip!(μ.iter_spikes(), γ1.iter_spikes_mut()) {
 let dv_x = v.differential(&δ.x);
-ρ.x = δ.x - &dv_x * (ρ.α.signum() * θτ);
+let g = &dv_x * (ρ.α.signum() * θ * τ);
-// Estimate Lipschitz factor of ∇v
+ρ.x = δ.x - g;
-let this_ℓ_v = (dv_x - v.differential(&ρ.x)).norm2();
+let n = g.norm2();
-*adaptive_ℓ_v = adaptive_ℓ_v.max(this_ℓ_v);
+if n >= F::EPSILON {
+// Estimate Lipschitz factor of ∇v
+let this_ℓ_v = (dv_x - v.differential(&ρ.x)).norm2() / n;
+*adaptive_ℓ_v = adaptive_ℓ_v.max(this_ℓ_v);
+θ = calculate_θ(*adaptive_ℓ_v, *max_transport);
+changes = true
+}
 }
-let new_θ = calculate_θ(*adaptive_ℓ_v / tconfig.adaptation, *max_transport);
+if !changes {
-if new_θ <= θ {
 break
 }
-θ = new_θ;
 }
 }
 }
 // 2. Adjust transport mass, if needed.
 let nr =γ1.norm(Radon);
 let n = τ * 2.0 * opAnorm * (opAapply(&*γ1)-opAapply(&*μ)).norm2();
 if n <= 0.0 || nr <= 0.0 {
 break
 }
-let reduction_needed = nr - (ε * tconfig.tolerance_dv / n);
+let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n);
 if reduction_needed <= 0.0 {
 break
 }
 let (min_nonzero, n_nonzero) = γ1.iter_masses()
 .map(|α| α.abs())
 let na = a.norm2();
 let n = τ * 2.0 * opAnorm * na;
 if n <= 0.0 || nr <= 0.0 {
 break
 }
-let reduction_needed = nr - (ε * tconfig.tolerance_dv / n);
+let reduction_needed = nr - (ε * tconfig.tolerance_mult_pri / n);
 if reduction_needed <= 0.0 {
 break
 }
 let mut max_d = 0.0;
 let mut max_d_ind = 0;
 pub(crate) fn aposteriori_transport<F, const N : usize>(
 γ1 : &mut RNDM<F, N>,
 μ : &mut RNDM<F, N>,
 μ_base_minus_γ0 : &mut RNDM<F, N>,
 μ_base_masses : &Vec<F>,
+extra : Option<F>,
 ε : F,
 tconfig : &TransportConfig<F>
 ) -> bool
 where F : Float + ToNalgebraRealField {
 // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z).
 //    through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ^k-(π_♯^1-π_♯^0)γ^{k+1},
 //    which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient.
 let nγ = γ1.norm(Radon);
-let nΔ = μ_base_minus_γ0.norm(Radon) + μ.dist_matching(&γ1);
+let nΔ = μ_base_minus_γ0.norm(Radon) + μ.dist_matching(&γ1) + extra.unwrap_or(0.0);
-let t = ε * tconfig.tolerance_ω;
+let t = ε * tconfig.tolerance_mult_pos;
 if nγ*nΔ > t {
 // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1,
 // this will guarantee that eventually ‖γ‖ decreases sufficiently that we
 // will not enter here.
 *γ1 *= tconfig.adaptation * t / ( nγ * nΔ );
 let calculate_θ = |ℓ_v, _| config.transport.θ0 / (τ*(ℓ + ℓ_v));
 let mut θ_or_adaptive = match opA.preadjoint().value_diff_unit_lipschitz_factor() {
 // We only estimate w (the uniform Lipschitz for of v), if we also estimate ℓ_v
 // (the uniform Lipschitz factor of ∇v).
 // We assume that the residual is decreasing.
-Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * residual.norm2(), 0.0)),
+Some(ℓ_v0) => TransportStepLength::Fixed(calculate_θ(ℓ_v0 * b.norm2(), 0.0)),
 None => TransportStepLength::FullyAdaptive {
-l : 0.0,
+l : 10.0 * F::EPSILON, // Start with something very small to estimate differentials
 max_transport : 0.0,
 g : calculate_θ
 },
 };
 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
 v, &config.transport,
 );
 // Solve finite-dimensional subproblem several times until the dual variable for the
 // regularisation term conforms to the assumptions made for the transport above.
-let (maybe_d, _within_tolerances, τv̆) = 'adapt_transport: loop {
+let (maybe_d, _within_tolerances, mut τv̆) = 'adapt_transport: loop {
 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b)
 let residual_μ̆ = calculate_residual2(&γ1, &μ_base_minus_γ0, opA, b);
 let mut τv̆ = opA.preadjoint().apply(residual_μ̆ * τ);
 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
 );
 // A posteriori transport adaptation.
 if aposteriori_transport(
 &mut γ1, &mut μ, &mut μ_base_minus_γ0, &μ_base_masses,
+None,
 ε, &config.transport
 ) {
 break 'adapt_transport (maybe_d, within_tolerances, τv̆)
 }
 };
 assert_eq!(μ_base_masses.len(), γ1.len());
 let (a, b) = stats.transport_error.unwrap_or((0.0, 0.0));
 (a + μ.dist_matching(&γ1), b + γ1.norm(Radon))
 });
-// // Merge spikes.
+// Merge spikes.
-// // This expects the prune below to prune γ.
+// This crucially expects the merge routine to be stable with respect to spike locations,
-// // TODO: This may not work correctly in all cases.
+// and not to performing any pruning. That is be to done below simultaneously for γ.
-// let ins = &config.insertion;
+let ins = &config.insertion;
-// if ins.merge_now(&state) {
+if ins.merge_now(&state) {
-//     if let SpikeMergingMethod::None = ins.merging {
+stats.merged += prox_penalty.merge_spikes(
-//     } else {
+&mut μ, &mut τv̆, &γ1, Some(&μ_base_minus_γ0), τ, ε, ins, &reg,
-//         stats.merged += μ.merge_spikes(ins.merging, |μ_candidate| {
+Some(|μ̃ : &RNDM<F, N>| L2Squared.calculate_fit_op(μ̃, opA, b)),
-//             let ν = μ_candidate.sub_matching(&γ1)-&μ_base_minus_γ0;
+);
-//             let mut d = &τv̆ + op𝒟.preapply(ν);
+}
-//             reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, ins)
-//         });
-//     }
-// }
 // Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the
 // latter needs to be pruned when μ is.
 // TODO: This could do with a two-vector Vec::retain to avoid copies.
 let μ_new = DiscreteMeasure::from_iter(μ.iter_spikes().filter(|δ| δ.α != F::ZERO).cloned());

comparison: src/sliding_fb.rs

src/sliding_fb.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_fb.rs

src/sliding_fb.rs