pointsource_algs: comparison src/forward

-:51dbb2cc25d9
+:8d1c959432c8
 /*!
 Solver for the point source localisation problem using a
 primal-dual proximal splitting with a forward step.
 */
+use crate::domain::{AdmissibleDomain, FullSpace};
 use crate::fb::*;
 use crate::measures::merging::SpikeMerging;
 use crate::measures::{DiscreteMeasure, RNDM};
 use crate::plot::Plotter;
 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
 use alg_tools::linops::{BoundedLinear, IdOp, SimplyAdjointable, 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_forward_pdps_pair`].
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 /// part $[σ_p L_z + σ_pσ_d‖K_z‖^2] < 1$. Then with these cohices, we solve
 /// $$
 ///     τ = τ_0 \frac{1 - σ_{p,0}}{(σ_d M (1-σ_p L_z) + (1 - σ_{p,0} L)}.
 /// $$
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_forward_pdps_pair<
+pub fn pointsource_forward_pdps_pair_dom<
 F,
 I,
 S,
 Dat,
 Reg,
 (μ0, mut z, mut y): (Option<RNDM<N, F>>, Z, Y),
 //opKμ : KOpM,
 opKz: &KOpZ,
 fnR: &R,
 fnH: &H,
+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>>,
 stats.inserted += μ.len() - μ_base.len();
 // 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 γ.
-// Merge spikes.
+let did_merge = config.insertion.if_merge_now(&state, |ins| {
-// 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| {
 stats.merged += prox_penalty.merge_spikes(
 &mut μ,
 &mut τv,
 &μ_base,
 τ,
 ε,
 ins,
 &reg,
-Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))+ reg.apply(μ̃)),
+Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)),
 );
-});
+stats.pruned += prune_with_stats(&mut μ);
-// Prune spikes with zero weight.
+if config.insertion.random_insertions {
-stats.pruned += prune_with_stats(&mut μ);
+if let Some(x) = domain.random_point() {
+μ += DeltaMeasure { x, α: 0.0 };
-for _ in 0..config.insertion.extra_weight_optimisation_steps {
+stats.inserted += 1;
-// Do extra weight optimisation step heuristic
+}
-let Pair(v, _z_tmp) = f.differential(Pair(&μ, &z));
+}
-τv = v * τ;
-prox_penalty.reweigh(
+for _ in 0..config.insertion.extra_weight_optimisation_steps {
-&mut μ,
+// Do extra weight optimisation step heuristic
-&mut τv,
+let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z));
-τ,
+τv = v * τ;
-ε,
+prox_penalty.reweigh(
-&config.insertion,
+&mut μ,
-&reg,
+&mut τv,
-&state,
+τ,
-&mut stats,
+ε,
-)?;
+&config.insertion,
-//(z, y) = zy_step(z, y, z_tmp);
+&reg,
-}
+&state,
+&mut stats,
-if config.insertion.extra_weight_optimisation_steps > 0 {
+)?;
-// Prune spikes with zero weight.
+// opKz_adj.apply_add(&mut z_tmp, &y);
+// z_tmp.axpy(1.0, &z, -σ_p);
+// let z_new = fnR.prox(σ_p, z_tmp);
+// z.axpy(1.0 + ω, &z_new, -ω);
+// //z = (z - &z_new) * (-ω) + &z_new;
+// opKz.gemv(&mut y, σ_d, z, 1.0);
+// y = starH.prox(σ_d, y);
+// z = z_new;
+}
+Ok(())
+})?;
+if did_merge.is_none() {
 stats.pruned += prune_with_stats(&mut μ);
 }
 // Update step length parameters
 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
 };
 μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v);
 μ.prune();
 Ok(Pair(μ, z))
+}
+#[deprecated]
+pub fn pointsource_forward_pdps_pair<
+F,
+I,
+S,
+Dat,
+Reg,
+P,
+Z,
+R,
+Y,
+/*KOpM, */ KOpZ,
+H,
+Plot,
+const N: usize,
+>(
+f: &Dat,
+reg: &Reg,
+prox_penalty: &P,
+config: &ForwardPDPSConfig<F>,
+iterator: I,
+plotter: Plot,
+init: (Option<RNDM<N, F>>, Z, Y),
+//opKμ : KOpM,
+opKz: &KOpZ,
+fnR: &R,
+fnH: &H,
+) -> 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>>,
+//Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex…
+S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
+RNDM<N, F>: SpikeMerging<F>,
+Reg: RegTerm<Loc<N, F>, F>,
+P: ProxPenalty<Loc<N, F>, S, Reg, F>,
+for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
++ GEMV<F, Z, Y>
++ SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>,
+KOpZ::SimpleAdjoint: GEMV<F, Y, Z>,
+Y: ClosedEuclidean<F>,
+for<'b> &'b Y: Instance<Y>,
+Z: ClosedEuclidean<F> + AXPY,
+for<'b> &'b Z: Instance<Z>,
+R: Prox<Z, Codomain = F>,
+H: Conjugable<Y, F, Codomain = F>,
+for<'b> H::Conjugate<'b>: Prox<Y>,
+Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+{
+pointsource_forward_pdps_pair_dom(
+f,
+reg,
+prox_penalty,
+config,
+iterator,
+plotter,
+init,
+opKz,
+fnR,
+fnH,
+&mut FullSpace,
+)
 }
 /// Iteratively solve the pointsource localisation with an additional variable
 /// using forward-backward splitting.
 ///
 /// The implementation uses [`pointsource_forward_pdps_pair`] with appropriate dummy
 /// variables, operators, and functions.
 #[replace_float_literals(F::cast_from(literal))]
-pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
+pub fn pointsource_fb_pair_dom<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
 f: &Dat,
 reg: &Reg,
 prox_penalty: &P,
 config: &FBConfig<F>,
 iterator: I,
 plotter: Plot,
 (μ0, z): (Option<RNDM<N, F>>, Z),
 //opKμ : KOpM,
 fnR: &R,
+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>>,
 // Convert config. We don't implement From (that could be done with the o2o crate), as σd0
 // needs to be chosen in a general case; for the problem of this fucntion, anything is valid.
 let &FBConfig { τ0, σp0, insertion } = config;
 let pdps_config = ForwardPDPSConfig { τ0, σp0, insertion, σd0: 0.0 };
-pointsource_forward_pdps_pair(
+pointsource_forward_pdps_pair_dom(
 f,
 reg,
 prox_penalty,
 &pdps_config,
 iterator,
 plotter,
 (μ0, z, Loc([])),
 &opKz,
 fnR,
 &fnH,
+domain,
 )
 }
+#[deprecated]
+#[replace_float_literals(F::cast_from(literal))]
+pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
+f: &Dat,
+reg: &Reg,
+prox_penalty: &P,
+config: &FBConfig<F>,
+iterator: I,
+plotter: Plot,
+init: (Option<RNDM<N, F>>, Z),
+//opKμ : KOpM,
+fnR: &R,
+) -> 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>>,
+S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
+RNDM<N, F>: SpikeMerging<F>,
+Reg: RegTerm<Loc<N, F>, F>,
+P: ProxPenalty<Loc<N, F>, S, Reg, F>,
+for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
+Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone,
+for<'b> &'b Z: Instance<Z>,
+R: Prox<Z, Codomain = F>,
+Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
+// We should not need to explicitly require this:
+for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
+{
+pointsource_fb_pair_dom(
+f,
+reg,
+prox_penalty,
+config,
+iterator,
+plotter,
+init,
+fnR,
+&mut FullSpace,
+)
+}

comparison: src/forward_pdps.rs

src/forward_pdps.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_pdps.rs

src/forward_pdps.rs