Thu, 08 Dec 2022 14:10:07 +0200
Save more CSV files when iteration-wise plotting is enabled.
This helps to generate TikZ illustrations for presentations.
0 | 1 | /*! |
2 | This module provides [`RunnableExperiment`] for running chosen algorithms on a chosen experiment. | |
3 | */ | |
4 | ||
5 | use numeric_literals::replace_float_literals; | |
6 | use colored::Colorize; | |
7 | use serde::{Serialize, Deserialize}; | |
8 | use serde_json; | |
9 | use nalgebra::base::DVector; | |
10 | use std::hash::Hash; | |
11 | use chrono::{DateTime, Utc}; | |
12 | use cpu_time::ProcessTime; | |
13 | use clap::ValueEnum; | |
14 | use std::collections::HashMap; | |
15 | use std::time::Instant; | |
16 | ||
17 | use rand::prelude::{ | |
18 | StdRng, | |
19 | SeedableRng | |
20 | }; | |
21 | use rand_distr::Distribution; | |
22 | ||
23 | use alg_tools::bisection_tree::*; | |
24 | use alg_tools::iterate::{ | |
25 | Timed, | |
26 | AlgIteratorOptions, | |
27 | Verbose, | |
28 | AlgIteratorFactory, | |
29 | }; | |
30 | use alg_tools::logger::Logger; | |
31 | use alg_tools::error::DynError; | |
32 | use alg_tools::tabledump::TableDump; | |
33 | use alg_tools::sets::Cube; | |
34 | use alg_tools::mapping::RealMapping; | |
35 | use alg_tools::nalgebra_support::ToNalgebraRealField; | |
36 | use alg_tools::euclidean::Euclidean; | |
37 | use alg_tools::norms::{Norm, L1}; | |
38 | use alg_tools::lingrid::lingrid; | |
39 | use alg_tools::sets::SetOrd; | |
40 | ||
41 | use crate::kernels::*; | |
42 | use crate::types::*; | |
43 | use crate::measures::*; | |
44 | use crate::measures::merging::SpikeMerging; | |
45 | use crate::forward_model::*; | |
46 | use crate::fb::{ | |
47 | FBConfig, | |
48 | pointsource_fb, | |
49 | FBMetaAlgorithm, FBGenericConfig, | |
50 | }; | |
51 | use crate::pdps::{ | |
52 | PDPSConfig, | |
53 | L2Squared, | |
54 | pointsource_pdps, | |
55 | }; | |
56 | use crate::frank_wolfe::{ | |
57 | FWConfig, | |
58 | FWVariant, | |
59 | pointsource_fw, | |
60 | prepare_optimise_weights, | |
61 | optimise_weights, | |
62 | }; | |
63 | use crate::subproblem::InnerSettings; | |
64 | use crate::seminorms::*; | |
65 | use crate::plot::*; | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
66 | use crate::{AlgorithmOverrides, CommandLineArgs}; |
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
67 | use crate::tolerance::Tolerance; |
0 | 68 | |
69 | /// Available algorithms and their configurations | |
70 | #[derive(Copy, Clone, Debug, Serialize, Deserialize)] | |
71 | pub enum AlgorithmConfig<F : Float> { | |
72 | FB(FBConfig<F>), | |
73 | FW(FWConfig<F>), | |
74 | PDPS(PDPSConfig<F>), | |
75 | } | |
76 | ||
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
77 | fn unpack_tolerance<F : Float>(v : &Vec<F>) -> Tolerance<F> { |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
78 | assert!(v.len() == 3); |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
79 | Tolerance::Power { initial : v[0], factor : v[1], exponent : v[2] } |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
80 | } |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
81 | |
0 | 82 | impl<F : ClapFloat> AlgorithmConfig<F> { |
83 | /// Override supported parameters based on the command line. | |
84 | pub fn cli_override(self, cli : &AlgorithmOverrides<F>) -> Self { | |
85 | let override_fb_generic = |g : FBGenericConfig<F>| { | |
86 | FBGenericConfig { | |
87 | bootstrap_insertions : cli.bootstrap_insertions | |
88 | .as_ref() | |
89 | .map_or(g.bootstrap_insertions, | |
90 | |n| Some((n[0], n[1]))), | |
91 | merge_every : cli.merge_every.unwrap_or(g.merge_every), | |
92 | merging : cli.merging.clone().unwrap_or(g.merging), | |
93 | final_merging : cli.final_merging.clone().unwrap_or(g.final_merging), | |
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
94 | tolerance: cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(g.tolerance), |
0 | 95 | .. g |
96 | } | |
97 | }; | |
98 | ||
99 | use AlgorithmConfig::*; | |
100 | match self { | |
101 | FB(fb) => FB(FBConfig { | |
102 | τ0 : cli.tau0.unwrap_or(fb.τ0), | |
103 | insertion : override_fb_generic(fb.insertion), | |
104 | .. fb | |
105 | }), | |
106 | PDPS(pdps) => PDPS(PDPSConfig { | |
107 | τ0 : cli.tau0.unwrap_or(pdps.τ0), | |
108 | σ0 : cli.sigma0.unwrap_or(pdps.σ0), | |
109 | acceleration : cli.acceleration.unwrap_or(pdps.acceleration), | |
110 | insertion : override_fb_generic(pdps.insertion), | |
111 | .. pdps | |
112 | }), | |
113 | FW(fw) => FW(FWConfig { | |
114 | merging : cli.merging.clone().unwrap_or(fw.merging), | |
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
115 | tolerance : cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(fw.tolerance), |
0 | 116 | .. fw |
117 | }) | |
118 | } | |
119 | } | |
120 | } | |
121 | ||
122 | /// Helper struct for tagging and [`AlgorithmConfig`] or [`Experiment`] with a name. | |
123 | #[derive(Clone, Debug, Serialize, Deserialize)] | |
124 | pub struct Named<Data> { | |
125 | pub name : String, | |
126 | #[serde(flatten)] | |
127 | pub data : Data, | |
128 | } | |
129 | ||
130 | /// Shorthand algorithm configurations, to be used with the command line parser | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
131 | #[derive(ValueEnum, Debug, Copy, Clone, Eq, PartialEq, Hash, Serialize, Deserialize)] |
0 | 132 | pub enum DefaultAlgorithm { |
133 | /// The μFB forward-backward method | |
134 | #[clap(name = "fb")] | |
135 | FB, | |
136 | /// The μFISTA inertial forward-backward method | |
137 | #[clap(name = "fista")] | |
138 | FISTA, | |
139 | /// The “fully corrective” conditional gradient method | |
140 | #[clap(name = "fw")] | |
141 | FW, | |
142 | /// The “relaxed conditional gradient method | |
143 | #[clap(name = "fwrelax")] | |
144 | FWRelax, | |
145 | /// The μPDPS primal-dual proximal splitting method | |
146 | #[clap(name = "pdps")] | |
147 | PDPS, | |
148 | } | |
149 | ||
150 | impl DefaultAlgorithm { | |
151 | /// Returns the algorithm configuration corresponding to the algorithm shorthand | |
152 | pub fn default_config<F : Float>(&self) -> AlgorithmConfig<F> { | |
153 | use DefaultAlgorithm::*; | |
154 | match *self { | |
155 | FB => AlgorithmConfig::FB(Default::default()), | |
156 | FISTA => AlgorithmConfig::FB(FBConfig{ | |
157 | meta : FBMetaAlgorithm::InertiaFISTA, | |
158 | .. Default::default() | |
159 | }), | |
160 | FW => AlgorithmConfig::FW(Default::default()), | |
161 | FWRelax => AlgorithmConfig::FW(FWConfig{ | |
162 | variant : FWVariant::Relaxed, | |
163 | .. Default::default() | |
164 | }), | |
165 | PDPS => AlgorithmConfig::PDPS(Default::default()), | |
166 | } | |
167 | } | |
168 | ||
169 | /// Returns the [`Named`] algorithm corresponding to the algorithm shorthand | |
170 | pub fn get_named<F : Float>(&self) -> Named<AlgorithmConfig<F>> { | |
171 | self.to_named(self.default_config()) | |
172 | } | |
173 | ||
174 | pub fn to_named<F : Float>(self, alg : AlgorithmConfig<F>) -> Named<AlgorithmConfig<F>> { | |
175 | let name = self.to_possible_value().unwrap().get_name().to_string(); | |
176 | Named{ name , data : alg } | |
177 | } | |
178 | } | |
179 | ||
180 | ||
181 | // // Floats cannot be hashed directly, so just hash the debug formatting | |
182 | // // for use as file identifier. | |
183 | // impl<F : Float> Hash for AlgorithmConfig<F> { | |
184 | // fn hash<H: Hasher>(&self, state: &mut H) { | |
185 | // format!("{:?}", self).hash(state); | |
186 | // } | |
187 | // } | |
188 | ||
189 | /// Plotting level configuration | |
190 | #[derive(Copy, Clone, Eq, PartialEq, Ord, PartialOrd, Serialize, ValueEnum, Debug)] | |
191 | pub enum PlotLevel { | |
192 | /// Plot nothing | |
193 | #[clap(name = "none")] | |
194 | None, | |
195 | /// Plot problem data | |
196 | #[clap(name = "data")] | |
197 | Data, | |
198 | /// Plot iterationwise state | |
199 | #[clap(name = "iter")] | |
200 | Iter, | |
201 | } | |
202 | ||
203 | type DefaultBT<F, const N : usize> = BT< | |
204 | DynamicDepth, | |
205 | F, | |
206 | usize, | |
207 | Bounds<F>, | |
208 | N | |
209 | >; | |
210 | type DefaultSeminormOp<F, K, const N : usize> = ConvolutionOp<F, K, DefaultBT<F, N>, N>; | |
211 | type DefaultSG<F, Sensor, Spread, const N : usize> = SensorGrid::< | |
212 | F, | |
213 | Sensor, | |
214 | Spread, | |
215 | DefaultBT<F, N>, | |
216 | N | |
217 | >; | |
218 | ||
219 | /// This is a dirty workaround to rust-csv not supporting struct flattening etc. | |
220 | #[derive(Serialize)] | |
221 | struct CSVLog<F> { | |
222 | iter : usize, | |
223 | cpu_time : f64, | |
224 | value : F, | |
225 | post_value : F, | |
226 | n_spikes : usize, | |
227 | inner_iters : usize, | |
228 | merged : usize, | |
229 | pruned : usize, | |
230 | this_iters : usize, | |
231 | } | |
232 | ||
233 | /// Collected experiment statistics | |
234 | #[derive(Clone, Debug, Serialize)] | |
235 | struct ExperimentStats<F : Float> { | |
236 | /// Signal-to-noise ratio in decibels | |
237 | ssnr : F, | |
238 | /// Proportion of noise in the signal as a number in $[0, 1]$. | |
239 | noise_ratio : F, | |
240 | /// When the experiment was run (UTC) | |
241 | when : DateTime<Utc>, | |
242 | } | |
243 | ||
244 | #[replace_float_literals(F::cast_from(literal))] | |
245 | impl<F : Float> ExperimentStats<F> { | |
246 | /// Calculate [`ExperimentStats`] based on a noisy `signal` and the separated `noise` signal. | |
247 | fn new<E : Euclidean<F>>(signal : &E, noise : &E) -> Self { | |
248 | let s = signal.norm2_squared(); | |
249 | let n = noise.norm2_squared(); | |
250 | let noise_ratio = (n / s).sqrt(); | |
251 | let ssnr = 10.0 * (s / n).log10(); | |
252 | ExperimentStats { | |
253 | ssnr, | |
254 | noise_ratio, | |
255 | when : Utc::now(), | |
256 | } | |
257 | } | |
258 | } | |
259 | /// Collected algorithm statistics | |
260 | #[derive(Clone, Debug, Serialize)] | |
261 | struct AlgorithmStats<F : Float> { | |
262 | /// Overall CPU time spent | |
263 | cpu_time : F, | |
264 | /// Real time spent | |
265 | elapsed : F | |
266 | } | |
267 | ||
268 | ||
269 | /// A wrapper for [`serde_json::to_writer_pretty`] that takes a filename as input | |
270 | /// and outputs a [`DynError`]. | |
271 | fn write_json<T : Serialize>(filename : String, data : &T) -> DynError { | |
272 | serde_json::to_writer_pretty(std::fs::File::create(filename)?, data)?; | |
273 | Ok(()) | |
274 | } | |
275 | ||
276 | ||
277 | /// Struct for experiment configurations | |
278 | #[derive(Debug, Clone, Serialize)] | |
279 | pub struct Experiment<F, NoiseDistr, S, K, P, const N : usize> | |
280 | where F : Float, | |
281 | [usize; N] : Serialize, | |
282 | NoiseDistr : Distribution<F>, | |
283 | S : Sensor<F, N>, | |
284 | P : Spread<F, N>, | |
285 | K : SimpleConvolutionKernel<F, N>, | |
286 | { | |
287 | /// Domain $Ω$. | |
288 | pub domain : Cube<F, N>, | |
289 | /// Number of sensors along each dimension | |
290 | pub sensor_count : [usize; N], | |
291 | /// Noise distribution | |
292 | pub noise_distr : NoiseDistr, | |
293 | /// Seed for random noise generation (for repeatable experiments) | |
294 | pub noise_seed : u64, | |
295 | /// Sensor $θ$; $θ * ψ$ forms the forward operator $𝒜$. | |
296 | pub sensor : S, | |
297 | /// Spread $ψ$; $θ * ψ$ forms the forward operator $𝒜$. | |
298 | pub spread : P, | |
299 | /// Kernel $ρ$ of $𝒟$. | |
300 | pub kernel : K, | |
301 | /// True point sources | |
302 | pub μ_hat : DiscreteMeasure<Loc<F, N>, F>, | |
303 | /// Regularisation parameter | |
304 | pub α : F, | |
305 | /// For plotting : how wide should the kernels be plotted | |
306 | pub kernel_plot_width : F, | |
307 | /// Data term | |
308 | pub dataterm : DataTerm, | |
309 | /// A map of default configurations for algorithms | |
310 | #[serde(skip)] | |
311 | pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>, | |
312 | } | |
313 | ||
314 | /// Trait for runnable experiments | |
315 | pub trait RunnableExperiment<F : ClapFloat> { | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
316 | /// Run all algorithms provided, or default algorithms if none provided, on the experiment. |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
317 | fn runall(&self, cli : &CommandLineArgs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
318 | algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError; |
0 | 319 | |
320 | /// Return algorithm default config | |
321 | fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>) | |
322 | -> Named<AlgorithmConfig<F>>; | |
323 | } | |
324 | ||
325 | impl<F, NoiseDistr, S, K, P, const N : usize> RunnableExperiment<F> for | |
326 | Named<Experiment<F, NoiseDistr, S, K, P, N>> | |
327 | where F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>, | |
328 | [usize; N] : Serialize, | |
329 | S : Sensor<F, N> + Copy + Serialize, | |
330 | P : Spread<F, N> + Copy + Serialize, | |
331 | Convolution<S, P>: Spread<F, N> + Bounded<F> + LocalAnalysis<F, Bounds<F>, N> + Copy, | |
332 | AutoConvolution<P> : BoundedBy<F, K>, | |
333 | K : SimpleConvolutionKernel<F, N> + LocalAnalysis<F, Bounds<F>, N> + Copy + Serialize, | |
334 | Cube<F, N>: P2Minimise<Loc<F, N>, F> + SetOrd, | |
335 | PlotLookup : Plotting<N>, | |
336 | DefaultBT<F, N> : SensorGridBT<F, S, P, N, Depth=DynamicDepth> + BTSearch<F, N>, | |
337 | BTNodeLookup: BTNode<F, usize, Bounds<F>, N>, | |
338 | DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>, | |
339 | NoiseDistr : Distribution<F> + Serialize { | |
340 | ||
341 | fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>) | |
342 | -> Named<AlgorithmConfig<F>> { | |
343 | alg.to_named( | |
344 | self.data | |
345 | .algorithm_defaults | |
346 | .get(&alg) | |
347 | .map_or_else(|| alg.default_config(), | |
348 | |config| config.clone()) | |
349 | .cli_override(cli) | |
350 | ) | |
351 | } | |
352 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
353 | fn runall(&self, cli : &CommandLineArgs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
354 | algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError { |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
355 | // Get experiment configuration |
0 | 356 | let &Named { |
357 | name : ref experiment_name, | |
358 | data : Experiment { | |
359 | domain, sensor_count, ref noise_distr, sensor, spread, kernel, | |
360 | ref μ_hat, α, kernel_plot_width, dataterm, noise_seed, | |
361 | .. | |
362 | } | |
363 | } = self; | |
364 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
365 | // Set up output directory |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
366 | let prefix = format!("{}/{}/", cli.outdir, self.name); |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
367 | |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
368 | // Set up algorithms |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
369 | let iterator_options = AlgIteratorOptions{ |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
370 | max_iter : cli.max_iter, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
371 | verbose_iter : cli.verbose_iter |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
372 | .map_or(Verbose::Logarithmic(10), |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
373 | |n| Verbose::Every(n)), |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
374 | quiet : cli.quiet, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
375 | }; |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
376 | let algorithms = match (algs, self.data.dataterm) { |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
377 | (Some(algs), _) => algs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
378 | (None, DataTerm::L2Squared) => vec![DefaultAlgorithm::FB.get_named()], |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
379 | (None, DataTerm::L1) => vec![DefaultAlgorithm::PDPS.get_named()], |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
380 | }; |
0 | 381 | |
382 | // Set up operators | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
383 | let depth = DynamicDepth(8); |
0 | 384 | let opA = DefaultSG::new(domain, sensor_count, sensor, spread, depth); |
385 | let op𝒟 = DefaultSeminormOp::new(depth, domain, kernel); | |
386 | ||
387 | // Set up random number generator. | |
388 | let mut rng = StdRng::seed_from_u64(noise_seed); | |
389 | ||
390 | // Generate the data and calculate SSNR statistic | |
391 | let b_hat = opA.apply(μ_hat); | |
392 | let noise = DVector::from_distribution(b_hat.len(), &noise_distr, &mut rng); | |
393 | let b = &b_hat + &noise; | |
394 | // Need to wrap calc_ssnr into a function to hide ultra-lame nalgebra::RealField | |
395 | // overloading log10 and conflicting with standard NumTraits one. | |
396 | let stats = ExperimentStats::new(&b, &noise); | |
397 | ||
398 | // Save experiment configuration and statistics | |
399 | let mkname_e = |t| format!("{prefix}{t}.json", prefix = prefix, t = t); | |
400 | std::fs::create_dir_all(&prefix)?; | |
401 | write_json(mkname_e("experiment"), self)?; | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
402 | write_json(mkname_e("config"), cli)?; |
0 | 403 | write_json(mkname_e("stats"), &stats)?; |
404 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
405 | plotall(cli, &prefix, &domain, &sensor, &kernel, &spread, |
0 | 406 | &μ_hat, &op𝒟, &opA, &b_hat, &b, kernel_plot_width)?; |
407 | ||
408 | // Run the algorithm(s) | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
409 | for named @ Named { name : alg_name, data : alg } in algorithms.iter() { |
0 | 410 | let this_prefix = format!("{}{}/", prefix, alg_name); |
411 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
412 | let running = || if !cli.quiet { |
0 | 413 | println!("{}\n{}\n{}", |
414 | format!("Running {} on experiment {}…", alg_name, experiment_name).cyan(), | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
415 | format!("{:?}", iterator_options).bright_black(), |
0 | 416 | format!("{:?}", alg).bright_black()); |
417 | }; | |
418 | ||
419 | // Create Logger and IteratorFactory | |
420 | let mut logger = Logger::new(); | |
421 | let findim_data = prepare_optimise_weights(&opA); | |
422 | let inner_config : InnerSettings<F> = Default::default(); | |
423 | let inner_it = inner_config.iterator_options; | |
424 | let logmap = |iter, Timed { cpu_time, data }| { | |
425 | let IterInfo { | |
426 | value, | |
427 | n_spikes, | |
428 | inner_iters, | |
429 | merged, | |
430 | pruned, | |
431 | postprocessing, | |
432 | this_iters, | |
433 | .. | |
434 | } = data; | |
435 | let post_value = match postprocessing { | |
436 | None => value, | |
437 | Some(mut μ) => { | |
438 | match dataterm { | |
439 | DataTerm::L2Squared => { | |
440 | optimise_weights( | |
441 | &mut μ, &opA, &b, α, &findim_data, &inner_config, | |
442 | inner_it | |
443 | ); | |
444 | dataterm.value_at_residual(opA.apply(&μ) - &b) + α * μ.norm(Radon) | |
445 | }, | |
446 | _ => value, | |
447 | } | |
448 | } | |
449 | }; | |
450 | CSVLog { | |
451 | iter, | |
452 | value, | |
453 | post_value, | |
454 | n_spikes, | |
455 | cpu_time : cpu_time.as_secs_f64(), | |
456 | inner_iters, | |
457 | merged, | |
458 | pruned, | |
459 | this_iters | |
460 | } | |
461 | }; | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
462 | let iterator = iterator_options.instantiate() |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
463 | .timed() |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
464 | .mapped(logmap) |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
465 | .into_log(&mut logger); |
0 | 466 | let plotgrid = lingrid(&domain, &[if N==1 { 1000 } else { 100 }; N]); |
467 | ||
468 | // Create plotter and directory if needed. | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
469 | let plot_count = if cli.plot >= PlotLevel::Iter { 2000 } else { 0 }; |
0 | 470 | let plotter = SeqPlotter::new(this_prefix, plot_count, plotgrid); |
471 | ||
472 | // Run the algorithm | |
473 | let start = Instant::now(); | |
474 | let start_cpu = ProcessTime::now(); | |
475 | let μ : DiscreteMeasure<Loc<F, N>, F> = match (alg, dataterm) { | |
476 | (AlgorithmConfig::FB(ref algconfig), DataTerm::L2Squared) => { | |
477 | running(); | |
478 | pointsource_fb(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter) | |
479 | }, | |
480 | (AlgorithmConfig::FW(ref algconfig), DataTerm::L2Squared) => { | |
481 | running(); | |
482 | pointsource_fw(&opA, &b, α, &algconfig, iterator, plotter) | |
483 | }, | |
484 | (AlgorithmConfig::PDPS(ref algconfig), DataTerm::L2Squared) => { | |
485 | running(); | |
486 | pointsource_pdps(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter, L2Squared) | |
487 | }, | |
488 | (AlgorithmConfig::PDPS(ref algconfig), DataTerm::L1) => { | |
489 | running(); | |
490 | pointsource_pdps(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter, L1) | |
491 | }, | |
492 | _ => { | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
493 | let msg = format!("Algorithm “{alg_name}” not implemented for \ |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
494 | dataterm {dataterm:?}. Skipping.").red(); |
0 | 495 | eprintln!("{}", msg); |
496 | continue | |
497 | } | |
498 | }; | |
499 | let elapsed = start.elapsed().as_secs_f64(); | |
500 | let cpu_time = start_cpu.elapsed().as_secs_f64(); | |
501 | ||
502 | println!("{}", format!("Elapsed {elapsed}s (CPU time {cpu_time}s)… ").yellow()); | |
503 | ||
504 | // Save results | |
505 | println!("{}", "Saving results…".green()); | |
506 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
507 | let mkname = |t| format!("{prefix}{alg_name}_{t}"); |
0 | 508 | |
509 | write_json(mkname("config.json"), &named)?; | |
510 | write_json(mkname("stats.json"), &AlgorithmStats { cpu_time, elapsed })?; | |
511 | μ.write_csv(mkname("reco.txt"))?; | |
512 | logger.write_csv(mkname("log.txt"))?; | |
513 | } | |
514 | ||
515 | Ok(()) | |
516 | } | |
517 | } | |
518 | ||
519 | /// Plot experiment setup | |
520 | #[replace_float_literals(F::cast_from(literal))] | |
521 | fn plotall<F, Sensor, Kernel, Spread, 𝒟, A, const N : usize>( | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
522 | cli : &CommandLineArgs, |
0 | 523 | prefix : &String, |
524 | domain : &Cube<F, N>, | |
525 | sensor : &Sensor, | |
526 | kernel : &Kernel, | |
527 | spread : &Spread, | |
528 | μ_hat : &DiscreteMeasure<Loc<F, N>, F>, | |
529 | op𝒟 : &𝒟, | |
530 | opA : &A, | |
531 | b_hat : &A::Observable, | |
532 | b : &A::Observable, | |
533 | kernel_plot_width : F, | |
534 | ) -> DynError | |
535 | where F : Float + ToNalgebraRealField, | |
536 | Sensor : RealMapping<F, N> + Support<F, N> + Clone, | |
537 | Spread : RealMapping<F, N> + Support<F, N> + Clone, | |
538 | Kernel : RealMapping<F, N> + Support<F, N>, | |
539 | Convolution<Sensor, Spread> : RealMapping<F, N> + Support<F, N>, | |
540 | 𝒟 : DiscreteMeasureOp<Loc<F, N>, F>, | |
541 | 𝒟::Codomain : RealMapping<F, N>, | |
542 | A : ForwardModel<Loc<F, N>, F>, | |
543 | A::PreadjointCodomain : RealMapping<F, N> + Bounded<F>, | |
544 | PlotLookup : Plotting<N>, | |
545 | Cube<F, N> : SetOrd { | |
546 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
547 | if cli.plot < PlotLevel::Data { |
0 | 548 | return Ok(()) |
549 | } | |
550 | ||
551 | let base = Convolution(sensor.clone(), spread.clone()); | |
552 | ||
553 | let resolution = if N==1 { 100 } else { 40 }; | |
554 | let pfx = |n| format!("{}{}", prefix, n); | |
555 | let plotgrid = lingrid(&[[-kernel_plot_width, kernel_plot_width]; N].into(), &[resolution; N]); | |
556 | ||
557 | PlotLookup::plot_into_file(sensor, plotgrid, pfx("sensor"), "sensor".to_string()); | |
558 | PlotLookup::plot_into_file(kernel, plotgrid, pfx("kernel"), "kernel".to_string()); | |
559 | PlotLookup::plot_into_file(spread, plotgrid, pfx("spread"), "spread".to_string()); | |
560 | PlotLookup::plot_into_file(&base, plotgrid, pfx("base_sensor"), "base_sensor".to_string()); | |
561 | ||
562 | let plotgrid2 = lingrid(&domain, &[resolution; N]); | |
563 | ||
564 | let ω_hat = op𝒟.apply(μ_hat); | |
565 | let noise = opA.preadjoint().apply(opA.apply(μ_hat) - b); | |
566 | PlotLookup::plot_into_file(&ω_hat, plotgrid2, pfx("omega_hat"), "ω̂".to_string()); | |
567 | PlotLookup::plot_into_file(&noise, plotgrid2, pfx("omega_noise"), | |
568 | "noise Aᵀ(Aμ̂ - b)".to_string()); | |
569 | ||
570 | let preadj_b = opA.preadjoint().apply(b); | |
571 | let preadj_b_hat = opA.preadjoint().apply(b_hat); | |
572 | //let bounds = preadj_b.bounds().common(&preadj_b_hat.bounds()); | |
573 | PlotLookup::plot_into_file_spikes( | |
574 | "Aᵀb".to_string(), &preadj_b, | |
575 | "Aᵀb̂".to_string(), Some(&preadj_b_hat), | |
576 | plotgrid2, None, &μ_hat, | |
577 | pfx("omega_b") | |
578 | ); | |
579 | ||
580 | // Save true solution and observables | |
581 | let pfx = |n| format!("{}{}", prefix, n); | |
582 | μ_hat.write_csv(pfx("orig.txt"))?; | |
583 | opA.write_observable(&b_hat, pfx("b_hat"))?; | |
584 | opA.write_observable(&b, pfx("b_noisy")) | |
585 | } | |
586 |