pointsource_algs: comparison src/sliding

-:4f468d35fa29
+:7a8a55fd41c0
 */
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
-use crate::measures::{DiscreteMeasure, Radon, RNDM};
+use crate::measures::{DiscreteMeasure, RNDM};
 use crate::plot::Plotter;
 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
 use crate::regularisation::SlidingRegTerm;
-use crate::sliding_fb::{
+use crate::sliding_fb::{SlidingFBConfig, Transport, TransportConfig, TransportStepLength};
-aposteriori_transport, initial_transport, SlidingFBConfig, TransportConfig, TransportStepLength,
-};
 use crate::types::*;
 use alg_tools::convex::{Conjugable, Prox, Zero};
 use alg_tools::direct_product::Pair;
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::ClosedEuclidean;
 use alg_tools::linops::{
 BoundedLinear, IdOp, SimplyAdjointable, StaticEuclideanOriginGenerator, ZeroOp, AXPY, GEMV,
 };
 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Instance};
 use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::{Norm, L2};
+use alg_tools::norms::L2;
 use anyhow::ensure;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
-//use colored::Colorize;
-//use nalgebra::{DVector, DMatrix};
-use std::iter::Iterator;
 /// Settings for [`pointsource_sliding_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 #[serde(default)]
 pub struct SlidingPDPSConfig<F: Float> {
 );*/
 config.transport.check()?;
 // Initialise iterates
 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
-let mut γ1 = DiscreteMeasure::new();
+let mut γ = Transport::new();
 //let zero_z = z.similar_origin();
 // Set up parameters
 // TODO: maybe this PairNorm doesn't make sense here?
 // let opAnorm = opA.opnorm_bound(PairNorm(Radon, L2, L2), L2);
 ensure!(β < 1.0);
 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_F + ℓ]) ∫ c_2 dγ for κ defined as:
 let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM);
 //  The factor two in the manuscript disappears due to the definition of 𝚹 being
 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2.
-let (maybe_ℓ_F, maybe_transport_lip) = f.curvature_bound_components(config.guess);
-let transport_lip = maybe_transport_lip?;
+let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) {
-let calculate_θ = |ℓ_F, max_transport| {
+(_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0),
-let ℓ_r = transport_lip * max_transport;
+(maybe_ℓ_F, Ok(transport_lip)) => {
-config.transport.θ0 / (τ * (ℓ + ℓ_F + ℓ_r) + κ * bigθ * max_transport)
+let calculate_θτ = move |ℓ_F, max_transport| {
-};
+let ℓ_r = transport_lip * max_transport;
-let mut θ_or_adaptive = match maybe_ℓ_F {
+config.transport.θ0 / ((ℓ + ℓ_F + ℓ_r) + κ * bigθ * max_transport / τ)
-// We assume that the residual is decreasing.
+};
-Ok(ℓ_F) => TransportStepLength::AdaptiveMax {
+match maybe_ℓ_F {
-l: ℓ_F, // TODO: could estimate computing the real reesidual
+Ok(ℓ_F) => TransportStepLength::AdaptiveMax {
-max_transport: 0.0,
+l: ℓ_F, // TODO: could estimate computing the real reesidual
-g: calculate_θ,
+max_transport: 0.0,
-},
+g: calculate_θτ,
-Err(_) => {
+},
-TransportStepLength::FullyAdaptive {
+Err(_) => TransportStepLength::FullyAdaptive {
-l: F::EPSILON, max_transport: 0.0, g: calculate_θ
+l: F::EPSILON, // Start with something very small to estimate differentials
+max_transport: 0.0,
+g: calculate_θτ,
+},
 }
 }
 };
 // Acceleration is not currently supported
 // let γ = dataterm.factor_of_strong_convexity();
 // This is much easier with K_μ = 0, which is the only reason why are enforcing it.
 // TODO: Write a version of initial_transport that can deal with K_μ ≠ 0.
 //dbg!(&μ);
-let (μ_base_masses, mut μ_base_minus_γ0) =
+γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v);
-initial_transport(&mut γ1, &mut μ, τ, &mut θ_or_adaptive, v);
+let mut attempts = 0;
 // 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, mut τv̆, z_new) = 'adapt_transport: loop {
+let (maybe_d, _within_tolerances, mut τv̆, z_new, μ̆) = 'adapt_transport: loop {
+// Set initial guess for μ=μ^{k+1}.
+γ.μ̆_into(&mut μ);
+let μ̆ = μ.clone();
 // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b)
-// let residual_μ̆ =
+let Pair(mut τv̆, τz̆) = f.differential(Pair(&μ̆, &z)) * τ;
-//     calculate_residual2(Pair(&γ1, &z), Pair(&μ_base_minus_γ0, &zero_z), opA, b);
-// let Pair(mut τv̆, τz̆) = opA.preadjoint().apply(residual_μ̆ * τ);
-// TODO: might be able to optimise the measure sum working as calculate_residual2 above.
-let Pair(mut τv̆, τz̆) = f.differential(Pair(&γ1 + &μ_base_minus_γ0, &z)) * τ;
 // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0);
 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
 let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh(
 &mut μ,
 &mut τv̆,
-&γ1,
-Some(&μ_base_minus_γ0),
 τ,
 ε,
 &config.insertion,
 &reg,
 &state,
 let mut z_new = τz̆;
 opKz_adj.gemv(&mut z_new, -σ_p, &y, -σ_p / τ);
 z_new = fnR.prox(σ_p, z_new + &z);
 // A posteriori transport adaptation.
-if aposteriori_transport(
+if γ.aposteriori_transport(
-&mut γ1,
+&μ,
-&mut μ,
+&μ̆,
-&mut μ_base_minus_γ0,
+&mut τv̆,
-&μ_base_masses,
 Some(z_new.dist2(&z)),
 ε,
 &config.transport,
+&mut attempts,
 ) {
-break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new);
+break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆);
 }
 };
-stats.untransported_fraction = Some({
+γ.get_transport_stats(&mut stats, &μ);
-assert_eq!(μ_base_masses.len(), γ1.len());
-let (a, b) = stats.untransported_fraction.unwrap_or((0.0, 0.0));
-let source = μ_base_masses.iter().map(|v| v.abs()).sum();
-(a + μ_base_minus_γ0.norm(Radon), b + source)
-});
-stats.transport_error = Some({
-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.
 // This crucially expects the merge routine to be stable with respect to spike locations,
 // and not to performing any pruning. That is be to done below simultaneously for γ.
-let ins = &config.insertion;
+if config.insertion.merge_now(&state) {
-if ins.merge_now(&state) {
 stats.merged += prox_penalty.merge_spikes_no_fitness(
 &mut μ,
 &mut τv̆,
-&γ1,
+&μ̆,
-Some(&μ_base_minus_γ0),
 τ,
 ε,
-ins,
+&config.insertion,
 &reg,
-//Some(|μ̃ : &RNDM<N, F>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()),
 );
 }
-// Prune spikes with zero weight. To maintain correct ordering between μ and γ1, also the
+γ.prune_compat(&mut μ, &mut stats);
-// 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());
-if μ_new.len() != μ.len() {
-let mut μ_iter = μ.iter_spikes();
-γ1.prune_by(|_| μ_iter.next().unwrap().α != F::ZERO);
-stats.pruned += μ.len() - μ_new.len();
-μ = μ_new;
-}
 // Do dual update
 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0);    // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
 opKz.gemv(&mut y, σ_d * (1.0 + ω), &z_new, 1.0);
 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs