Mercurial > repos > pointsource_algs / file comparison

comparison: src/fb.rs

src/fb.rs

branch
dev
changeset 37
c5d8bd1a7728
parent 35
b087e3eab191
child 39
6316d68b58af
equal deleted inserted replaced
36:fb911f72e698 37:c5d8bd1a7728
78 */ 78 */
79 79
80 use numeric_literals::replace_float_literals; 80 use numeric_literals::replace_float_literals;
81 use serde::{Serialize, Deserialize}; 81 use serde::{Serialize, Deserialize};
82 use colored::Colorize; 82 use colored::Colorize;
83 use nalgebra::DVector; 83
84 84 use alg_tools::iterate::AlgIteratorFactory;
85 use alg_tools::iterate::{
86 AlgIteratorFactory,
87 AlgIteratorIteration,
88 AlgIterator,
89 };
90 use alg_tools::euclidean::Euclidean; 85 use alg_tools::euclidean::Euclidean;
91 use alg_tools::linops::{Mapping, GEMV}; 86 use alg_tools::linops::{Mapping, GEMV};
92 use alg_tools::sets::Cube;
93 use alg_tools::loc::Loc;
94 use alg_tools::bisection_tree::{
95 BTFN,
96 PreBTFN,
97 Bounds,
98 BTNodeLookup,
99 BTNode,
100 BTSearch,
101 P2Minimise,
102 SupportGenerator,
103 LocalAnalysis,
104 BothGenerators,
105 };
106 use alg_tools::mapping::RealMapping; 87 use alg_tools::mapping::RealMapping;
107 use alg_tools::nalgebra_support::ToNalgebraRealField; 88 use alg_tools::nalgebra_support::ToNalgebraRealField;
108 use alg_tools::instance::Instance; 89 use alg_tools::instance::Instance;
109 use alg_tools::norms::Linfinity;
110 90
111 use crate::types::*; 91 use crate::types::*;
112 use crate::measures::{ 92 use crate::measures::{
113 DiscreteMeasure, 93 DiscreteMeasure,
114 RNDM, 94 RNDM,
115 DeltaMeasure,
116 Radon,
117 }; 95 };
118 use crate::measures::merging::{ 96 use crate::measures::merging::{
119 SpikeMergingMethod, 97 SpikeMergingMethod,
120 SpikeMerging, 98 SpikeMerging,
121 }; 99 };
122 use crate::forward_model::{ 100 use crate::forward_model::{
123 ForwardModel, 101 ForwardModel,
124 AdjointProductBoundedBy 102 AdjointProductBoundedBy,
125 }; 103 };
126 use crate::seminorms::DiscreteMeasureOp;
127 use crate::subproblem::{
128 InnerSettings,
129 InnerMethod,
130 };
131 use crate::tolerance::Tolerance;
132 use crate::plot::{ 104 use crate::plot::{
133 SeqPlotter, 105 SeqPlotter,
134 Plotting, 106 Plotting,
135 PlotLookup 107 PlotLookup
136 }; 108 };
137 use crate::regularisation::RegTerm; 109 use crate::regularisation::RegTerm;
138 use crate::dataterm::{ 110 use crate::dataterm::{
139 calculate_residual, 111 calculate_residual,
140 L2Squared, 112 L2Squared,
141 DataTerm, 113 DataTerm,
114 };
115 pub use crate::prox_penalty::{
116 FBGenericConfig,
117 ProxPenalty
142 }; 118 };
143 119
144 /// Settings for [`pointsource_fb_reg`]. 120 /// Settings for [`pointsource_fb_reg`].
145 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] 121 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
146 #[serde(default)] 122 #[serde(default)]
149 pub τ0 : F, 125 pub τ0 : F,
150 /// Generic parameters 126 /// Generic parameters
151 pub generic : FBGenericConfig<F>, 127 pub generic : FBGenericConfig<F>,
152 } 128 }
153 129
154 /// Settings for the solution of the stepwise optimality condition.
155 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
156 #[serde(default)]
157 pub struct FBGenericConfig<F : Float> {
158 /// Tolerance for point insertion.
159 pub tolerance : Tolerance<F>,
160
161 /// Stop looking for predual maximum (where to isert a new point) below
162 /// `tolerance` multiplied by this factor.
163 ///
164 /// Not used by [`super::radon_fb`].
165 pub insertion_cutoff_factor : F,
166
167 /// Settings for branch and bound refinement when looking for predual maxima
168 pub refinement : RefinementSettings<F>,
169
170 /// Maximum insertions within each outer iteration
171 ///
172 /// Not used by [`super::radon_fb`].
173 pub max_insertions : usize,
174
175 /// Pair `(n, m)` for maximum insertions `m` on first `n` iterations.
176 ///
177 /// Not used by [`super::radon_fb`].
178 pub bootstrap_insertions : Option<(usize, usize)>,
179
180 /// Inner method settings
181 pub inner : InnerSettings<F>,
182
183 /// Spike merging method
184 pub merging : SpikeMergingMethod<F>,
185
186 /// Tolerance multiplier for merges
187 pub merge_tolerance_mult : F,
188
189 /// Spike merging method after the last step
190 pub final_merging : SpikeMergingMethod<F>,
191
192 /// Iterations between merging heuristic tries
193 pub merge_every : usize,
194
195 // /// Save $μ$ for postprocessing optimisation
196 // pub postprocessing : bool
197 }
198
199 #[replace_float_literals(F::cast_from(literal))] 130 #[replace_float_literals(F::cast_from(literal))]
200 impl<F : Float> Default for FBConfig<F> { 131 impl<F : Float> Default for FBConfig<F> {
201 fn default() -> Self { 132 fn default() -> Self {
202 FBConfig { 133 FBConfig {
203 τ0 : 0.99, 134 τ0 : 0.99,
204 generic : Default::default(), 135 generic : Default::default(),
205 } 136 }
206 } 137 }
207 }
208
209 #[replace_float_literals(F::cast_from(literal))]
210 impl<F : Float> Default for FBGenericConfig<F> {
211 fn default() -> Self {
212 FBGenericConfig {
213 tolerance : Default::default(),
214 insertion_cutoff_factor : 1.0,
215 refinement : Default::default(),
216 max_insertions : 100,
217 //bootstrap_insertions : None,
218 bootstrap_insertions : Some((10, 1)),
219 inner : InnerSettings {
220 method : InnerMethod::Default,
221 .. Default::default()
222 },
223 merging : SpikeMergingMethod::None,
224 //merging : Default::default(),
225 final_merging : Default::default(),
226 merge_every : 10,
227 merge_tolerance_mult : 2.0,
228 // postprocessing : false,
229 }
230 }
231 }
232
233 impl<F : Float> FBGenericConfig<F> {
234 /// Check if merging should be attempted this iteration
235 pub fn merge_now<I : AlgIterator>(&self, state : &AlgIteratorIteration<I>) -> bool {
236 state.iteration() % self.merge_every == 0
237 }
238 }
239
240 /// TODO: document.
241 /// `μ_base + ν_delta` is the base point, where `μ` and `μ_base` are assumed to have the same spike
242 /// locations, while `ν_delta` may have different locations.
243 #[replace_float_literals(F::cast_from(literal))]
244 pub(crate) fn insert_and_reweigh<
245 'a, F, GA, 𝒟, BTA, G𝒟, S, K, Reg, I, const N : usize
246 >(
247 μ : &mut RNDM<F, N>,
248 τv : &BTFN<F, GA, BTA, N>,
249 μ_base : &RNDM<F, N>,
250 ν_delta: Option<&RNDM<F, N>>,
251 op𝒟 : &'a 𝒟,
252 op𝒟norm : F,
253 τ : F,
254 ε : F,
255 config : &FBGenericConfig<F>,
256 reg : &Reg,
257 state : &AlgIteratorIteration<I>,
258 stats : &mut IterInfo<F, N>,
259 ) -> (BTFN<F, BothGenerators<GA, G𝒟>, BTA, N>, bool)
260 where F : Float + ToNalgebraRealField,
261 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone,
262 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>,
263 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone,
264 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>,
265 𝒟::Codomain : RealMapping<F, N>,
266 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
267 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>,
268 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
269 Reg : RegTerm<F, N>,
270 I : AlgIterator {
271
272 // Maximum insertion count and measure difference calculation depend on insertion style.
273 let (max_insertions, warn_insertions) = match (state.iteration(), config.bootstrap_insertions) {
274 (i, Some((l, k))) if i <= l => (k, false),
275 _ => (config.max_insertions, !state.is_quiet()),
276 };
277
278 let ω0 = match ν_delta {
279 None => op𝒟.apply(μ_base),
280 Some(ν) => op𝒟.apply(μ_base + ν),
281 };
282
283 // Add points to support until within error tolerance or maximum insertion count reached.
284 let mut count = 0;
285 let (within_tolerances, d) = 'insertion: loop {
286 if μ.len() > 0 {
287 // Form finite-dimensional subproblem. The subproblem references to the original μ^k
288 // from the beginning of the iteration are all contained in the immutable c and g.
289 // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional
290 // problems have not yet been updated to sign change.
291 let à = op𝒟.findim_matrix(μ.iter_locations());
292 let g̃ = DVector::from_iterator(μ.len(),
293 μ.iter_locations()
294 .map(|ζ| ω0.apply(ζ) - τv.apply(ζ))
295 .map(F::to_nalgebra_mixed));
296 let mut x = μ.masses_dvector();
297
298 // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃.
299 // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩
300 // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ
301 // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2
302 // = n |𝒟| |x|_2, where n is the number of points. Therefore
303 let Ã_normest = op𝒟norm * F::cast_from(μ.len());
304
305 // Solve finite-dimensional subproblem.
306 stats.inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config);
307
308 // Update masses of μ based on solution of finite-dimensional subproblem.
309 μ.set_masses_dvector(&x);
310 }
311
312 // Form d = τv + 𝒟μ - ω0 = τv + 𝒟(μ - μ^k) for checking the proximate optimality
313 // conditions in the predual space, and finding new points for insertion, if necessary.
314 let mut d = τv + match ν_delta {
315 None => op𝒟.preapply(μ.sub_matching(μ_base)),
316 Some(ν) => op𝒟.preapply(μ.sub_matching(μ_base) - ν)
317 };
318
319 // If no merging heuristic is used, let's be more conservative about spike insertion,
320 // and skip it after first round. If merging is done, being more greedy about spike
321 // insertion also seems to improve performance.
322 let skip_by_rough_check = if let SpikeMergingMethod::None = config.merging {
323 false
324 } else {
325 count > 0
326 };
327
328 // Find a spike to insert, if needed
329 let (ξ, _v_ξ, in_bounds) = match reg.find_tolerance_violation(
330 &mut d, τ, ε, skip_by_rough_check, config
331 ) {
332 None => break 'insertion (true, d),
333 Some(res) => res,
334 };
335
336 // Break if maximum insertion count reached
337 if count >= max_insertions {
338 break 'insertion (in_bounds, d)
339 }
340
341 // No point in optimising the weight here; the finite-dimensional algorithm is fast.
342 *μ += DeltaMeasure { x : ξ, α : 0.0 };
343 count += 1;
344 stats.inserted += 1;
345 };
346
347 if !within_tolerances && warn_insertions {
348 // Complain (but continue) if we failed to get within tolerances
349 // by inserting more points.
350 let err = format!("Maximum insertions reached without achieving \
351 subproblem solution tolerance");
352 println!("{}", err.red());
353 }
354
355 (d, within_tolerances)
356 } 138 }
357 139
358 pub(crate) fn prune_with_stats<F : Float, const N : usize>( 140 pub(crate) fn prune_with_stats<F : Float, const N : usize>(
359 μ : &mut RNDM<F, N>, 141 μ : &mut RNDM<F, N>,
360 ) -> usize { 142 ) -> usize {
407 /// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions. 189 /// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions.
408 /// 190 ///
409 /// Returns the final iterate. 191 /// Returns the final iterate.
410 #[replace_float_literals(F::cast_from(literal))] 192 #[replace_float_literals(F::cast_from(literal))]
411 pub fn pointsource_fb_reg< 193 pub fn pointsource_fb_reg<
412 'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Reg, const N : usize 194 F, I, A, Reg, P, const N : usize
413 >( 195 >(
414 opA : &'a A, 196 opA : &A,
415 b : &A::Observable, 197 b : &A::Observable,
416 reg : Reg, 198 reg : Reg,
417 op𝒟 : &'a 𝒟, 199 prox_penalty : &P,
418 fbconfig : &FBConfig<F>, 200 fbconfig : &FBConfig<F>,
419 iterator : I, 201 iterator : I,
420 mut plotter : SeqPlotter<F, N>, 202 mut plotter : SeqPlotter<F, N>,
421 ) -> RNDM<F, N> 203 ) -> RNDM<F, N>
422 where F : Float + ToNalgebraRealField, 204 where
423 I : AlgIteratorFactory<IterInfo<F, N>>, 205 F : Float + ToNalgebraRealField,
424 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>, 206 I : AlgIteratorFactory<IterInfo<F, N>>,
425 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 207 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
426 A : ForwardModel<RNDM<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>> 208 A : ForwardModel<RNDM<F, N>, F>
427 + AdjointProductBoundedBy<RNDM<F, N>, 𝒟, FloatType=F>, 209 + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType=F>,
428 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, 210 A::PreadjointCodomain : RealMapping<F, N>,
429 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone, 211 PlotLookup : Plotting<N>,
430 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>, 212 RNDM<F, N> : SpikeMerging<F>,
431 𝒟::Codomain : RealMapping<F, N>, 213 Reg : RegTerm<F, N>,
432 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 214 P : ProxPenalty<F, A::PreadjointCodomain, Reg, N>,
433 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 215 {
434 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
435 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
436 PlotLookup : Plotting<N>,
437 RNDM<F, N> : SpikeMerging<F>,
438 Reg : RegTerm<F, N> {
439 216
440 // Set up parameters 217 // Set up parameters
441 let config = &fbconfig.generic; 218 let config = &fbconfig.generic;
442 let op𝒟norm = op𝒟.opnorm_bound(Radon, Linfinity); 219 let τ = fbconfig.τ0/opA.adjoint_product_bound(prox_penalty).unwrap();
443 let τ = fbconfig.τ0/opA.adjoint_product_bound(&op𝒟).unwrap();
444 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled 220 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
445 // by τ compared to the conditional gradient approach. 221 // by τ compared to the conditional gradient approach.
446 let tolerance = config.tolerance * τ * reg.tolerance_scaling(); 222 let tolerance = config.tolerance * τ * reg.tolerance_scaling();
447 let mut ε = tolerance.initial(); 223 let mut ε = tolerance.initial();
448 224
463 let mut stats = IterInfo::new(); 239 let mut stats = IterInfo::new();
464 240
465 // Run the algorithm 241 // Run the algorithm
466 for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, stats.clone())) { 242 for state in iterator.iter_init(|| full_stats(&residual, &μ, ε, stats.clone())) {
467 // Calculate smooth part of surrogate model. 243 // Calculate smooth part of surrogate model.
468 let τv = opA.preadjoint().apply(residual * τ); 244 let mut τv = opA.preadjoint().apply(residual * τ);
469 245
470 // Save current base point 246 // Save current base point
471 let μ_base = μ.clone(); 247 let μ_base = μ.clone();
472 248
473 // Insert and reweigh 249 // Insert and reweigh
474 let (d, _within_tolerances) = insert_and_reweigh( 250 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
475 &mut μ, &τv, &μ_base, None, 251 &mut μ, &mut τv, &μ_base, None,
476 op𝒟, op𝒟norm,
477 τ, ε, 252 τ, ε,
478 config, &reg, &state, &mut stats 253 config, &reg, &state, &mut stats
479 ); 254 );
480 255
481 // Prune and possibly merge spikes 256 // Prune and possibly merge spikes
482 if config.merge_now(&state) { 257 if config.merge_now(&state) {
483 stats.merged += μ.merge_spikes(config.merging, |μ_candidate| { 258 stats.merged += prox_penalty.merge_spikes(&mut μ, &mut τv, &μ_base, τ, ε, config, &reg);
484 let mut d = &τv + op𝒟.preapply(μ_candidate.sub_matching(&μ_base));
485 reg.verify_merge_candidate(&mut d, μ_candidate, τ, ε, &config)
486 });
487 } 259 }
260
488 stats.pruned += prune_with_stats(&mut μ); 261 stats.pruned += prune_with_stats(&mut μ);
489 262
490 // Update residual 263 // Update residual
491 residual = calculate_residual(&μ, opA, b); 264 residual = calculate_residual(&μ, opA, b);
492 265
493 let iter = state.iteration(); 266 let iter = state.iteration();
494 stats.this_iters += 1; 267 stats.this_iters += 1;
495 268
496 // Give statistics if needed 269 // Give statistics if needed
497 state.if_verbose(|| { 270 state.if_verbose(|| {
498 plotter.plot_spikes(iter, Some(&d), Some(&τv), &μ); 271 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ);
499 full_stats(&residual, &μ, ε, std::mem::replace(&mut stats, IterInfo::new())) 272 full_stats(&residual, &μ, ε, std::mem::replace(&mut stats, IterInfo::new()))
500 }); 273 });
501 274
502 // Update main tolerance for next iteration 275 // Update main tolerance for next iteration
503 ε = tolerance.update(ε, iter); 276 ε = tolerance.update(ε, iter);
524 /// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions. 297 /// of [`std::sync::Arc`] to only update relevant parts of the bisection tree when adding functions.
525 /// 298 ///
526 /// Returns the final iterate. 299 /// Returns the final iterate.
527 #[replace_float_literals(F::cast_from(literal))] 300 #[replace_float_literals(F::cast_from(literal))]
528 pub fn pointsource_fista_reg< 301 pub fn pointsource_fista_reg<
529 'a, F, I, A, GA, 𝒟, BTA, G𝒟, S, K, Reg, const N : usize 302 F, I, A, Reg, P, const N : usize
530 >( 303 >(
531 opA : &'a A, 304 opA : &A,
532 b : &A::Observable, 305 b : &A::Observable,
533 reg : Reg, 306 reg : Reg,
534 op𝒟 : &'a 𝒟, 307 prox_penalty : &P,
535 fbconfig : &FBConfig<F>, 308 fbconfig : &FBConfig<F>,
536 iterator : I, 309 iterator : I,
537 mut plotter : SeqPlotter<F, N>, 310 mut plotter : SeqPlotter<F, N>,
538 ) -> RNDM<F, N> 311 ) -> RNDM<F, N>
539 where F : Float + ToNalgebraRealField, 312 where
540 I : AlgIteratorFactory<IterInfo<F, N>>, 313 F : Float + ToNalgebraRealField,
541 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>, 314 I : AlgIteratorFactory<IterInfo<F, N>>,
542 GA : SupportGenerator<F, N, SupportType = S, Id = usize> + Clone, 315 for<'b> &'b A::Observable : std::ops::Neg<Output=A::Observable>,
543 A : ForwardModel<RNDM<F, N>, F, PreadjointCodomain = BTFN<F, GA, BTA, N>> 316 A : ForwardModel<RNDM<F, N>, F>
544 + AdjointProductBoundedBy<RNDM<F, N>, 𝒟, FloatType=F>, 317 + AdjointProductBoundedBy<RNDM<F, N>, P, FloatType=F>,
545 BTA : BTSearch<F, N, Data=usize, Agg=Bounds<F>>, 318 A::PreadjointCodomain : RealMapping<F, N>,
546 G𝒟 : SupportGenerator<F, N, SupportType = K, Id = usize> + Clone, 319 PlotLookup : Plotting<N>,
547 𝒟 : DiscreteMeasureOp<Loc<F, N>, F, PreCodomain = PreBTFN<F, G𝒟, N>>, 320 RNDM<F, N> : SpikeMerging<F>,
548 𝒟::Codomain : RealMapping<F, N>, 321 Reg : RegTerm<F, N>,
549 S: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 322 P : ProxPenalty<F, A::PreadjointCodomain, Reg, N>,
550 K: RealMapping<F, N> + LocalAnalysis<F, Bounds<F>, N>, 323 {
551 BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
552 Cube<F, N>: P2Minimise<Loc<F, N>, F>,
553 PlotLookup : Plotting<N>,
554 RNDM<F, N> : SpikeMerging<F>,
555 Reg : RegTerm<F, N> {
556 324
557 // Set up parameters 325 // Set up parameters
558 let config = &fbconfig.generic; 326 let config = &fbconfig.generic;
559 let op𝒟norm = op𝒟.opnorm_bound(Radon, Linfinity); 327 let τ = fbconfig.τ0/opA.adjoint_product_bound(prox_penalty).unwrap();
560 let τ = fbconfig.τ0/opA.adjoint_product_bound(&op𝒟).unwrap();
561 let mut λ = 1.0; 328 let mut λ = 1.0;
562 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled 329 // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled
563 // by τ compared to the conditional gradient approach. 330 // by τ compared to the conditional gradient approach.
564 let tolerance = config.tolerance * τ * reg.tolerance_scaling(); 331 let tolerance = config.tolerance * τ * reg.tolerance_scaling();
565 let mut ε = tolerance.initial(); 332 let mut ε = tolerance.initial();
581 let mut stats = IterInfo::new(); 348 let mut stats = IterInfo::new();
582 349
583 // Run the algorithm 350 // Run the algorithm
584 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { 351 for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) {
585 // Calculate smooth part of surrogate model. 352 // Calculate smooth part of surrogate model.
586 let τv = opA.preadjoint().apply(residual * τ); 353 let mut τv = opA.preadjoint().apply(residual * τ);
587 354
588 // Save current base point 355 // Save current base point
589 let μ_base = μ.clone(); 356 let μ_base = μ.clone();
590 357
591 // Insert new spikes and reweigh 358 // Insert new spikes and reweigh
592 let (d, _within_tolerances) = insert_and_reweigh( 359 let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh(
593 &mut μ, &τv, &μ_base, None, 360 &mut μ, &mut τv, &μ_base, None,
594 op𝒟, op𝒟norm,
595 τ, ε, 361 τ, ε,
596 config, &reg, &state, &mut stats 362 config, &reg, &state, &mut stats
597 ); 363 );
598 364
599 // (Do not) merge spikes. 365 // (Do not) merge spikes.
630 let iter = state.iteration(); 396 let iter = state.iteration();
631 stats.this_iters += 1; 397 stats.this_iters += 1;
632 398
633 // Give statistics if needed 399 // Give statistics if needed
634 state.if_verbose(|| { 400 state.if_verbose(|| {
635 plotter.plot_spikes(iter, Some(&d), Some(&τv), &μ_prev); 401 plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ_prev);
636 full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new())) 402 full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new()))
637 }); 403 });
638 404
639 // Update main tolerance for next iteration 405 // Update main tolerance for next iteration
640 ε = tolerance.update(ε, iter); 406 ε = tolerance.update(ε, iter);

Mercurial Repository: pointsource_algs

mercurial