Mercurial > repos > pointsource_algs / file comparison

comparison: src/frank_wolfe.rs

src/frank_wolfe.rs

changeset 93
8f1e57821d76
parent 63
7a8a55fd41c0
equal deleted inserted replaced
92:759634edb536 93:8f1e57821d76
18 use serde::{Deserialize, Serialize}; 18 use serde::{Deserialize, Serialize};
19 //use colored::Colorize; 19 //use colored::Colorize;
20 use crate::dataterm::QuadraticDataTerm; 20 use crate::dataterm::QuadraticDataTerm;
21 use crate::forward_model::ForwardModel; 21 use crate::forward_model::ForwardModel;
22 use crate::measures::merging::{SpikeMerging, SpikeMergingMethod}; 22 use crate::measures::merging::{SpikeMerging, SpikeMergingMethod};
23 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM}; 23 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon};
24 use crate::plot::Plotter; 24 use crate::plot::Plotter;
25 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, RegTerm}; 25 use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, RegTerm};
26 #[allow(unused_imports)] // Used in documentation 26 #[allow(unused_imports)] // Used in documentation
27 use crate::subproblem::{ 27 use crate::subproblem::{
28 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained, InnerMethod, InnerSettings, 28 nonneg::quadratic_nonneg, unconstrained::quadratic_unconstrained, InnerMethod, InnerSettings,
29 }; 29 };
30 use crate::tolerance::Tolerance; 30 use crate::tolerance::Tolerance;
31 use crate::types::*; 31 use crate::types::*;
32 use alg_tools::bisection_tree::P2Minimise;
33 use alg_tools::bounds::MinMaxMapping; 32 use alg_tools::bounds::MinMaxMapping;
34 use alg_tools::error::DynResult; 33 use alg_tools::error::DynResult;
35 use alg_tools::euclidean::Euclidean; 34 use alg_tools::euclidean::Euclidean;
36 use alg_tools::instance::Instance; 35 use alg_tools::instance::Instance;
37 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorOptions, ValueIteratorFactory}; 36 use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorOptions, ValueIteratorFactory};
38 use alg_tools::linops::Mapping; 37 use alg_tools::linops::{Mapping, Space};
39 use alg_tools::loc::Loc;
40 use alg_tools::nalgebra_support::ToNalgebraRealField; 38 use alg_tools::nalgebra_support::ToNalgebraRealField;
41 use alg_tools::norms::Norm; 39 use alg_tools::norms::Norm;
42 use alg_tools::norms::L2; 40 use alg_tools::norms::L2;
43 use alg_tools::sets::Cube;
44 41
45 /// Settings for [`pointsource_fw_reg`]. 42 /// Settings for [`pointsource_fw_reg`].
46 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 43 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
47 #[serde(default)] 44 #[serde(default)]
48 pub struct FWConfig<F: Float> { 45 pub struct FWConfig<F: Float> {
104 } 101 }
105 102
106 /// Trait for finite dimensional weight optimisation. 103 /// Trait for finite dimensional weight optimisation.
107 pub trait WeightOptim< 104 pub trait WeightOptim<
108 F: Float + ToNalgebraRealField, 105 F: Float + ToNalgebraRealField,
109 A: ForwardModel<RNDM<N, F>, F>, 106 U: Space + Clone + PartialEq,
107 A: ForwardModel<DiscreteMeasure<U, F>, F>,
110 I: AlgIteratorFactory<F>, 108 I: AlgIteratorFactory<F>,
111 const N: usize,
112 > 109 >
113 { 110 {
114 /// Return a pre-initialisation struct for [`Self::optimise_weights`]. 111 /// Return a pre-initialisation struct for [`Self::optimise_weights`].
115 /// 112 ///
116 /// The parameter `opA` is the forward operator $A$. 113 /// The parameter `opA` is the forward operator $A$.
133 /// prepared using [`Self::prepare_optimise_weights`]: 130 /// prepared using [`Self::prepare_optimise_weights`]:
134 /// 131 ///
135 /// Returns the number of iterations taken by the method configured in `inner`. 132 /// Returns the number of iterations taken by the method configured in `inner`.
136 fn optimise_weights<'a>( 133 fn optimise_weights<'a>(
137 &self, 134 &self,
138 μ: &mut RNDM<N, F>, 135 μ: &mut DiscreteMeasure<U, F>,
139 opA: &'a A, 136 opA: &'a A,
140 b: &A::Observable, 137 b: &A::Observable,
141 findim_data: &FindimData<F>, 138 findim_data: &FindimData<F>,
142 inner: &InnerSettings<F>, 139 inner: &InnerSettings<F>,
143 iterator: I, 140 iterator: I,
145 } 142 }
146 143
147 /// Trait for regularisation terms supported by [`pointsource_fw_reg`]. 144 /// Trait for regularisation terms supported by [`pointsource_fw_reg`].
148 pub trait RegTermFW< 145 pub trait RegTermFW<
149 F: Float + ToNalgebraRealField, 146 F: Float + ToNalgebraRealField,
150 A: ForwardModel<RNDM<N, F>, F>, 147 U: Space + Clone + PartialEq,
148 A: ForwardModel<DiscreteMeasure<U, F>, F>,
151 I: AlgIteratorFactory<F>, 149 I: AlgIteratorFactory<F>,
152 const N: usize, 150 >: RegTerm<U, F> + WeightOptim<F, U, A, I> + Mapping<DiscreteMeasure<U, F>, Codomain = F>
153 >: RegTerm<Loc<N, F>, F> + WeightOptim<F, A, I, N> + Mapping<RNDM<N, F>, Codomain = F>
154 { 151 {
155 /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted 152 /// With $g = A\_\*(Aμ-b)$, returns $(x, g(x))$ for $x$ a new point to be inserted
156 /// into $μ$, as determined by the regulariser. 153 /// into $μ$, as determined by the regulariser.
157 /// 154 ///
158 /// The parameters `refinement_tolerance` and `max_steps` are passed to relevant 155 /// The parameters `refinement_tolerance` and `max_steps` are passed to relevant
160 fn find_insertion( 157 fn find_insertion(
161 &self, 158 &self,
162 g: &mut A::PreadjointCodomain, 159 g: &mut A::PreadjointCodomain,
163 refinement_tolerance: F, 160 refinement_tolerance: F,
164 max_steps: usize, 161 max_steps: usize,
165 ) -> (Loc<N, F>, F); 162 ) -> (U, F);
166 163
167 /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen. 164 /// Insert point `ξ` into `μ` for the relaxed algorithm from Bredies–Pikkarainen.
168 fn relaxed_insert<'a>( 165 fn relaxed_insert<'a>(
169 &self, 166 &self,
170 μ: &mut RNDM<N, F>, 167 μ: &mut DiscreteMeasure<U, F>,
171 g: &A::PreadjointCodomain, 168 g: &A::PreadjointCodomain,
172 opA: &'a A, 169 opA: &'a A,
173 ξ: Loc<N, F>, 170 ξ: U,
174 v_ξ: F, 171 v_ξ: F,
175 findim_data: &FindimData<F>, 172 findim_data: &FindimData<F>,
176 ); 173 );
177 } 174 }
178 175
179 #[replace_float_literals(F::cast_from(literal))] 176 #[replace_float_literals(F::cast_from(literal))]
180 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N> 177 impl<F: Float + ToNalgebraRealField, U, A, I> WeightOptim<F, U, A, I> for RadonRegTerm<F>
181 for RadonRegTerm<F>
182 where 178 where
179 U: Space + Clone + PartialEq,
183 I: AlgIteratorFactory<F>, 180 I: AlgIteratorFactory<F>,
184 A: FindimQuadraticModel<Loc<N, F>, F>, 181 A: FindimQuadraticModel<U, F>,
185 { 182 {
186 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> { 183 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
187 Ok(FindimData { 184 Ok(FindimData {
188 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2), 185 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
189 m0: b.norm2_squared() / (2.0 * self.α()), 186 m0: b.norm2_squared() / (2.0 * self.α()),
190 }) 187 })
191 } 188 }
192 189
193 fn optimise_weights<'a>( 190 fn optimise_weights<'a>(
194 &self, 191 &self,
195 μ: &mut RNDM<N, F>, 192 μ: &mut DiscreteMeasure<U, F>,
196 opA: &'a A, 193 opA: &'a A,
197 b: &A::Observable, 194 b: &A::Observable,
198 findim_data: &FindimData<F>, 195 findim_data: &FindimData<F>,
199 inner: &InnerSettings<F>, 196 inner: &InnerSettings<F>,
200 iterator: I, 197 iterator: I,
219 iters 216 iters
220 } 217 }
221 } 218 }
222 219
223 #[replace_float_literals(F::cast_from(literal))] 220 #[replace_float_literals(F::cast_from(literal))]
224 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N> for RadonRegTerm<F> 221 impl<F: Float + ToNalgebraRealField, U, A, I> RegTermFW<F, U, A, I> for RadonRegTerm<F>
225 where 222 where
226 Cube<N, F>: P2Minimise<Loc<N, F>, F>, 223 U: Space + Clone + PartialEq,
224 RadonRegTerm<F>: RegTerm<U, F>,
227 I: AlgIteratorFactory<F>, 225 I: AlgIteratorFactory<F>,
228 A: FindimQuadraticModel<Loc<N, F>, F>, 226 A: FindimQuadraticModel<U, F>,
229 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>, 227 A::PreadjointCodomain: MinMaxMapping<U, F>,
230 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>, 228 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
231 // FIXME: the following *should not* be needed, they are already implied 229 // FIXME: the following *should not* be needed, they are already implied
232 RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>, 230 DiscreteMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
233 DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>, 231 DeltaMeasure<U, F>: Copy + Mapping<A::PreadjointCodomain, Codomain = F>,
234 { 232 {
235 fn find_insertion( 233 fn find_insertion(
236 &self, 234 &self,
237 g: &mut A::PreadjointCodomain, 235 g: &mut A::PreadjointCodomain,
238 refinement_tolerance: F, 236 refinement_tolerance: F,
239 max_steps: usize, 237 max_steps: usize,
240 ) -> (Loc<N, F>, F) { 238 ) -> (U, F) {
241 let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps); 239 let (ξmax, v_ξmax) = g.maximise(refinement_tolerance, max_steps);
242 let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps); 240 let (ξmin, v_ξmin) = g.minimise(refinement_tolerance, max_steps);
243 if v_ξmin < 0.0 && -v_ξmin > v_ξmax { 241 if v_ξmin < 0.0 && -v_ξmin > v_ξmax {
244 (ξmin, v_ξmin) 242 (ξmin, v_ξmin)
245 } else { 243 } else {
247 } 245 }
248 } 246 }
249 247
250 fn relaxed_insert<'a>( 248 fn relaxed_insert<'a>(
251 &self, 249 &self,
252 μ: &mut RNDM<N, F>, 250 μ: &mut DiscreteMeasure<U, F>,
253 g: &A::PreadjointCodomain, 251 g: &A::PreadjointCodomain,
254 opA: &'a A, 252 opA: &'a A,
255 ξ: Loc<N, F>, 253 ξ: U,
256 v_ξ: F, 254 v_ξ: F,
257 findim_data: &FindimData<F>, 255 findim_data: &FindimData<F>,
258 ) { 256 ) {
259 let α = self.0; 257 let α = self.0;
260 let m0 = findim_data.m0; 258 let m0 = findim_data.m0;
283 *μ += δ * s; 281 *μ += δ * s;
284 } 282 }
285 } 283 }
286 284
287 #[replace_float_literals(F::cast_from(literal))] 285 #[replace_float_literals(F::cast_from(literal))]
288 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> WeightOptim<F, A, I, N> 286 impl<F: Float + ToNalgebraRealField, A, U, I> WeightOptim<F, U, A, I> for NonnegRadonRegTerm<F>
289 for NonnegRadonRegTerm<F>
290 where 287 where
288 U: Space + Clone + PartialEq,
291 I: AlgIteratorFactory<F>, 289 I: AlgIteratorFactory<F>,
292 A: FindimQuadraticModel<Loc<N, F>, F>, 290 A: FindimQuadraticModel<U, F>,
293 { 291 {
294 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> { 292 fn prepare_optimise_weights(&self, opA: &A, b: &A::Observable) -> DynResult<FindimData<F>> {
295 Ok(FindimData { 293 Ok(FindimData {
296 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2), 294 opAnorm_squared: opA.opnorm_bound(Radon, L2)?.powi(2),
297 m0: b.norm2_squared() / (2.0 * self.α()), 295 m0: b.norm2_squared() / (2.0 * self.α()),
298 }) 296 })
299 } 297 }
300 298
301 fn optimise_weights<'a>( 299 fn optimise_weights<'a>(
302 &self, 300 &self,
303 μ: &mut RNDM<N, F>, 301 μ: &mut DiscreteMeasure<U, F>,
304 opA: &'a A, 302 opA: &'a A,
305 b: &A::Observable, 303 b: &A::Observable,
306 findim_data: &FindimData<F>, 304 findim_data: &FindimData<F>,
307 inner: &InnerSettings<F>, 305 inner: &InnerSettings<F>,
308 iterator: I, 306 iterator: I,
327 iters 325 iters
328 } 326 }
329 } 327 }
330 328
331 #[replace_float_literals(F::cast_from(literal))] 329 #[replace_float_literals(F::cast_from(literal))]
332 impl<F: Float + ToNalgebraRealField, A, I, const N: usize> RegTermFW<F, A, I, N> 330 impl<F: Float + ToNalgebraRealField, U, A, I> RegTermFW<F, U, A, I> for NonnegRadonRegTerm<F>
333 for NonnegRadonRegTerm<F>
334 where 331 where
335 Cube<N, F>: P2Minimise<Loc<N, F>, F>, 332 U: Space + Clone + PartialEq,
333 NonnegRadonRegTerm<F>: RegTerm<U, F>,
336 I: AlgIteratorFactory<F>, 334 I: AlgIteratorFactory<F>,
337 A: FindimQuadraticModel<Loc<N, F>, F>, 335 A: FindimQuadraticModel<U, F>,
338 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>, 336 A::PreadjointCodomain: MinMaxMapping<U, F>,
339 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>, 337 for<'a> &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
340 // FIXME: the following *should not* be needed, they are already implied 338 // FIXME: the following *should not* be needed, they are already implied
341 RNDM<N, F>: Mapping<A::PreadjointCodomain, Codomain = F>, 339 DiscreteMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
342 DeltaMeasure<Loc<N, F>, F>: Mapping<A::PreadjointCodomain, Codomain = F>, 340 DeltaMeasure<U, F>: Mapping<A::PreadjointCodomain, Codomain = F>,
343 { 341 {
344 fn find_insertion( 342 fn find_insertion(
345 &self, 343 &self,
346 g: &mut A::PreadjointCodomain, 344 g: &mut A::PreadjointCodomain,
347 refinement_tolerance: F, 345 refinement_tolerance: F,
348 max_steps: usize, 346 max_steps: usize,
349 ) -> (Loc<N, F>, F) { 347 ) -> (U, F) {
350 g.maximise(refinement_tolerance, max_steps) 348 g.maximise(refinement_tolerance, max_steps)
351 } 349 }
352 350
353 fn relaxed_insert<'a>( 351 fn relaxed_insert<'a>(
354 &self, 352 &self,
355 μ: &mut RNDM<N, F>, 353 μ: &mut DiscreteMeasure<U, F>,
356 g: &A::PreadjointCodomain, 354 g: &A::PreadjointCodomain,
357 opA: &'a A, 355 opA: &'a A,
358 ξ: Loc<N, F>, 356 ξ: U,
359 v_ξ: F, 357 v_ξ: F,
360 findim_data: &FindimData<F>, 358 findim_data: &FindimData<F>,
361 ) { 359 ) {
362 // This is just a verbatim copy of RadonRegTerm::relaxed_insert. 360 // This is just a verbatim copy of RadonRegTerm::relaxed_insert.
363 let α = self.0; 361 let α = self.0;
399 /// The `opA` parameter is the forward operator $A$, while `b`$ is as in the 397 /// The `opA` parameter is the forward operator $A$, while `b`$ is as in the
400 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while 398 /// objective above. The method parameter are set in `config` (see [`FWConfig`]), while
401 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to 399 /// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
402 /// save intermediate iteration states as images. 400 /// save intermediate iteration states as images.
403 #[replace_float_literals(F::cast_from(literal))] 401 #[replace_float_literals(F::cast_from(literal))]
404 pub fn pointsource_fw_reg<'a, F, I, A, Reg, Plot, const N: usize>( 402 pub fn pointsource_fw_reg<'a, F, U, I, A, Reg, Plot>(
405 f: &'a QuadraticDataTerm<F, RNDM<N, F>, A>, 403 f: &'a QuadraticDataTerm<F, DiscreteMeasure<U, F>, A>,
406 reg: &Reg, 404 reg: &Reg,
407 //domain : Cube<N, F>,
408 config: &FWConfig<F>, 405 config: &FWConfig<F>,
409 iterator: I, 406 iterator: I,
410 mut plotter: Plot, 407 mut plotter: Plot,
411 μ0: Option<RNDM<N, F>>, 408 μ0: Option<DiscreteMeasure<U, F>>,
412 ) -> DynResult<RNDM<N, F>> 409 ) -> DynResult<DiscreteMeasure<U, F>>
413 where 410 where
414 F: Float + ToNalgebraRealField, 411 F: Float + ToNalgebraRealField,
412 U: Space + Clone + PartialEq,
415 I: AlgIteratorFactory<IterInfo<F>>, 413 I: AlgIteratorFactory<IterInfo<F>>,
416 A: ForwardModel<RNDM<N, F>, F>, 414 A: ForwardModel<DiscreteMeasure<U, F>, F>,
417 A::PreadjointCodomain: MinMaxMapping<Loc<N, F>, F>, 415 A::PreadjointCodomain: MinMaxMapping<U, F>,
418 &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>, 416 &'a A::PreadjointCodomain: Instance<A::PreadjointCodomain>,
419 Cube<N, F>: P2Minimise<Loc<N, F>, F>, 417 //Cube<N, F>: P2Minimise<U, F>,
420 RNDM<N, F>: SpikeMerging<F>, 418 DiscreteMeasure<U, F>: SpikeMerging<F>,
421 Reg: RegTermFW<F, A, ValueIteratorFactory<F, AlgIteratorOptions>, N>, 419 Reg: RegTermFW<F, U, A, ValueIteratorFactory<F, AlgIteratorOptions>>,
422 Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, RNDM<N, F>>, 420 Plot: Plotter<A::PreadjointCodomain, A::PreadjointCodomain, DiscreteMeasure<U, F>>,
423 { 421 {
424 let opA = f.operator(); 422 let opA = f.operator();
425 let b = f.data(); 423 let b = f.data();
426 424
427 // Set up parameters 425 // Set up parameters
436 // Initialise iterates 434 // Initialise iterates
437 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); 435 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new());
438 let mut residual = f.residual(&μ); 436 let mut residual = f.residual(&μ);
439 437
440 // Statistics 438 // Statistics
441 let full_stats = |residual: &A::Observable, ν: &RNDM<N, F>, ε, stats| IterInfo { 439 let full_stats = |residual: &A::Observable, ν: &DiscreteMeasure<U, F>, ε, stats| IterInfo {
442 value: residual.norm2_squared_div2() + reg.apply(ν), 440 value: residual.norm2_squared_div2() + reg.apply(ν),
443 n_spikes: ν.len(), 441 n_spikes: ν.len(),
444 ε, 442 ε,
445 ..stats 443 ..stats
446 }; 444 };

Mercurial Repository: pointsource_algs

mercurial