--- a/src/forward_pdps.rs Fri Aug 14 10:38:27 2026 -0500
+++ b/src/forward_pdps.rs Fri Oct 02 12:03:25 2026 -0500
@@ -11,6 +11,7 @@
use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
use crate::regularisation::RegTerm;
use crate::types::*;
+use alg_tools::bounds::MinMaxMapping;
use alg_tools::convex::{Conjugable, Prox, Zero};
use alg_tools::direct_product::Pair;
use alg_tools::error::DynResult;
@@ -22,11 +23,12 @@
};
use alg_tools::mapping::{DifferentiableMapping, Instance};
use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::norms::L2;
+use alg_tools::norms::{Norm, L2};
use anyhow::ensure;
-use measures::DeltaMeasure;
+use measures::{DeltaMeasure, Radon};
use numeric_literals::replace_float_literals;
use serde::{Deserialize, Serialize};
+//use std::ops::Mul;
/// Settings for [`pointsource_forward_pdps_pair`].
#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
@@ -123,9 +125,9 @@
U: Space + Clone,
I: AlgIteratorFactory<IterInfo<F>>,
Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
- S: ClosedMul<F>,
+ S: ClosedMul<F> + MinMaxMapping<U, F>,
for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
- DiscreteMeasure<U, F>: SpikeMerging<F>,
+ DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
Reg: RegTerm<U, F>,
P: ProxPenalty<U, S, Reg, F>,
KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
@@ -134,8 +136,8 @@
KOpZ::SimpleAdjoint: GEMV<F, Y>,
Y: ClosedEuclidean<F>,
for<'b> &'b Y: Instance<Y>,
- Z: ClosedEuclidean<F> + AXPY,
- for<'b> &'b Z: Instance<Z>,
+ Z: ClosedEuclidean<F> + AXPY + Clone,
+ for<'b> &'b Z: Instance<Z>, // compiler overflow + Mul<F, Output = Z>,
R: Prox<Z, Codomain = F>,
H: Conjugable<Y, F, Codomain = F>,
for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -158,9 +160,21 @@
let idOpZ = IdOp::new();
let opKz_adj = opKz.adjoint();
let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
+ let l_z_prime = if l_z == 0.0 {
+ 1.0
+ } else {
+ ensure!(config.σp0 < 1.0);
+ l_z
+ };
+ let l_prime = if l == 0.0 {
+ 1.0
+ } else {
+ ensure!(config.τ0 < 1.0);
+ l
+ };
let (τ, σ_p, σ_d) = if is_fb {
- (config.τ0 / l, config.σp0 / l_z, 0.0)
+ (config.τ0 / l_prime, config.σp0 / l_z_prime, 0.0)
} else {
ensure!(
//config.τ0 < 1.0 &&
@@ -168,10 +182,10 @@
"Invalid step length parameters"
);
let σ_d = config.σd0 / nKz;
- let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
+ let σ_p = config.σp0 / (l_z_prime + config.σd0 * nKz);
let φ = 1.0 - config.σp0;
- let a = 1.0 - σ_p * l_z;
- let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
+ let a = 1.0 - σ_p * l_z_prime;
+ let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l_prime);
(τ, σ_p, σ_d)
};
@@ -187,24 +201,37 @@
let starH = fnH.conjugate();
// Statistics
- let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
- value: f.apply(Pair(μ, z))
+ let full_stats = |μ: &DiscreteMeasure<U, F>, v: &mut S, z: &Z, dz: &Z, ε, stats| {
+ let value = f.apply(Pair(μ, z))
+ fnR.apply(z)
+ reg.apply(μ)
- + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
- n_spikes: μ.len(),
- ε,
- // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
- ..stats
+ + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z));
+ let discrepancy = reg.discrepancy(v, τ, ε, μ, &config.insertion);
+ let aux_discrepancy = is_fb.then(|| z.dist2(fnR.prox(σ_p, Z::clone(dz) * (-σ_p) + z)));
+ let joint_discrepancy = aux_discrepancy
+ .map(|d| (discrepancy.powi(2) + τ * (2.0 / σ_p - l_z) * d.powi(2)).sqrt());
+ IterInfo {
+ value,
+ n_spikes: μ.len(),
+ discrepancy: Some(discrepancy),
+ μ_norm: μ.norm(Radon),
+ // compiler overflow with the necesssary trait bounds:
+ //aux_discrepancy: is_fb.then(|| z.dist2(fnR.prox(σ_p, dz * (-σ_p) + z))),
+ aux_discrepancy,
+ joint_discrepancy,
+ ε,
+ // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
+ ..stats
+ }
};
let mut stats = IterInfo::new();
- let zy_step = |mut z, mut y, mut z_tmp| {
+ let zy_step = |mut z, mut y, mut dz| {
// Do z variable primal update
- opKz_adj.apply_add(&mut z_tmp, &y);
- z_tmp.axpy(1.0, &z, -σ_p);
- let z_new = fnR.prox(σ_p, z_tmp);
- //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z);
+ opKz_adj.apply_add(&mut dz, &y);
+ dz.axpy(1.0, &z, -σ_p);
+ let z_new = fnR.prox(σ_p, dz);
+ //let z_new = fnR.prox(σ_p, dz * (-σ_p) + &z);
// Do dual update
z.axpy(1.0 + ω, &z_new, -ω);
//z = (z - &z_new) * (-ω) + &z_new;
@@ -213,11 +240,13 @@
(z_new, y)
};
+ // Calculate differential of smooth part
+ let Pair(mut v, mut dz) = f.differential(Pair(&μ, &z));
+ let mut τv = v * τ;
+
// Run the algorithm
- for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) {
+ for state in iterator.iter_init(|| full_stats(&μ, &mut τv, &z, &dz, ε, stats.clone())) {
// Calculate initial transport
- let Pair(v, z_tmp) = f.differential(Pair(&μ, &z));
- let mut τv = v * τ;
let μ_base = μ.clone();
let μ_base_len = μ_base.len();
@@ -233,7 +262,7 @@
&mut stats,
)?;
- (z, y) = zy_step(z, y, z_tmp);
+ (z, y) = zy_step(z, y, dz);
stats.inserted += μ.len() - μ_base.len();
@@ -263,7 +292,7 @@
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));
+ let Pair(v, mut _dz) = f.differential(Pair(&μ, &z));
τv = v * τ;
prox_penalty.reweigh(
&mut μ,
@@ -276,9 +305,9 @@
&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 dz, &y);
+ // dz.axpy(1.0, &z, -σ_p);
+ // let z_new = fnR.prox(σ_p, dz);
// z.axpy(1.0 + ω, &z_new, -ω);
// //z = (z - &z_new) * (-ω) + &z_new;
// opKz.gemv(&mut y, σ_d, z, 1.0);
@@ -300,9 +329,20 @@
let iter = state.iteration();
stats.this_iters += 1;
+ // Update differential of smooth part
+ Pair(v, dz) = f.differential(Pair(&μ, &z));
+ τv = v * τ;
+
state.if_verbose(|| {
plotter.plot_spikes(iter, None, Some(&τv), &μ);
- full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new()))
+ full_stats(
+ &μ,
+ &mut τv,
+ &z,
+ &dz,
+ ε,
+ std::mem::replace(&mut stats, IterInfo::new()),
+ )
});
// Update main tolerance for next iteration
@@ -354,10 +394,10 @@
//for<'a> &'a U: Instance<U>,
I: AlgIteratorFactory<IterInfo<F>>,
Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
- S: ClosedMul<F>,
+ S: ClosedMul<F> + MinMaxMapping<U, F>,
for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
//Pair<S, Z>: ClosedMul<F>,
- DiscreteMeasure<U, F>: SpikeMerging<F>,
+ DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
Reg: RegTerm<U, F>,
P: ProxPenalty<U, S, Reg, F>,
KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
@@ -366,8 +406,8 @@
KOpZ::SimpleAdjoint: GEMV<F, Y>,
Y: ClosedEuclidean<F>,
for<'b> &'b Y: Instance<Y>,
- Z: ClosedEuclidean<F> + AXPY,
- for<'b> &'b Z: Instance<Z>,
+ Z: ClosedEuclidean<F> + AXPY + Clone, // clone not neede if Mul worked below
+ for<'b> &'b Z: Instance<Z>, // compiler overflow + Mul<F, Output = Z>,
R: Prox<Z, Codomain = F>,
H: Conjugable<Y, F, Codomain = F>,
for<'b> H::Conjugate<'b>: Prox<Y>,
@@ -412,13 +452,13 @@
for<'a> &'a U: Instance<U>,
I: AlgIteratorFactory<IterInfo<F>>,
Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
- S: ClosedMul<F>,
- DiscreteMeasure<U, F>: SpikeMerging<F>,
+ S: ClosedMul<F> + MinMaxMapping<U, F>,
+ DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
Reg: RegTerm<U, F>,
P: ProxPenalty<U, S, Reg, F>,
for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
Z: ClosedEuclidean<F> + AXPY + Clone,
- for<'b> &'b Z: Instance<Z>,
+ for<'b> &'b Z: Instance<Z>, // compiler overflow + Mul<F, Output = Z>,
R: Prox<Z, Codomain = F>,
Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
{
@@ -464,13 +504,13 @@
for<'a> &'a U: Instance<U>,
I: AlgIteratorFactory<IterInfo<F>>,
Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
- S: ClosedMul<F>,
- DiscreteMeasure<U, F>: SpikeMerging<F>,
+ S: ClosedMul<F> + MinMaxMapping<U, F>,
+ DiscreteMeasure<U, F>: SpikeMerging<F> + Norm<Radon, F>,
Reg: RegTerm<U, F>,
P: ProxPenalty<U, S, Reg, F>,
for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
Z: ClosedEuclidean<F> + AXPY + Clone,
- for<'b> &'b Z: Instance<Z>,
+ for<'b> &'b Z: Instance<Z>, // compiler overflow + Mul<F, Output = Z>,
R: Prox<Z, Codomain = F>,
Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
{