pointsource_algs: comparison src/forward

-:f91f382ba4ae
+:b9818eb24b47
 use crate::measures::DiscreteMeasure;
 use crate::plot::Plotter;
 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;
 use alg_tools::euclidean::ClosedEuclidean;
 use alg_tools::iterate::AlgIteratorFactory;
 BoundedLinear, IdOp, SimplyAdjointable, Space, StaticEuclideanOriginGenerator, ZeroOp, AXPY,
 GEMV,
 };
 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)]
 #[serde(default)]
 pub struct ForwardPDPSConfig<F: Float> {
 where
 F: Float + ToNalgebraRealField,
 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>
 + GEMV<F, Z>
 + SimplyAdjointable<Z, Y, AdjointCodomain = Z>,
 KOpZ::SimpleAdjoint: GEMV<F, Y>,
 Y: ClosedEuclidean<F>,
 for<'b> &'b Y: Instance<Y>,
-Z: ClosedEuclidean<F> + AXPY,
+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>,
 H: Conjugable<Y, F, Codomain = F>,
 for<'b> H::Conjugate<'b>: Prox<Y>,
 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {
 let nKz = opKz.opnorm_bound(L2, L2)?;
 let is_fb = nKz == 0.0;
 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 &&
 config.σp0 < 1.0 && config.σd0 > 0.0 && config.σp0 * config.σd0 <= 1.0,
 "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 a = 1.0 - σ_p * l_z_prime;
-let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
+let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l_prime);
 (τ, σ_p, σ_d)
 };
 // Acceleration is not currently supported
 // let γ = dataterm.factor_of_strong_convexity();
 let mut ε = tolerance.initial();
 let starH = fnH.conjugate();
 // Statistics
-let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
+let full_stats = |μ: &DiscreteMeasure<U, F>, v: &mut S, z: &Z, dz: &Z, ε, stats| {
-value: f.apply(Pair(μ, z))
+let value = f.apply(Pair(μ, z))
 + fnR.apply(z)
 + reg.apply(μ)
-+ fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
++ fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z));
-n_spikes: μ.len(),
+let discrepancy = reg.discrepancy(v, τ, ε, μ, &config.insertion);
-ε,
+let aux_discrepancy = is_fb.then(|| z.dist2(fnR.prox(σ_p, Z::clone(dz) * (-σ_p) + z)));
-// postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
+let joint_discrepancy = aux_discrepancy
-..stats
+.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);
+opKz_adj.apply_add(&mut dz, &y);
-z_tmp.axpy(1.0, &z, -σ_p);
+dz.axpy(1.0, &z, -σ_p);
-let z_new = fnR.prox(σ_p, z_tmp);
+let z_new = fnR.prox(σ_p, dz);
-//let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z);
+//let z_new = fnR.prox(σ_p, dz * (-σ_p) + &z);
 // Do dual update
 z.axpy(1.0 + ω, &z_new, -ω);
 //z = (z - &z_new) * (-ω) + &z_new;
 opKz.gemv(&mut y, σ_d, z, 1.0);
 let y = starH.prox(σ_d, y);
 (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();
 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
 prox_penalty.insert_and_reweigh(
 &reg,
 &state,
 &mut stats,
 )?;
-(z, y) = zy_step(z, y, z_tmp);
+(z, y) = zy_step(z, y, dz);
 stats.inserted += μ.len() - μ_base.len();
 // Merge spikes.
 // This crucially expects the merge routine to be stable with respect to spike locations,
 }
 }
 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 μ,
 &mut τv,
 τ,
 &reg,
 &state,
 &mut stats,
 )?;
-// opKz_adj.apply_add(&mut z_tmp, &y);
+// opKz_adj.apply_add(&mut dz, &y);
-// z_tmp.axpy(1.0, &z, -σ_p);
+// dz.axpy(1.0, &z, -σ_p);
-// let z_new = fnR.prox(σ_p, z_tmp);
+// 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);
 // y = starH.prox(σ_d, y);
 // z = z_new;
 // Give statistics if requested
 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
 ε = tolerance.update(ε, iter);
 }
 F: Float + ToNalgebraRealField,
 U: Space + Clone,
 //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>
 + GEMV<F, Z>
 + SimplyAdjointable<Z, Y, AdjointCodomain = Z>,
 KOpZ::SimpleAdjoint: GEMV<F, Y>,
 Y: ClosedEuclidean<F>,
 for<'b> &'b Y: Instance<Y>,
-Z: ClosedEuclidean<F> + AXPY,
+Z: ClosedEuclidean<F> + AXPY + Clone, // clone not neede if Mul worked below
-for<'b> &'b Z: Instance<Z>,
+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>,
 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
 {
 F: Float + ToNalgebraRealField,
 U: Space + Clone,
 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>,
-DiscreteMeasure<U, F>: SpikeMerging<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>>,
 {
 let opKz: ZeroOp<Z, Loc<0, F>, _, _, F> =
 ZeroOp::new_dualisable(StaticEuclideanOriginGenerator, z.dual_origin());
 F: Float + ToNalgebraRealField,
 U: Space + Clone,
 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>,
-DiscreteMeasure<U, F>: SpikeMerging<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>>,
 {
 pointsource_fb_pair_dom(
 f,

comparison: src/forward_pdps.rs

src/forward_pdps.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_pdps.rs

src/forward_pdps.rs