--- a/src/forward_pdps.rs Tue Aug 04 13:46:29 2026 -0500
+++ b/src/forward_pdps.rs Tue Aug 04 14:12:18 2026 -0500
@@ -3,6 +3,7 @@
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};
@@ -20,6 +21,7 @@
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};
@@ -85,7 +87,7 @@
/// τ = τ_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,
@@ -111,6 +113,7 @@
opKz: &KOpZ,
fnR: &R,
fnH: &H,
+ domain: &mut impl AdmissibleDomain<Loc<N, F>>,
) -> DynResult<MeasureZ<F, Z, N>>
where
F: Float + ToNalgebraRealField,
@@ -237,10 +240,7 @@
// 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.
- // 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| {
stats.merged += prox_penalty.merge_spikes(
&mut μ,
&mut τv,
@@ -249,32 +249,45 @@
ε,
ins,
®,
- 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.
- stats.pruned += prune_with_stats(&mut μ);
+ 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 {
- // Do extra weight optimisation step heuristic
- let Pair(v, _z_tmp) = f.differential(Pair(&μ, &z));
- τv = v * τ;
- prox_penalty.reweigh(
- &mut μ,
- &mut τv,
- τ,
- ε,
- &config.insertion,
- ®,
- &state,
- &mut stats,
- )?;
- //(z, y) = zy_step(z, y, z_tmp);
- }
+ for _ in 0..config.insertion.extra_weight_optimisation_steps {
+ // Do extra weight optimisation step heuristic
+ let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z));
+ τv = v * τ;
+ prox_penalty.reweigh(
+ &mut μ,
+ &mut τv,
+ τ,
+ ε,
+ &config.insertion,
+ ®,
+ &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 μ);
}
@@ -305,13 +318,79 @@
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,
@@ -321,6 +400,7 @@
(μ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,
@@ -345,7 +425,7 @@
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,
@@ -356,5 +436,48 @@
&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,
+ )
+}