pointsource_algs: comparison src/sliding

-:51dbb2cc25d9
+:8d1c959432c8
 /*!
 Solver for the point source localisation problem using a sliding
 primal-dual proximal splitting method.
 */
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::forward_model::{BoundedCurvature, BoundedCurvatureGuess};
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DiscreteMeasure, RNDM};
 use crate::plot::Plotter;
 use crate::regularisation::SlidingRegTerm;
 use crate::sliding_fb::{SlidingFBConfig, Transport, TransportConfig, TransportProxPenalty};
 use crate::types::*;
 use alg_tools::bounds::MinMaxMapping;
 use alg_tools::convex::{Conjugable, Prox, Zero};
-use alg_tools::norms::L2;
 use alg_tools::direct_product::Pair;
 use alg_tools::error::DynResult;
 use alg_tools::euclidean::ClosedEuclidean;
 use alg_tools::iterate::AlgIteratorFactory;
 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::L2;
 use anyhow::ensure;
+use measures::DeltaMeasure;
 use numeric_literals::replace_float_literals;
 use serde::{Deserialize, Serialize};
 /// Settings for [`pointsource_sliding_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 (μ0, mut z, mut y): (Option<RNDM<N, F>>, Z, Y),
 //opKμ : KOpM,
 opKz: &KOpZ,
 fnR: &R,
 fnH: &H,
-domain: Option<&P::Dom>,
+domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
 F: Float + ToNalgebraRealField,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 z = z_new;
 // 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 γ.
-config.insertion.if_merge_now(&state, |ins| {
+let did_merge = config.insertion.if_merge_now(&state, |ins| {
-let m = prox_penalty.merge_spikes(
+stats.merged += prox_penalty.merge_spikes(
 &mut μ,
 &mut τv̆,
 &μ̆,
 τ,
 ε,
 ins,
 &reg,
 Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)),
 );
-if m > 0 {
+γ.prune_compat(&mut μ, &mut stats);
-stats.merged += m;
+if config.insertion.random_insertions {
+if let Some(x) = domain.random_point() {
+μ += DeltaMeasure { x, α: 0.0 };
+stats.inserted += 1;
+}
 }
-});
+for _ in 0..config.insertion.extra_weight_optimisation_steps {
-γ.prune_compat(&mut μ, &mut stats);
+// Do extra weight optimisation step heuristic
+let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z));
-for _ in 0..config.insertion.extra_weight_optimisation_steps {
+τv̆ = v * τ;
-// Do extra weight optimisation step heuristic
+prox_penalty.reweigh(
-let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z));
+&mut μ,
-τv̆ = v * τ;
+&mut τv̆,
-prox_penalty.reweigh(
+τ,
-&mut μ,
+ε,
-&mut τv̆,
+&config.insertion,
-τ,
+&reg,
-ε,
+&state,
-&config.insertion,
+&mut stats,
-&reg,
+)?;
-&state,
-&mut stats,
+// opKz_adj.apply_add(&mut z_tmp, &y);
-)?;
+// z_tmp.axpy(1.0, &z, -σ_p);
+// let z_new = fnR.prox(σ_p, z_tmp);
-// opKz_adj.apply_add(&mut z_tmp, &y);
+// z.axpy(1.0 + ω, &z_new, -ω);
-// z_tmp.axpy(1.0, &z, -σ_p);
+// //z = (z - &z_new) * (-ω) + &z_new;
-// let z_new = fnR.prox(σ_p, z_tmp);
+// opKz.gemv(&mut y, σ_d, z, 1.0);
-// z.axpy(1.0 + ω, &z_new, -ω);
+// y = starH.prox(σ_d, y);
-// //z = (z - &z_new) * (-ω) + &z_new;
+// z = z_new;
-// opKz.gemv(&mut y, σ_d, z, 1.0);
+}
-// y = starH.prox(σ_d, y);
+Ok(())
-// z = z_new;
+})?;
-}
+if did_merge.is_none() {
-if config.insertion.extra_weight_optimisation_steps > 0 {
 γ.prune_compat(&mut μ, &mut stats);
 }
 // Update step length parameters
 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
 iterator: I,
 plotter: Plot,
 (μ0, z): (Option<RNDM<N, F>>, Z),
 //opKμ : KOpM,
 fnR: &R,
-domain: Option<&P::Dom>
+domain: &mut impl AdmissibleDomain<Loc<N, F>>,
 ) -> DynResult<MeasureZ<F, Z, N>>
 where
 F: Float + ToNalgebraRealField,
 I: AlgIteratorFactory<IterInfo<F>>,
 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>
 plotter,
 (μ0, z, Loc([])),
 &opKz,
 fnR,
 &fnH,
-domain
+domain,
 )
 }
 #[deprecated]
 /// Version of [`pointsource_sliding_pdps_pair`] without domain specification.
 #[replace_float_literals(F::cast_from(literal))]
 pub fn pointsource_sliding_pdps_pair<
 plotter,
 init,
 opKz,
 fnR,
 fnH,
-None,
+&mut FullSpace,
 )
 }
 #[deprecated]
 /// Version of [`pointsource_sliding_fb_pair_dom`] without domain specification
 plotter,
 (μ0, z, Loc([])),
 &opKz,
 fnR,
 &fnH,
-None
+&mut FullSpace,
 )
 }

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/sliding_pdps.rs

src/sliding_pdps.rs