src/pdps.rs

Sun, 11 Dec 2022 23:19:17 +0200

author
Tuomo Valkonen <tuomov@iki.fi>
date
Sun, 11 Dec 2022 23:19:17 +0200
changeset 23
9869fa1e0ccd
parent 13
bdc57366d4f5
child 24
d29d1fcf5423
permissions
-rw-r--r--

Print out experiment information when running it

0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
1 /*!
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
2 Solver for the point source localisation problem with primal-dual proximal splitting.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
3
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
4 This corresponds to the manuscript
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
5
13
bdc57366d4f5 arXiv links, README beautification
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
6 * Valkonen T. - _Proximal methods for point source localisation_,
bdc57366d4f5 arXiv links, README beautification
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
7 [arXiv:2212.02991](https://arxiv.org/abs/2212.02991).
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
8
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
9 The main routine is [`pointsource_pdps`]. It is based on specilisatinn of
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
10 [`generic_pointsource_fb`] through relevant [`FBSpecialisation`] implementations.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
11 Both norm-2-squared and norm-1 data terms are supported. That is, implemented are solvers for
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
12 <div>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
13 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
14 \min_{μ ∈ ℳ(Ω)}~ F_0(Aμ - b) + α \|μ\|_{ℳ(Ω)} + δ_{≥ 0}(μ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
15 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
16 for both $F_0(y)=\frac{1}{2}\|y\|_2^2$ and $F_0(y)=\|y\|_1$ with the forward operator
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
17 $A \in 𝕃(ℳ(Ω); ℝ^n)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
18 </div>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
19
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
20 ## Approach
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
21
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
22 <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
23 The problem above can be written as
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
24 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
25 \min_μ \max_y G(μ) + ⟨y, Aμ-b⟩ - F_0^*(μ),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
26 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
27 where $G(μ) = α \|μ\|_{ℳ(Ω)} + δ_{≥ 0}(μ)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
28 The Fenchel–Rockafellar optimality conditions, employing the predual in $ℳ(Ω)$, are
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
29 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
30 0 ∈ A_*y + ∂G(μ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
31 \quad\text{and}\quad
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
32 Aμ - b ∈ ∂ F_0^*(y).
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
33 $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
34 The solution of the first part is as for forward-backward, treated in the manuscript.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
35 This is the task of <code>generic_pointsource_fb</code>, where we use <code>FBSpecialisation</code>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
36 to replace the specific residual $Aμ-b$ by $y$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
37 For $F_0(y)=\frac{1}{2}\|y\|_2^2$ the second part reads $y = Aμ -b$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
38 For $F_0(y)=\|y\|_1$ the second part reads $y ∈ ∂\|·\|_1(Aμ - b)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
39 </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
40
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
41 Based on zero initialisation for $μ$, we use the [`Subdifferentiable`] trait to make an
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
42 initialisation corresponding to the second part of the optimality conditions.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
43 In the algorithm itself, standard proximal steps are taking with respect to $F\_0^* + ⟨b, ·⟩$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
44 */
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
45
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
46 use numeric_literals::replace_float_literals;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
47 use serde::{Serialize, Deserialize};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
48 use nalgebra::DVector;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
49 use clap::ValueEnum;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
50
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
51 use alg_tools::iterate:: AlgIteratorFactory;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
52 use alg_tools::sets::Cube;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
53 use alg_tools::loc::Loc;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
54 use alg_tools::euclidean::Euclidean;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
55 use alg_tools::norms::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
56 L1, Linfinity,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
57 Projection, Norm,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
58 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
59 use alg_tools::bisection_tree::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
60 BTFN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
61 PreBTFN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
62 Bounds,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
63 BTNodeLookup,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
64 BTNode,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
65 BTSearch,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
66 P2Minimise,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
67 SupportGenerator,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
68 LocalAnalysis,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
69 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
70 use alg_tools::mapping::RealMapping;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
71 use alg_tools::nalgebra_support::ToNalgebraRealField;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
72 use alg_tools::linops::AXPY;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
73
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
74 use crate::types::*;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
75 use crate::measures::DiscreteMeasure;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
76 use crate::measures::merging::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
77 SpikeMerging,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
78 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
79 use crate::forward_model::ForwardModel;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
80 use crate::seminorms::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
81 DiscreteMeasureOp, Lipschitz
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
82 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
83 use crate::plot::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
84 SeqPlotter,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
85 Plotting,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
86 PlotLookup
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
87 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
88 use crate::fb::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
89 FBGenericConfig,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
90 FBSpecialisation,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
91 generic_pointsource_fb
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
92 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
93
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
94 /// Acceleration
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
95 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, ValueEnum, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
96 pub enum Acceleration {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
97 /// No acceleration
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
98 #[clap(name = "none")]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
99 None,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
100 /// Partial acceleration, $ω = 1/\sqrt{1+σ}$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
101 #[clap(name = "partial", help = "Partial acceleration, ω = 1/√(1+σ)")]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
102 Partial,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
103 /// Full acceleration, $ω = 1/\sqrt{1+2σ}$; no gap convergence guaranteed
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
104 #[clap(name = "full", help = "Full acceleration, ω = 1/√(1+2σ); no gap convergence guaranteed")]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
105 Full
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
106 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
107
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
108 /// Settings for [`pointsource_pdps`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
109 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
110 #[serde(default)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
111 pub struct PDPSConfig<F : Float> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
112 /// Primal step length scaling. We must have `τ0 * σ0 < 1`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
113 pub τ0 : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
114 /// Dual step length scaling. We must have `τ0 * σ0 < 1`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
115 pub σ0 : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
116 /// Accelerate if available
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
117 pub acceleration : Acceleration,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
118 /// Generic parameters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
119 pub insertion : FBGenericConfig<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
120 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
121
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
122 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
123 impl<F : Float> Default for PDPSConfig<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
124 fn default() -> Self {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
125 let τ0 = 0.5;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
126 PDPSConfig {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
127 τ0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
128 σ0 : 0.99/τ0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
129 acceleration : Acceleration::Partial,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
130 insertion : Default::default()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
131 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
132 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
133 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
134
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
135 /// Trait for subdifferentiable objects
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
136 pub trait Subdifferentiable<F : Float, V, U=V> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
137 /// Calculate some subdifferential at `x`
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
138 fn some_subdifferential(&self, x : V) -> U;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
139 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
140
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
141 /// Type for indicating norm-2-squared data fidelity.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
142 pub struct L2Squared;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
143
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
144 impl<F : Float, V : Euclidean<F>> Subdifferentiable<F, V> for L2Squared {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
145 fn some_subdifferential(&self, x : V) -> V { x }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
146 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
147
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
148 impl<F : Float + nalgebra::RealField> Subdifferentiable<F, DVector<F>> for L1 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
149 fn some_subdifferential(&self, mut x : DVector<F>) -> DVector<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
150 // nalgebra sucks for providing second copies of the same stuff that's elsewhere as well.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
151 x.iter_mut()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
152 .for_each(|v| if *v != F::ZERO { *v = *v/<F as NumTraitsFloat>::abs(*v) });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
153 x
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
154 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
155 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
156
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
157 /// Specialisation of [`generic_pointsource_fb`] to PDPS.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
158 pub struct PDPS<
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
159 'a,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
160 F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
161 A : ForwardModel<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
162 D,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
163 const N : usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
164 > {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
165 /// The data
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
166 b : &'a A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
167 /// The forward operator
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
168 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
169 /// Primal step length
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
170 τ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
171 // Dual step length
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
172 σ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
173 /// Whether acceleration should be applied (if data term supports)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
174 acceleration : Acceleration,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
175 /// The dataterm. Only used by the type system.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
176 _dataterm : D,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
177 /// Previous dual iterate.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
178 y_prev : A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
179 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
180
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
181 /// Implementation of [`FBSpecialisation`] for μPDPS with norm-2-squared data fidelity.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
182 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
183 impl<
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
184 'a,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
185 F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
186 A : ForwardModel<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
187 const N : usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
188 > FBSpecialisation<F, A::Observable, N> for PDPS<'a, F, A, L2Squared, N>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
189 where for<'b> &'b A::Observable : std::ops::Add<A::Observable, Output=A::Observable> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
190
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
191 fn update(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
192 &mut self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
193 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
194 μ_base : &DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
195 ) -> (A::Observable, Option<F>) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
196 let σ = self.σ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
197 let τ = self.τ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
198 let ω = match self.acceleration {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
199 Acceleration::None => 1.0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
200 Acceleration::Partial => {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
201 let ω = 1.0 / (1.0 + σ).sqrt();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
202 self.σ = σ * ω;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
203 self.τ = τ / ω;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
204 ω
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
205 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
206 Acceleration::Full => {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
207 let ω = 1.0 / (1.0 + 2.0 * σ).sqrt();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
208 self.σ = σ * ω;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
209 self.τ = τ / ω;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
210 ω
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
211 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
212 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
213
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
214 μ.prune();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
215
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
216 let mut y = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
217 self.opA.gemv(&mut y, 1.0 + ω, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
218 self.opA.gemv(&mut y, -ω, μ_base, 1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
219 y.axpy(1.0 / (1.0 + σ), &self.y_prev, σ / (1.0 + σ));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
220 self.y_prev.copy_from(&y);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
221
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
222 (y, Some(self.τ))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
223 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
224
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
225 fn calculate_fit(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
226 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
227 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
228 _y : &A::Observable
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
229 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
230 self.calculate_fit_simple(μ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
231 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
232
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
233 fn calculate_fit_simple(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
234 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
235 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
236 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
237 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
238 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
239 residual.norm2_squared_div2()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
240 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
241 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
242
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
243 /// Implementation of [`FBSpecialisation`] for μPDPS with norm-1 data fidelity.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
244 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
245 impl<
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
246 'a,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
247 F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
248 A : ForwardModel<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
249 const N : usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
250 > FBSpecialisation<F, A::Observable, N> for PDPS<'a, F, A, L1, N>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
251 where A::Observable : Projection<F, Linfinity> + Norm<F, L1>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
252 for<'b> &'b A::Observable : std::ops::Add<A::Observable, Output=A::Observable> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
253 fn update(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
254 &mut self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
255 μ : &mut DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
256 μ_base : &DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
257 ) -> (A::Observable, Option<F>) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
258 let σ = self.σ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
259
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
260 μ.prune();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
261
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
262 //let ȳ = self.opA.apply(μ) * 2.0 - self.opA.apply(μ_base);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
263 //*y = proj_{[-1,1]}(&self.y_prev + (ȳ - self.b) * σ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
264 let mut y = self.y_prev.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
265 self.opA.gemv(&mut y, 2.0 * σ, μ, 1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
266 self.opA.gemv(&mut y, -σ, μ_base, 1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
267 y.axpy(-σ, self.b, 1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
268 y.proj_ball_mut(1.0, Linfinity);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
269 self.y_prev.copy_from(&y);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
270
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
271 (y, None)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
272 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
273
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
274 fn calculate_fit(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
275 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
276 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
277 _y : &A::Observable
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
278 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
279 self.calculate_fit_simple(μ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
280 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
281
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
282 fn calculate_fit_simple(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
283 &self,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
284 μ : &DiscreteMeasure<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
285 ) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
286 let mut residual = self.b.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
287 self.opA.gemv(&mut residual, 1.0, μ, -1.0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
288 residual.norm(L1)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
289 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
290 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
291
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
292 /// Iteratively solve the pointsource localisation problem using primal-dual proximal splitting.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
293 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
294 /// The `dataterm` should be either [`L1`] for norm-1 data term or [`L2Squared`] for norm-2-squared.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
295 /// The settings in `config` have their [respective documentation](PDPSConfig). `opA` is the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
296 /// forward operator $A$, $b$ the observable, and $\lambda$ the regularisation weight.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
297 /// The operator `op𝒟` is used for forming the proximal term. Typically it is a convolution
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
298 /// operator. Finally, the `iterator` is an outer loop verbosity and iteration count control
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
299 /// as documented in [`alg_tools::iterate`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
300 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
301 /// For the mathematical formulation, see the [module level](self) documentation and the manuscript.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
302 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
303 /// Returns the final iterate.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
304 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
305 pub fn pointsource_pdps<'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, D, const N : usize>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
306 opA : &'a A,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
307 b : &'a A::Observable,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
308 α : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
309 op𝒟 : &'a 𝒟,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
310 config : &PDPSConfig<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
311 iterator : I,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
312 plotter : SeqPlotter<F, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
313 dataterm : D,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
314 ) -> DiscreteMeasure<Loc<F, N>, F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
315 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
316 I : AlgIteratorFactory<IterInfo<F, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
317 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
318 + std::ops::Add<A::Observable, Output=A::Observable>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
319 //+ std::ops::Mul<F, Output=A::Observable>, // <-- FIXME: compiler overflow
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
320 A::Observable : std::ops::MulAssign<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
321 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
322 A : ForwardModel<Loc<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
323 + Lipschitz<𝒟, FloatType=F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
324 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
325 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
326 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
327 𝒟::Codomain : RealMapping<F, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
328 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
329 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
330 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
331 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
332 PlotLookup : Plotting<N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
333 DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
334 PDPS<'a, F, A, D, N> : FBSpecialisation<F, A::Observable, N>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
335 D : Subdifferentiable<F, A::Observable> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
336
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
337 let y = dataterm.some_subdifferential(-b);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
338 let l = opA.lipschitz_factor(&op𝒟).unwrap().sqrt();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
339 let τ = config.τ0 / l;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
340 let σ = config.σ0 / l;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
341
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
342 let pdps = PDPS {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
343 b,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
344 opA,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
345 τ,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
346 σ,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
347 acceleration : config.acceleration,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
348 _dataterm : dataterm,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
349 y_prev : y.clone(),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
350 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
351
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
352 generic_pointsource_fb(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
353 opA, α, op𝒟, τ, &config.insertion, iterator, plotter, y,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
354 pdps
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
355 )
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
356 }

mercurial