Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_pdps.rs

src/forward_pdps.rs

changeset 93
8f1e57821d76
parent 91
47d495348e38
child 95
d9318317ccac
equal deleted inserted replaced
92:759634edb536 93:8f1e57821d76
4 */ 4 */
5 5
6 use crate::domain::{AdmissibleDomain, FullSpace}; 6 use crate::domain::{AdmissibleDomain, FullSpace};
7 use crate::fb::*; 7 use crate::fb::*;
8 use crate::measures::merging::SpikeMerging; 8 use crate::measures::merging::SpikeMerging;
9 use crate::measures::{DiscreteMeasure, RNDM}; 9 use crate::measures::DiscreteMeasure;
10 use crate::plot::Plotter; 10 use crate::plot::Plotter;
11 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair}; 11 use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair};
12 use crate::regularisation::RegTerm; 12 use crate::regularisation::RegTerm;
13 use crate::types::*; 13 use crate::types::*;
14 use alg_tools::convex::{Conjugable, Prox, Zero}; 14 use alg_tools::convex::{Conjugable, Prox, Zero};
15 use alg_tools::direct_product::Pair; 15 use alg_tools::direct_product::Pair;
16 use alg_tools::error::DynResult; 16 use alg_tools::error::DynResult;
17 use alg_tools::euclidean::ClosedEuclidean; 17 use alg_tools::euclidean::ClosedEuclidean;
18 use alg_tools::iterate::AlgIteratorFactory; 18 use alg_tools::iterate::AlgIteratorFactory;
19 use alg_tools::linops::{BoundedLinear, IdOp, SimplyAdjointable, ZeroOp, AXPY, GEMV}; 19 use alg_tools::linops::{BoundedLinear, IdOp, SimplyAdjointable, Space, ZeroOp, AXPY, GEMV};
20 use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Instance}; 20 use alg_tools::mapping::{DifferentiableMapping, Instance};
21 use alg_tools::nalgebra_support::ToNalgebraRealField; 21 use alg_tools::nalgebra_support::ToNalgebraRealField;
22 use alg_tools::norms::L2; 22 use alg_tools::norms::L2;
23 use anyhow::ensure; 23 use anyhow::ensure;
24 use measures::DeltaMeasure; 24 use measures::DeltaMeasure;
25 use numeric_literals::replace_float_literals; 25 use numeric_literals::replace_float_literals;
46 fn default() -> Self { 46 fn default() -> Self {
47 ForwardPDPSConfig { τ0: 0.99, σd0: 0.05, σp0: 0.99, insertion: Default::default() } 47 ForwardPDPSConfig { τ0: 0.99, σd0: 0.05, σp0: 0.99, insertion: Default::default() }
48 } 48 }
49 } 49 }
50 50
51 type MeasureZ<F, Z, const N: usize> = Pair<RNDM<N, F>, Z>; 51 type MeasureZ<F, Z, U> = Pair<DiscreteMeasure<U, F>, Z>;
52 52
53 /// Iteratively solve the pointsource localisation with an additional variable 53 /// Iteratively solve the pointsource localisation with an additional variable
54 /// using primal-dual proximal splitting with a forward step. 54 /// using primal-dual proximal splitting with a forward step.
55 /// 55 ///
56 /// The problem is 56 /// The problem is
87 /// τ = τ_0 \frac{1 - σ_{p,0}}{(σ_d M (1-σ_p L_z) + (1 - σ_{p,0} L)}. 87 /// τ = τ_0 \frac{1 - σ_{p,0}}{(σ_d M (1-σ_p L_z) + (1 - σ_{p,0} L)}.
88 /// $$ 88 /// $$
89 #[replace_float_literals(F::cast_from(literal))] 89 #[replace_float_literals(F::cast_from(literal))]
90 pub fn pointsource_forward_pdps_pair_dom< 90 pub fn pointsource_forward_pdps_pair_dom<
91 F, 91 F,
92 U,
92 I, 93 I,
93 S, 94 S,
94 Dat, 95 Dat,
95 Reg, 96 Reg,
96 P, 97 P,
98 R, 99 R,
99 Y, 100 Y,
100 /*KOpM, */ KOpZ, 101 /*KOpM, */ KOpZ,
101 H, 102 H,
102 Plot, 103 Plot,
103 const N: usize,
104 >( 104 >(
105 f: &Dat, 105 f: &Dat,
106 reg: &Reg, 106 reg: &Reg,
107 prox_penalty: &P, 107 prox_penalty: &P,
108 config: &ForwardPDPSConfig<F>, 108 config: &ForwardPDPSConfig<F>,
109 iterator: I, 109 iterator: I,
110 mut plotter: Plot, 110 mut plotter: Plot,
111 (μ0, mut z, mut y): (Option<RNDM<N, F>>, Z, Y), 111 (μ0, mut z, mut y): (Option<DiscreteMeasure<U, F>>, Z, Y),
112 //opKμ : KOpM, 112 //opKμ : KOpM,
113 opKz: &KOpZ, 113 opKz: &KOpZ,
114 fnR: &R, 114 fnR: &R,
115 fnH: &H, 115 fnH: &H,
116 domain: &mut impl AdmissibleDomain<Loc<N, F>>, 116 domain: &mut impl AdmissibleDomain<U>,
117 ) -> DynResult<MeasureZ<F, Z, N>> 117 ) -> DynResult<MeasureZ<F, Z, U>>
118 where 118 where
119 U: Space + Clone,
119 F: Float + ToNalgebraRealField, 120 F: Float + ToNalgebraRealField,
120 I: AlgIteratorFactory<IterInfo<F>>, 121 I: AlgIteratorFactory<IterInfo<F>>,
121 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>, 122 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
122 //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex… 123 //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex…
123 S: DifferentiableRealMapping<N, F> + ClosedMul<F>, 124 S: ClosedMul<F>,
124 RNDM<N, F>: SpikeMerging<F>, 125 DiscreteMeasure<U, F>: SpikeMerging<F>,
125 Reg: RegTerm<Loc<N, F>, F>, 126 Reg: RegTerm<U, F>,
126 P: ProxPenalty<Loc<N, F>, S, Reg, F>, 127 P: ProxPenalty<U, S, Reg, F>,
127 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, 128 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
128 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y> 129 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
129 + GEMV<F, Z, Y> 130 + GEMV<F, Z, Y>
130 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>, 131 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>,
131 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>, 132 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>,
134 Z: ClosedEuclidean<F> + AXPY, 135 Z: ClosedEuclidean<F> + AXPY,
135 for<'b> &'b Z: Instance<Z>, 136 for<'b> &'b Z: Instance<Z>,
136 R: Prox<Z, Codomain = F>, 137 R: Prox<Z, Codomain = F>,
137 H: Conjugable<Y, F, Codomain = F>, 138 H: Conjugable<Y, F, Codomain = F>,
138 for<'b> H::Conjugate<'b>: Prox<Y>, 139 for<'b> H::Conjugate<'b>: Prox<Y>,
139 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, 140 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
140 { 141 {
141 // Check parameters 142 // Check parameters
142 ensure!( 143 ensure!(
143 config.τ0 > 0.0 144 config.τ0 > 0.0
144 && config.τ0 < 1.0 145 && config.τ0 < 1.0
186 let mut ε = tolerance.initial(); 187 let mut ε = tolerance.initial();
187 188
188 let starH = fnH.conjugate(); 189 let starH = fnH.conjugate();
189 190
190 // Statistics 191 // Statistics
191 let full_stats = |μ: &RNDM<N, F>, z: &Z, ε, stats| IterInfo { 192 let full_stats = |μ: &DiscreteMeasure<U, F>, z: &Z, ε, stats| IterInfo {
192 value: f.apply(Pair(μ, z)) 193 value: f.apply(Pair(μ, z))
193 + fnR.apply(z) 194 + fnR.apply(z)
194 + reg.apply(μ) 195 + reg.apply(μ)
195 + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)), 196 + fnH.apply(/* opKμ.apply(μ) + */ opKz.apply(z)),
196 n_spikes: μ.len(), 197 n_spikes: μ.len(),
248 &μ_base, 249 &μ_base,
249 τ, 250 τ,
250 ε, 251 ε,
251 ins, 252 ins,
252 &reg, 253 &reg,
253 Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)), 254 Some(|μ̃: &DiscreteMeasure<U, F>| f.apply(Pair(μ̃, &z)) + reg.apply(μ̃)),
254 ); 255 );
255 // Prune. Do not prune new insertions on first round. 256 // Prune. Do not prune new insertions on first round.
256 stats.pruned += prune_with_stats(&mut μ, μ_base_len); 257 stats.pruned += prune_with_stats(&mut μ, μ_base_len);
257 258
258 if config.insertion.random_insertions { 259 if config.insertion.random_insertions {
308 309
309 // Update main tolerance for next iteration 310 // Update main tolerance for next iteration
310 ε = tolerance.update(ε, iter); 311 ε = tolerance.update(ε, iter);
311 } 312 }
312 313
313 let fit = |μ̃: &RNDM<N, F>| { 314 let fit = |μ̃: &DiscreteMeasure<U, F>| {
314 f.apply(Pair(μ̃, &z)) /*+ fnR.apply(z) + reg.apply(μ)*/ 315 f.apply(Pair(μ̃, &z)) /*+ fnR.apply(z) + reg.apply(μ)*/
315 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z)) 316 + fnH.apply(/* opKμ.apply(&μ̃) + */ opKz.apply(&z))
316 + reg.apply(μ̃) 317 + reg.apply(μ̃)
317 }; 318 };
318 319
322 } 323 }
323 324
324 #[deprecated] 325 #[deprecated]
325 pub fn pointsource_forward_pdps_pair< 326 pub fn pointsource_forward_pdps_pair<
326 F, 327 F,
328 U,
327 I, 329 I,
328 S, 330 S,
329 Dat, 331 Dat,
330 Reg, 332 Reg,
331 P, 333 P,
333 R, 335 R,
334 Y, 336 Y,
335 /*KOpM, */ KOpZ, 337 /*KOpM, */ KOpZ,
336 H, 338 H,
337 Plot, 339 Plot,
338 const N: usize,
339 >( 340 >(
340 f: &Dat, 341 f: &Dat,
341 reg: &Reg, 342 reg: &Reg,
342 prox_penalty: &P, 343 prox_penalty: &P,
343 config: &ForwardPDPSConfig<F>, 344 config: &ForwardPDPSConfig<F>,
344 iterator: I, 345 iterator: I,
345 plotter: Plot, 346 plotter: Plot,
346 init: (Option<RNDM<N, F>>, Z, Y), 347 init: (Option<DiscreteMeasure<U, F>>, Z, Y),
347 //opKμ : KOpM, 348 //opKμ : KOpM,
348 opKz: &KOpZ, 349 opKz: &KOpZ,
349 fnR: &R, 350 fnR: &R,
350 fnH: &H, 351 fnH: &H,
351 ) -> DynResult<MeasureZ<F, Z, N>> 352 ) -> DynResult<MeasureZ<F, Z, U>>
352 where 353 where
354 U: Space + Clone,
353 F: Float + ToNalgebraRealField, 355 F: Float + ToNalgebraRealField,
354 I: AlgIteratorFactory<IterInfo<F>>, 356 I: AlgIteratorFactory<IterInfo<F>>,
355 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>, 357 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
356 //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex… 358 //Pair<S, Z>: ClosedMul<F>, // Doesn't really need to be closed, if make this signature more complex…
357 S: DifferentiableRealMapping<N, F> + ClosedMul<F>, 359 S: ClosedMul<F>,
358 RNDM<N, F>: SpikeMerging<F>, 360 DiscreteMeasure<U, F>: SpikeMerging<F>,
359 Reg: RegTerm<Loc<N, F>, F>, 361 Reg: RegTerm<U, F>,
360 P: ProxPenalty<Loc<N, F>, S, Reg, F>, 362 P: ProxPenalty<U, S, Reg, F>,
361 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, 363 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
362 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y> 364 KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y>
363 + GEMV<F, Z, Y> 365 + GEMV<F, Z, Y>
364 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>, 366 + SimplyAdjointable<Z, Y, Codomain = Y, AdjointCodomain = Z>,
365 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>, 367 KOpZ::SimpleAdjoint: GEMV<F, Y, Z>,
368 Z: ClosedEuclidean<F> + AXPY, 370 Z: ClosedEuclidean<F> + AXPY,
369 for<'b> &'b Z: Instance<Z>, 371 for<'b> &'b Z: Instance<Z>,
370 R: Prox<Z, Codomain = F>, 372 R: Prox<Z, Codomain = F>,
371 H: Conjugable<Y, F, Codomain = F>, 373 H: Conjugable<Y, F, Codomain = F>,
372 for<'b> H::Conjugate<'b>: Prox<Y>, 374 for<'b> H::Conjugate<'b>: Prox<Y>,
373 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, 375 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
374 { 376 {
375 pointsource_forward_pdps_pair_dom( 377 pointsource_forward_pdps_pair_dom(
376 f, 378 f,
377 reg, 379 reg,
378 prox_penalty, 380 prox_penalty,
391 /// using forward-backward splitting. 393 /// using forward-backward splitting.
392 /// 394 ///
393 /// The implementation uses [`pointsource_forward_pdps_pair`] with appropriate dummy 395 /// The implementation uses [`pointsource_forward_pdps_pair`] with appropriate dummy
394 /// variables, operators, and functions. 396 /// variables, operators, and functions.
395 #[replace_float_literals(F::cast_from(literal))] 397 #[replace_float_literals(F::cast_from(literal))]
396 pub fn pointsource_fb_pair_dom<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>( 398 pub fn pointsource_fb_pair_dom<F, U, I, S, Dat, Reg, P, Z, R, Plot>(
397 f: &Dat, 399 f: &Dat,
398 reg: &Reg, 400 reg: &Reg,
399 prox_penalty: &P, 401 prox_penalty: &P,
400 config: &FBConfig<F>, 402 config: &FBConfig<F>,
401 iterator: I, 403 iterator: I,
402 plotter: Plot, 404 plotter: Plot,
403 (μ0, z): (Option<RNDM<N, F>>, Z), 405 (μ0, z): (Option<DiscreteMeasure<U, F>>, Z),
404 //opKμ : KOpM, 406 //opKμ : KOpM,
405 fnR: &R, 407 fnR: &R,
406 domain: &mut impl AdmissibleDomain<Loc<N, F>>, 408 domain: &mut impl AdmissibleDomain<U>,
407 ) -> DynResult<MeasureZ<F, Z, N>> 409 ) -> DynResult<MeasureZ<F, Z, U>>
408 where 410 where
411 U: Space + Clone,
409 F: Float + ToNalgebraRealField, 412 F: Float + ToNalgebraRealField,
410 I: AlgIteratorFactory<IterInfo<F>>, 413 I: AlgIteratorFactory<IterInfo<F>>,
411 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>, 414 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
412 S: DifferentiableRealMapping<N, F> + ClosedMul<F>, 415 S: ClosedMul<F>,
413 RNDM<N, F>: SpikeMerging<F>, 416 DiscreteMeasure<U, F>: SpikeMerging<F>,
414 Reg: RegTerm<Loc<N, F>, F>, 417 Reg: RegTerm<U, F>,
415 P: ProxPenalty<Loc<N, F>, S, Reg, F>, 418 P: ProxPenalty<U, S, Reg, F>,
416 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, 419 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
417 Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone, 420 Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone,
418 for<'b> &'b Z: Instance<Z>, 421 for<'b> &'b Z: Instance<Z>,
419 R: Prox<Z, Codomain = F>, 422 R: Prox<Z, Codomain = F>,
420 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, 423 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
421 // We should not need to explicitly require this: 424 // We should not need to explicitly require this:
422 for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>, 425 for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
423 { 426 {
424 let opKz = ZeroOp::new_dualisable(Loc([]), z.dual_origin()); 427 let opKz = ZeroOp::new_dualisable(Loc([]), z.dual_origin());
425 let fnH = Zero::new(); 428 let fnH = Zero::new();
443 ) 446 )
444 } 447 }
445 448
446 #[deprecated] 449 #[deprecated]
447 #[replace_float_literals(F::cast_from(literal))] 450 #[replace_float_literals(F::cast_from(literal))]
448 pub fn pointsource_fb_pair<F, I, S, Dat, Reg, P, Z, R, Plot, const N: usize>( 451 pub fn pointsource_fb_pair<F, U, I, S, Dat, Reg, P, Z, R, Plot>(
449 f: &Dat, 452 f: &Dat,
450 reg: &Reg, 453 reg: &Reg,
451 prox_penalty: &P, 454 prox_penalty: &P,
452 config: &FBConfig<F>, 455 config: &FBConfig<F>,
453 iterator: I, 456 iterator: I,
454 plotter: Plot, 457 plotter: Plot,
455 init: (Option<RNDM<N, F>>, Z), 458 init: (Option<DiscreteMeasure<U, F>>, Z),
456 //opKμ : KOpM, 459 //opKμ : KOpM,
457 fnR: &R, 460 fnR: &R,
458 ) -> DynResult<MeasureZ<F, Z, N>> 461 ) -> DynResult<MeasureZ<F, Z, U>>
459 where 462 where
463 U: Space + Clone,
460 F: Float + ToNalgebraRealField, 464 F: Float + ToNalgebraRealField,
461 I: AlgIteratorFactory<IterInfo<F>>, 465 I: AlgIteratorFactory<IterInfo<F>>,
462 Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>>, 466 Dat: DifferentiableMapping<MeasureZ<F, Z, U>, Codomain = F, DerivativeDomain = Pair<S, Z>>,
463 S: DifferentiableRealMapping<N, F> + ClosedMul<F>, 467 S: DifferentiableMapping<U, Codomain = F> + ClosedMul<F>,
464 RNDM<N, F>: SpikeMerging<F>, 468 DiscreteMeasure<U, F>: SpikeMerging<F>,
465 Reg: RegTerm<Loc<N, F>, F>, 469 Reg: RegTerm<U, F>,
466 P: ProxPenalty<Loc<N, F>, S, Reg, F>, 470 P: ProxPenalty<U, S, Reg, F>,
467 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, 471 for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>,
468 Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone, 472 Z: ClosedEuclidean<F> + AXPY<Field = F> + Clone,
469 for<'b> &'b Z: Instance<Z>, 473 for<'b> &'b Z: Instance<Z>,
470 R: Prox<Z, Codomain = F>, 474 R: Prox<Z, Codomain = F>,
471 Plot: Plotter<P::ReturnMapping, S, RNDM<N, F>>, 475 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>,
472 // We should not need to explicitly require this: 476 // We should not need to explicitly require this:
473 for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>, 477 for<'b> &'b Loc<0, F>: Instance<Loc<0, F>>,
474 { 478 {
475 pointsource_fb_pair_dom( 479 pointsource_fb_pair_dom(
476 f, 480 f,

Mercurial Repository: pointsource_algs

mercurial