Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_pdps.rs

src/forward_pdps.rs

changeset 70
ed16d0f10d08
parent 66
fe47ad484deb
equal deleted inserted replaced
58:6099ba025aac 70:ed16d0f10d08
1 /*! 1 /*!
2 Solver for the point source localisation problem using a 2 Solver for the point source localisation problem using a
3 primal-dual proximal splitting with a forward step. 3 primal-dual proximal splitting with a forward step.
4 */ 4 */
5 5
6 use crate::fb::*;
7 use crate::measures::merging::SpikeMerging;
8 use crate::measures::{DiscreteMeasure, RNDM};
9 use crate::plot::Plotter;
10 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
11 use crate::regularisation::RegTerm;
12 use crate::types::*;
13 use alg_tools::convex::{Conjugable, Prox, Zero};
14 use alg_tools::direct_product::Pair;
15 use alg_tools::error::DynResult;
16 use alg_tools::euclidean::ClosedEuclidean;
17 use alg_tools::iterate::AlgIteratorFactory;
18 use alg_tools::linops::{BoundedLinear, IdOp, SimplyAdjointable, ZeroOp, AXPY, GEMV};
19 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Instance};
20 use alg_tools::nalgebra_support::ToNalgebraRealField;
21 use alg_tools::norms::L2;
22 use anyhow::ensure;
6 use numeric_literals::replace_float_literals; 23 use numeric_literals::replace_float_literals;
7 use serde::{Serialize, Deserialize}; 24 use serde::{Deserialize, Serialize};
8
9 use alg_tools::iterate::AlgIteratorFactory;
10 use alg_tools::euclidean::Euclidean;
11 use alg_tools::mapping::{Mapping, DifferentiableRealMapping, Instance};
12 use alg_tools::norms::Norm;
13 use alg_tools::direct_product::Pair;
14 use alg_tools::nalgebra_support::ToNalgebraRealField;
15 use alg_tools::linops::{
16 BoundedLinear, AXPY, GEMV, Adjointable, IdOp,
17 };
18 use alg_tools::convex::{Conjugable, Prox};
19 use alg_tools::norms::{L2, PairNorm};
20
21 use crate::types::*;
22 use crate::measures::{DiscreteMeasure, Radon, RNDM};
23 use crate::measures::merging::SpikeMerging;
24 use crate::forward_model::{
25 ForwardModel,
26 AdjointProductPairBoundedBy,
27 };
28 use crate::plot::{
29 SeqPlotter,
30 Plotting,
31 PlotLookup
32 };
33 use crate::fb::*;
34 use crate::regularisation::RegTerm;
35 use crate::dataterm::calculate_residual;
36 25
37 /// Settings for [`pointsource_forward_pdps_pair`]. 26 /// Settings for [`pointsource_forward_pdps_pair`].
38 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 27 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
39 #[serde(default)] 28 #[serde(default)]
40 pub struct ForwardPDPSConfig<F : Float> { 29 pub struct ForwardPDPSConfig<F: Float> {
41 /// Primal step length scaling. 30 /// Overall primal step length scaling.
42 pub τ0 : F, 31 pub τ0: F,
43 /// Primal step length scaling. 32 /// Primal step length scaling for additional variable.
44 pub σp0 : F, 33 pub σp0: F,
45 /// Dual step length scaling. 34 /// Dual step length scaling for additional variable.
46 pub σd0 : F, 35 ///
36 /// Taken zero for [`pointsource_fb_pair`].
37 pub σd0: F,
47 /// Generic parameters 38 /// Generic parameters
48 pub insertion : FBGenericConfig<F>, 39 pub insertion: InsertionConfig<F>,
49 } 40 }
50 41
51 #[replace_float_literals(F::cast_from(literal))] 42 #[replace_float_literals(F::cast_from(literal))]
52 impl<F : Float> Default for ForwardPDPSConfig<F> { 43 impl<F: Float> Default for ForwardPDPSConfig<F> {
53 fn default() -> Self { 44 fn default() -> Self {
54 ForwardPDPSConfig { 45 ForwardPDPSConfig { τ0: 0.99, σd0: 0.05, σp0: 0.99, insertion: Default::default() }
55 τ0 : 0.99,
56 σd0 : 0.05,
57 σp0 : 0.99,
58 insertion : Default::default()
59 }
60 } 46 }
61 } 47 }
62 48
63 type MeasureZ<F, Z, const N : usize> = Pair<RNDM<F, N>, Z>; 49 type MeasureZ<F, Z, const N: usize> = Pair<RNDM<N, F>, Z>;
64 50
65 /// Iteratively solve the pointsource localisation with an additional variable 51 /// Iteratively solve the pointsource localisation with an additional variable
66 /// using primal-dual proximal splitting with a forward step. 52 /// using primal-dual proximal splitting with a forward step.
53 ///
54 /// The problem is
55 /// $$
56 /// \min_{μ, z}~ F(μ, z) + R(z) + H(K_z z) + Q(μ),
57 /// $$
58 /// where
59 /// * The data term $F$ is given in `f`,
60 /// * the measure (Radon or positivity-constrained Radon) regulariser in $Q$ is given in `reg`,
61 /// * the functions $R$ and $H$ are given in `fnR` and `fnH`, and
62 /// * the operator $K_z$ in `opKz`.
63 ///
64 /// This is dualised to
65 /// $$
66 /// \min_{μ, z}\max_y~ F(μ, z) + R(z) + ⟨K_z z, y⟩ + Q(μ) - H^*(y).
67 /// $$
68 ///
69 /// The algorithm is controlled by:
70 /// * the proximal penalty in `prox_penalty`.
71 /// * the initial iterates in `z`, `y`
72 /// * The configuration in `config`.
73 /// * The `iterator` that controls stopping and reporting.
74 /// Moreover, plotting is performed by `plotter`.
75 ///
76 /// The step lengths need to satisfy
77 /// $$
78 /// τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1
79 /// $$ ^^^^^^^^^^^^^^^^^^^^^^^^^
80 /// with $1 > σ_p L_z$ and $1 > τ L$.
81 /// Since we are given “scalings” $τ_0$, $σ_{p,0}$, and $σ_{d,0}$ in `config`, we take
82 /// $σ_d=σ_{d,0}/‖K_z‖$, and $σ_p = σ_{p,0} / (L_z σ_d‖K_z‖)$. This satisfies the
83 /// part $[σ_p L_z + σ_pσ_d‖K_z‖^2] < 1$. Then with these cohices, we solve
84 /// $$
85 /// τ = τ_0 \frac{1 - σ_{p,0}}{(σ_d M (1-σ_p L_z) + (1 - σ_{p,0} L)}.
86 /// $$
67 #[replace_float_literals(F::cast_from(literal))] 87 #[replace_float_literals(F::cast_from(literal))]
68 pub fn pointsource_forward_pdps_pair< 88 pub fn pointsource_forward_pdps_pair<
69 F, I, A, S, Reg, P, Z, R, Y, /*KOpM, */ KOpZ, H, const N : usize 89 F,
90 I,
91 S,
92 Dat,
93 Reg,
94 P,
95 Z,
96 R,
97 Y,
98 /*KOpM, */ KOpZ,
99 H,
100 Plot,
101 const N: usize,
70 >( 102 >(
71 opA : &A, 103 f: &Dat,
72 b : &A::Observable, 104 reg: &Reg,
73 reg : Reg, 105 prox_penalty: &P,
74 prox_penalty : &P, 106 config: &ForwardPDPSConfig<F>,
75 config : &ForwardPDPSConfig<F>, 107 iterator: I,
76 iterator : I, 108 mut plotter: Plot,
77 mut plotter : SeqPlotter<F, N>, 109 (μ0, mut z, mut y): (Option<RNDM<N, F>>, Z, Y),
78 //opKμ : KOpM, 110 //opKμ : KOpM,
79 opKz : &KOpZ, 111 opKz: &KOpZ,
80 fnR : &R, 112 fnR: &R,
81 fnH : &H, 113 fnH: &H,
82 mut z : Z, 114 ) -> DynResult<MeasureZ<F, Z, N>>
83 mut y : Y,
84 ) -> MeasureZ<F, Z, N>
85 where 115 where
86 F : Float + ToNalgebraRealField, 116 F: Float + ToNalgebraRealField,
87 I : AlgIteratorFactory<IterInfo<F, N>>, 117 I: AlgIteratorFactory<IterInfo<F>>,
88 A : ForwardModel< 118 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
89 MeasureZ<F, Z, N>, 119 //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex…
90 F, 120 S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
91 PairNorm<Radon, L2, L2>, 121 RNDM<N, F>: SpikeMerging<F>,
92 PreadjointCodomain = Pair<S, Z>, 122 Reg: RegTerm<Loc<N, F>, F>,
93 > 123 P: ProxPenalty<Loc<N, F>, S, Reg, F>,
94 + AdjointProductPairBoundedBy<MeasureZ<F, Z, N>, P, IdOp<Z>, FloatType=F>, 124 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
95 S: DifferentiableRealMapping<F, N>, 125 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
96 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable> + Instance<A::Observable>, 126 + GEMV<F, Z, Y>
97 PlotLookup : Plotting<N>, 127 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>,
98 RNDM<F, N> : SpikeMerging<F>, 128 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>,
99 Reg : RegTerm<F, N>, 129 Y: ClosedEuclidean<F>,
100 P : ProxPenalty<F, S, Reg, N>, 130 for<'b> &'b Y: Instance<Y>,
101 KOpZ : BoundedLinear<Z, L2, L2, F, Codomain=Y> 131 Z: ClosedEuclidean<F>,
102 + GEMV<F, Z> 132 for<'b> &'b Z: Instance<Z>,
103 + Adjointable<Z, Y, AdjointCodomain = Z>, 133 R: Prox<Z, Codomain = F>,
104 for<'b> KOpZ::Adjoint<'b> : GEMV<F, Y>, 134 H: Conjugable<Y, F, Codomain = F>,
105 Y : AXPY<F> + Euclidean<F, Output=Y> + Clone + ClosedAdd, 135 for<'b> H::Conjugate<'b>: Prox<Y>,
106 for<'b> &'b Y : Instance<Y>, 136 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
107 Z : AXPY<F, Owned=Z> + Euclidean<F, Output=Z> + Clone + Norm<F, L2>,
108 for<'b> &'b Z : Instance<Z>,
109 R : Prox<Z, Codomain=F>,
110 H : Conjugable<Y, F, Codomain=F>,
111 for<'b> H::Conjugate<'b> : Prox<Y>,
112 { 137 {
113
114 // Check parameters 138 // Check parameters
115 assert!(config.τ0 > 0.0 && 139 ensure!(
116 config.τ0 < 1.0 && 140 config.τ0 > 0.0
117 config.σp0 > 0.0 && 141 && config.τ0 < 1.0
118 config.σp0 < 1.0 && 142 && config.σp0 > 0.0
119 config.σd0 > 0.0 && 143 && config.σp0 < 1.0
120 config.σp0 * config.σd0 <= 1.0, 144 && config.σd0 >= 0.0
121 "Invalid step length parameters"); 145 && config.σp0 * config.σd0 <= 1.0,
146 "Invalid step length parameters"
147 );
122 148
123 // Initialise iterates 149 // Initialise iterates
124 let mut μ = DiscreteMeasure::new(); 150 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
125 let mut residual = calculate_residual(Pair(&μ, &z), opA, b);
126 151
127 // Set up parameters 152 // Set up parameters
128 let bigM = 0.0; //opKμ.adjoint_product_bound(prox_penalty).unwrap().sqrt(); 153 let bigM = 0.0; //opKμ.adjoint_product_bound(prox_penalty).unwrap().sqrt();
129 let nKz = opKz.opnorm_bound(L2, L2); 154 let nKz = opKz.opnorm_bound(L2, L2)?;
130 let opIdZ = IdOp::new(); 155 let is_fb = nKz == 0.0;
131 let (l, l_z) = opA.adjoint_product_pair_bound(prox_penalty, &opIdZ).unwrap(); 156 let idOpZ = IdOp::new();
157 let opKz_adj = opKz.adjoint();
158 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?;
132 // We need to satisfy 159 // We need to satisfy
133 // 160 //
134 // τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1 161 // τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1
135 // ^^^^^^^^^^^^^^^^^^^^^^^^^ 162 // ^^^^^^^^^^^^^^^^^^^^^^^^^
136 // with 1 > σ_p L_z and 1 > τ L. 163 // with 1 > σ_p L_z and 1 > τ L.
137 // 164 //
138 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition 165 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition
139 // ^^^^^ < 1. then we solve τ from the rest. 166 // ^^^^^ < 1. then we solve τ from the rest.
140 let σ_d = config.σd0 / nKz; 167 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below.
168 let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz };
141 let σ_p = config.σp0 / (l_z + config.σd0 * nKz); 169 let σ_p = config.σp0 / (l_z + config.σd0 * nKz);
142 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0} 170 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0}
143 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L) 171 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L)
144 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0}) 172 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0})
145 let φ = 1.0 - config.σp0; 173 let φ = 1.0 - config.σp0;
146 let a = 1.0 - σ_p * l_z; 174 let a = 1.0 - σ_p * l_z;
147 let τ = config.τ0 * φ / ( σ_d * bigM * a + φ * l ); 175 let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l);
148 // Acceleration is not currently supported 176 // Acceleration is not currently supported
149 // let γ = dataterm.factor_of_strong_convexity(); 177 // let γ = dataterm.factor_of_strong_convexity();
150 let ω = 1.0; 178 let ω = 1.0;
151 179
152 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled 180 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
155 let mut ε = tolerance.initial(); 183 let mut ε = tolerance.initial();
156 184
157 let starH = fnH.conjugate(); 185 let starH = fnH.conjugate();
158 186
159 // Statistics 187 // Statistics
160 let full_stats = |residual : &A::Observable, μ : &RNDM<F, N>, z : &Z, ε, stats| IterInfo { 188 let full_stats = |μ: &RNDM<N, F>, z: &Z, ε, stats| IterInfo {
161 value : residual.norm2_squared_div2() + fnR.apply(z) 189 value: f.apply(Pair(μ, z))
162 + reg.apply(μ) + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)), 190 + fnR.apply(z)
163 n_spikes : μ.len(), 191 + reg.apply(μ)
192 + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
193 n_spikes: μ.len(),
164 ε, 194 ε,
165 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()), 195 // postprocessing: config.insertion.postprocessing.then(|| μ.clone()),
166 .. stats 196 ..stats
167 }; 197 };
168 let mut stats = IterInfo::new(); 198 let mut stats = IterInfo::new();
169 199
170 // Run the algorithm 200 // Run the algorithm
171 for state in iterator.iter_init(|| full_stats(&residual, &μ, &z, ε, stats.clone())) { 201 for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) {
172 // Calculate initial transport 202 // Calculate initial transport
173 let Pair(mut τv, τz) = opA.preadjoint().apply(residual * τ); 203 let Pair(mut τv, τz) = f.differential(Pair(&μ, &z)) * τ;
174 let μ_base = μ.clone(); 204 let μ_base = μ.clone();
175 205
176 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. 206 // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes.
177 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( 207 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
178 &mut μ, &mut τv, &μ_base, None, 208 &mut μ,
179 τ, ε, &config.insertion, 209 &mut τv,
180 &reg, &state, &mut stats, 210 τ,
181 ); 211 ε,
212 &config.insertion,
213 &reg,
214 &state,
215 &mut stats,
216 )?;
217
218 stats.inserted += μ.len() - μ_base.len();
182 219
183 // Merge spikes. 220 // Merge spikes.
184 // This crucially expects the merge routine to be stable with respect to spike locations, 221 // This crucially expects the merge routine to be stable with respect to spike locations,
185 // and not to performing any pruning. That is be to done below simultaneously for γ. 222 // and not to performing any pruning. That is be to done below simultaneously for γ.
186 // Merge spikes. 223 // Merge spikes.
187 // This crucially expects the merge routine to be stable with respect to spike locations, 224 // This crucially expects the merge routine to be stable with respect to spike locations,
188 // and not to performing any pruning. That is be to done below simultaneously for γ. 225 // and not to performing any pruning. That is be to done below simultaneously for γ.
189 let ins = &config.insertion; 226 let ins = &config.insertion;
190 if ins.merge_now(&state) { 227 if ins.merge_now(&state) {
191 stats.merged += prox_penalty.merge_spikes_no_fitness( 228 stats.merged += prox_penalty.merge_spikes(
192 &mut μ, &mut τv, &μ_base, None, τ, ε, ins, &reg, 229 &mut μ,
193 //Some(|μ̃ : &RNDM<F, N>| calculate_residual(Pair(μ̃, &z), opA, b).norm2_squared_div2()), 230 &mut τv,
231 &μ_base,
232 τ,
233 ε,
234 ins,
235 &reg,
236 is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))),
194 ); 237 );
195 } 238 }
196 239
197 // Prune spikes with zero weight. 240 // Prune spikes with zero weight.
198 stats.pruned += prune_with_stats(&mut μ); 241 stats.pruned += prune_with_stats(&mut μ);
199 242
200 // Do z variable primal update 243 // Do z variable primal update
201 let mut z_new = τz; 244 let mut z_new = τz;
202 opKz.adjoint().gemv(&mut z_new, -σ_p, &y, -σ_p/τ); 245 opKz_adj.gemv(&mut z_new, -σ_p, &y, -σ_p / τ);
203 z_new = fnR.prox(σ_p, z_new + &z); 246 z_new = fnR.prox(σ_p, z_new + &z);
204 // Do dual update 247 // Do dual update
205 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] 248 // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}]
206 opKz.gemv(&mut y, σ_d*(1.0 + ω), &z_new, 1.0); 249 opKz.gemv(&mut y, σ_d * (1.0 + ω), &z_new, 1.0);
207 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b 250 // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b
208 opKz.gemv(&mut y, -σ_d*ω, z, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b 251 opKz.gemv(&mut y, -σ_d * ω, z, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b
209 y = starH.prox(σ_d, y); 252 y = starH.prox(σ_d, y);
210 z = z_new; 253 z = z_new;
211
212 // Update residual
213 residual = calculate_residual(Pair(&μ, &z), opA, b);
214 254
215 // Update step length parameters 255 // Update step length parameters
216 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); 256 // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ);
217 257
218 // Give statistics if requested 258 // Give statistics if requested
219 let iter = state.iteration(); 259 let iter = state.iteration();
220 stats.this_iters += 1; 260 stats.this_iters += 1;
221 261
222 state.if_verbose(|| { 262 state.if_verbose(|| {
223 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ); 263 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ);
224 full_stats(&residual, &μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) 264 full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new()))
225 }); 265 });
226 266
227 // Update main tolerance for next iteration 267 // Update main tolerance for next iteration
228 ε = tolerance.update(ε, iter); 268 ε = tolerance.update(ε, iter);
229 } 269 }
230 270
231 let fit = |μ̃ : &RNDM<F, N>| { 271 let fit = |μ̃: &RNDM<N, F>| {
232 (opA.apply(Pair(μ̃, &z))-b).norm2_squared_div2() 272 f.apply(Pair(μ̃, &z)) /*+ fnR.apply(z) + reg.apply(μ)*/
233 //+ fnR.apply(z) + reg.apply(μ)
234 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z)) 273 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z))
235 }; 274 };
236 275
237 μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v); 276 μ.merge_spikes_fitness(config.insertion.final_merging_method(), fit, |&v| v);
238 μ.prune(); 277 μ.prune();
239 Pair(μ, z) 278 Ok(Pair(μ, z))
240 } 279 }
280
281 /// Iteratively solve the pointsource localisation with an additional variable
282 /// using forward-backward splitting.
283 ///
284 /// The implementation uses [`pointsource_forward_pdps_pair`] with appropriate dummy
285 /// variables, operators, and functions.
286 #[replace_float_literals(F::cast_from(literal))]
287 pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>(
288 f: &Dat,
289 reg: &Reg,
290 prox_penalty: &P,
291 config: &FBConfig<F>,
292 iterator: I,
293 plotter: Plot,
294 (μ0, z): (Option<RNDM<N, F>>, Z),
295 //opKμ : KOpM,
296 fnR: &R,
297 ) -> DynResult<MeasureZ<F, Z, N>>
298 where
299 F: Float + ToNalgebraRealField,
300 I: AlgIteratorFactory<IterInfo<F>>,
301 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
302 S: DifferentiableRealMapping<N, F> + ClosedMul<F>,
303 RNDM<N, F>: SpikeMerging<F>,
304 Reg: RegTerm<Loc<N, F>, F>,
305 P: ProxPenalty<Loc<N, F>, S, Reg, F>,
306 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
307 Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone,
308 for<'b> &'b Z: Instance<Z>,
309 R: Prox<Z, Codomain = F>,
310 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>,
311 // We should not need to explicitly require this:
312 for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
313 {
314 let opKz = ZeroOp::new_dualisable(Loc([]), z.dual_origin());
315 let fnH = Zero::new();
316 // Convert config. We don't implement From (that could be done with the o2o crate), as σd0
317 // needs to be chosen in a general case; for the problem of this fucntion, anything is valid.
318 let &FBConfig { τ0, σp0, insertion } = config;
319 let pdps_config = ForwardPDPSConfig { τ0, σp0, insertion, σd0: 0.0 };
320
321 pointsource_forward_pdps_pair(
322 f,
323 reg,
324 prox_penalty,
325 &pdps_config,
326 iterator,
327 plotter,
328 (μ0, z, Loc([])),
329 &opKz,
330 fnR,
331 &fnH,
332 )
333 }

Mercurial Repository: pointsource_algs

mercurial