Sun, 11 Dec 2022 23:25:53 +0200
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
* Fixes the conditional gradient methods that were incorrectly solving
positivity constrained subproblems although the infinite-dimensional
versions do not include such constraints.
* Introduces the `ExperimentV2` struct that has `regularisation` in place
of `α`. The `Experiment` struct is now deprecated.
* The L^2-squared experiments were switch to be unconstrained, as the
Franke-Wolfe implementations do not support constraints. (This would
be easy to add for the “fully corrective” variant, but is not
immediate for the “relaxed” variant.)
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; | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
37 | use alg_tools::norms::L1; |
0 | 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, | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
48 | pointsource_fb_reg, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
49 | FBMetaAlgorithm, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
50 | FBGenericConfig, |
0 | 51 | }; |
52 | use crate::pdps::{ | |
53 | PDPSConfig, | |
54 | L2Squared, | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
55 | pointsource_pdps_reg, |
0 | 56 | }; |
57 | use crate::frank_wolfe::{ | |
58 | FWConfig, | |
59 | FWVariant, | |
60 | pointsource_fw, | |
61 | prepare_optimise_weights, | |
62 | optimise_weights, | |
63 | }; | |
64 | use crate::subproblem::InnerSettings; | |
65 | use crate::seminorms::*; | |
66 | use crate::plot::*; | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
67 | use crate::{AlgorithmOverrides, CommandLineArgs}; |
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
68 | use crate::tolerance::Tolerance; |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
69 | use crate::regularisation::{Regularisation, RadonRegTerm, NonnegRadonRegTerm}; |
0 | 70 | |
71 | /// Available algorithms and their configurations | |
72 | #[derive(Copy, Clone, Debug, Serialize, Deserialize)] | |
73 | pub enum AlgorithmConfig<F : Float> { | |
74 | FB(FBConfig<F>), | |
75 | FW(FWConfig<F>), | |
76 | PDPS(PDPSConfig<F>), | |
77 | } | |
78 | ||
20
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
79 | 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
|
80 | assert!(v.len() == 3); |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
81 | 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
|
82 | } |
90f77ad9a98d
Added command line option for (power) tolerance
Tuomo Valkonen <tuomov@iki.fi>
parents:
9
diff
changeset
|
83 | |
0 | 84 | impl<F : ClapFloat> AlgorithmConfig<F> { |
85 | /// Override supported parameters based on the command line. | |
86 | pub fn cli_override(self, cli : &AlgorithmOverrides<F>) -> Self { | |
87 | let override_fb_generic = |g : FBGenericConfig<F>| { | |
88 | FBGenericConfig { | |
89 | bootstrap_insertions : cli.bootstrap_insertions | |
90 | .as_ref() | |
91 | .map_or(g.bootstrap_insertions, | |
92 | |n| Some((n[0], n[1]))), | |
93 | merge_every : cli.merge_every.unwrap_or(g.merge_every), | |
94 | merging : cli.merging.clone().unwrap_or(g.merging), | |
95 | 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
|
96 | tolerance: cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(g.tolerance), |
0 | 97 | .. g |
98 | } | |
99 | }; | |
100 | ||
101 | use AlgorithmConfig::*; | |
102 | match self { | |
103 | FB(fb) => FB(FBConfig { | |
104 | τ0 : cli.tau0.unwrap_or(fb.τ0), | |
105 | insertion : override_fb_generic(fb.insertion), | |
106 | .. fb | |
107 | }), | |
108 | PDPS(pdps) => PDPS(PDPSConfig { | |
109 | τ0 : cli.tau0.unwrap_or(pdps.τ0), | |
110 | σ0 : cli.sigma0.unwrap_or(pdps.σ0), | |
111 | acceleration : cli.acceleration.unwrap_or(pdps.acceleration), | |
112 | insertion : override_fb_generic(pdps.insertion), | |
113 | .. pdps | |
114 | }), | |
115 | FW(fw) => FW(FWConfig { | |
116 | 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
|
117 | tolerance : cli.tolerance.as_ref().map(unpack_tolerance).unwrap_or(fw.tolerance), |
0 | 118 | .. fw |
119 | }) | |
120 | } | |
121 | } | |
122 | } | |
123 | ||
124 | /// Helper struct for tagging and [`AlgorithmConfig`] or [`Experiment`] with a name. | |
125 | #[derive(Clone, Debug, Serialize, Deserialize)] | |
126 | pub struct Named<Data> { | |
127 | pub name : String, | |
128 | #[serde(flatten)] | |
129 | pub data : Data, | |
130 | } | |
131 | ||
132 | /// 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
|
133 | #[derive(ValueEnum, Debug, Copy, Clone, Eq, PartialEq, Hash, Serialize, Deserialize)] |
0 | 134 | pub enum DefaultAlgorithm { |
135 | /// The μFB forward-backward method | |
136 | #[clap(name = "fb")] | |
137 | FB, | |
138 | /// The μFISTA inertial forward-backward method | |
139 | #[clap(name = "fista")] | |
140 | FISTA, | |
141 | /// The “fully corrective” conditional gradient method | |
142 | #[clap(name = "fw")] | |
143 | FW, | |
144 | /// The “relaxed conditional gradient method | |
145 | #[clap(name = "fwrelax")] | |
146 | FWRelax, | |
147 | /// The μPDPS primal-dual proximal splitting method | |
148 | #[clap(name = "pdps")] | |
149 | PDPS, | |
150 | } | |
151 | ||
152 | impl DefaultAlgorithm { | |
153 | /// Returns the algorithm configuration corresponding to the algorithm shorthand | |
154 | pub fn default_config<F : Float>(&self) -> AlgorithmConfig<F> { | |
155 | use DefaultAlgorithm::*; | |
156 | match *self { | |
157 | FB => AlgorithmConfig::FB(Default::default()), | |
158 | FISTA => AlgorithmConfig::FB(FBConfig{ | |
159 | meta : FBMetaAlgorithm::InertiaFISTA, | |
160 | .. Default::default() | |
161 | }), | |
162 | FW => AlgorithmConfig::FW(Default::default()), | |
163 | FWRelax => AlgorithmConfig::FW(FWConfig{ | |
164 | variant : FWVariant::Relaxed, | |
165 | .. Default::default() | |
166 | }), | |
167 | PDPS => AlgorithmConfig::PDPS(Default::default()), | |
168 | } | |
169 | } | |
170 | ||
171 | /// Returns the [`Named`] algorithm corresponding to the algorithm shorthand | |
172 | pub fn get_named<F : Float>(&self) -> Named<AlgorithmConfig<F>> { | |
173 | self.to_named(self.default_config()) | |
174 | } | |
175 | ||
176 | pub fn to_named<F : Float>(self, alg : AlgorithmConfig<F>) -> Named<AlgorithmConfig<F>> { | |
177 | let name = self.to_possible_value().unwrap().get_name().to_string(); | |
178 | Named{ name , data : alg } | |
179 | } | |
180 | } | |
181 | ||
182 | ||
183 | // // Floats cannot be hashed directly, so just hash the debug formatting | |
184 | // // for use as file identifier. | |
185 | // impl<F : Float> Hash for AlgorithmConfig<F> { | |
186 | // fn hash<H: Hasher>(&self, state: &mut H) { | |
187 | // format!("{:?}", self).hash(state); | |
188 | // } | |
189 | // } | |
190 | ||
191 | /// Plotting level configuration | |
192 | #[derive(Copy, Clone, Eq, PartialEq, Ord, PartialOrd, Serialize, ValueEnum, Debug)] | |
193 | pub enum PlotLevel { | |
194 | /// Plot nothing | |
195 | #[clap(name = "none")] | |
196 | None, | |
197 | /// Plot problem data | |
198 | #[clap(name = "data")] | |
199 | Data, | |
200 | /// Plot iterationwise state | |
201 | #[clap(name = "iter")] | |
202 | Iter, | |
203 | } | |
204 | ||
205 | type DefaultBT<F, const N : usize> = BT< | |
206 | DynamicDepth, | |
207 | F, | |
208 | usize, | |
209 | Bounds<F>, | |
210 | N | |
211 | >; | |
212 | type DefaultSeminormOp<F, K, const N : usize> = ConvolutionOp<F, K, DefaultBT<F, N>, N>; | |
213 | type DefaultSG<F, Sensor, Spread, const N : usize> = SensorGrid::< | |
214 | F, | |
215 | Sensor, | |
216 | Spread, | |
217 | DefaultBT<F, N>, | |
218 | N | |
219 | >; | |
220 | ||
221 | /// This is a dirty workaround to rust-csv not supporting struct flattening etc. | |
222 | #[derive(Serialize)] | |
223 | struct CSVLog<F> { | |
224 | iter : usize, | |
225 | cpu_time : f64, | |
226 | value : F, | |
227 | post_value : F, | |
228 | n_spikes : usize, | |
229 | inner_iters : usize, | |
230 | merged : usize, | |
231 | pruned : usize, | |
232 | this_iters : usize, | |
233 | } | |
234 | ||
235 | /// Collected experiment statistics | |
236 | #[derive(Clone, Debug, Serialize)] | |
237 | struct ExperimentStats<F : Float> { | |
238 | /// Signal-to-noise ratio in decibels | |
239 | ssnr : F, | |
240 | /// Proportion of noise in the signal as a number in $[0, 1]$. | |
241 | noise_ratio : F, | |
242 | /// When the experiment was run (UTC) | |
243 | when : DateTime<Utc>, | |
244 | } | |
245 | ||
246 | #[replace_float_literals(F::cast_from(literal))] | |
247 | impl<F : Float> ExperimentStats<F> { | |
248 | /// Calculate [`ExperimentStats`] based on a noisy `signal` and the separated `noise` signal. | |
249 | fn new<E : Euclidean<F>>(signal : &E, noise : &E) -> Self { | |
250 | let s = signal.norm2_squared(); | |
251 | let n = noise.norm2_squared(); | |
252 | let noise_ratio = (n / s).sqrt(); | |
253 | let ssnr = 10.0 * (s / n).log10(); | |
254 | ExperimentStats { | |
255 | ssnr, | |
256 | noise_ratio, | |
257 | when : Utc::now(), | |
258 | } | |
259 | } | |
260 | } | |
261 | /// Collected algorithm statistics | |
262 | #[derive(Clone, Debug, Serialize)] | |
263 | struct AlgorithmStats<F : Float> { | |
264 | /// Overall CPU time spent | |
265 | cpu_time : F, | |
266 | /// Real time spent | |
267 | elapsed : F | |
268 | } | |
269 | ||
270 | ||
271 | /// A wrapper for [`serde_json::to_writer_pretty`] that takes a filename as input | |
272 | /// and outputs a [`DynError`]. | |
273 | fn write_json<T : Serialize>(filename : String, data : &T) -> DynError { | |
274 | serde_json::to_writer_pretty(std::fs::File::create(filename)?, data)?; | |
275 | Ok(()) | |
276 | } | |
277 | ||
278 | ||
279 | /// Struct for experiment configurations | |
280 | #[derive(Debug, Clone, Serialize)] | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
281 | pub struct ExperimentV2<F, NoiseDistr, S, K, P, const N : usize> |
0 | 282 | where F : Float, |
283 | [usize; N] : Serialize, | |
284 | NoiseDistr : Distribution<F>, | |
285 | S : Sensor<F, N>, | |
286 | P : Spread<F, N>, | |
287 | K : SimpleConvolutionKernel<F, N>, | |
288 | { | |
289 | /// Domain $Ω$. | |
290 | pub domain : Cube<F, N>, | |
291 | /// Number of sensors along each dimension | |
292 | pub sensor_count : [usize; N], | |
293 | /// Noise distribution | |
294 | pub noise_distr : NoiseDistr, | |
295 | /// Seed for random noise generation (for repeatable experiments) | |
296 | pub noise_seed : u64, | |
297 | /// Sensor $θ$; $θ * ψ$ forms the forward operator $𝒜$. | |
298 | pub sensor : S, | |
299 | /// Spread $ψ$; $θ * ψ$ forms the forward operator $𝒜$. | |
300 | pub spread : P, | |
301 | /// Kernel $ρ$ of $𝒟$. | |
302 | pub kernel : K, | |
303 | /// True point sources | |
304 | pub μ_hat : DiscreteMeasure<Loc<F, N>, F>, | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
305 | /// Regularisation term and parameter |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
306 | pub regularisation : Regularisation<F>, |
0 | 307 | /// For plotting : how wide should the kernels be plotted |
308 | pub kernel_plot_width : F, | |
309 | /// Data term | |
310 | pub dataterm : DataTerm, | |
311 | /// A map of default configurations for algorithms | |
312 | #[serde(skip)] | |
313 | pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>, | |
314 | } | |
315 | ||
316 | /// Trait for runnable experiments | |
317 | 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
|
318 | /// 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
|
319 | fn runall(&self, cli : &CommandLineArgs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
320 | algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError; |
0 | 321 | |
322 | /// Return algorithm default config | |
323 | fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>) | |
324 | -> Named<AlgorithmConfig<F>>; | |
325 | } | |
326 | ||
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
327 | // *** macro boilerplate *** |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
328 | macro_rules! impl_experiment { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
329 | ($type:ident, $reg_field:ident, $reg_convert:path) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
330 | // *** macro *** |
0 | 331 | impl<F, NoiseDistr, S, K, P, const N : usize> RunnableExperiment<F> for |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
332 | Named<$type<F, NoiseDistr, S, K, P, N>> |
0 | 333 | where F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>, |
334 | [usize; N] : Serialize, | |
23
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
335 | S : Sensor<F, N> + Copy + Serialize + std::fmt::Debug, |
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
336 | P : Spread<F, N> + Copy + Serialize + std::fmt::Debug, |
0 | 337 | Convolution<S, P>: Spread<F, N> + Bounded<F> + LocalAnalysis<F, Bounds<F>, N> + Copy, |
338 | AutoConvolution<P> : BoundedBy<F, K>, | |
23
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
339 | K : SimpleConvolutionKernel<F, N> + LocalAnalysis<F, Bounds<F>, N> |
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
340 | + Copy + Serialize + std::fmt::Debug, |
0 | 341 | Cube<F, N>: P2Minimise<Loc<F, N>, F> + SetOrd, |
342 | PlotLookup : Plotting<N>, | |
343 | DefaultBT<F, N> : SensorGridBT<F, S, P, N, Depth=DynamicDepth> + BTSearch<F, N>, | |
344 | BTNodeLookup: BTNode<F, usize, Bounds<F>, N>, | |
345 | DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>, | |
23
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
346 | NoiseDistr : Distribution<F> + Serialize + std::fmt::Debug { |
0 | 347 | |
348 | fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>) | |
349 | -> Named<AlgorithmConfig<F>> { | |
350 | alg.to_named( | |
351 | self.data | |
352 | .algorithm_defaults | |
353 | .get(&alg) | |
354 | .map_or_else(|| alg.default_config(), | |
355 | |config| config.clone()) | |
356 | .cli_override(cli) | |
357 | ) | |
358 | } | |
359 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
360 | fn runall(&self, cli : &CommandLineArgs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
361 | 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
|
362 | // Get experiment configuration |
0 | 363 | let &Named { |
364 | name : ref experiment_name, | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
365 | data : $type { |
0 | 366 | domain, sensor_count, ref noise_distr, sensor, spread, kernel, |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
367 | ref μ_hat, /*regularisation,*/ kernel_plot_width, dataterm, noise_seed, |
0 | 368 | .. |
369 | } | |
370 | } = self; | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
371 | #[allow(deprecated)] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
372 | let regularisation = $reg_convert(self.data.$reg_field); |
0 | 373 | |
23
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
374 | println!("{}\n{}", |
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
375 | format!("Performing experiment {}…", experiment_name).cyan(), |
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
376 | format!("{:?}", &self.data).bright_black()); |
9869fa1e0ccd
Print out experiment information when running it
Tuomo Valkonen <tuomov@iki.fi>
parents:
20
diff
changeset
|
377 | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
378 | // Set up output directory |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
379 | 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
|
380 | |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
381 | // Set up algorithms |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
382 | let iterator_options = AlgIteratorOptions{ |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
383 | max_iter : cli.max_iter, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
384 | verbose_iter : cli.verbose_iter |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
385 | .map_or(Verbose::Logarithmic(10), |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
386 | |n| Verbose::Every(n)), |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
387 | quiet : cli.quiet, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
388 | }; |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
389 | 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
|
390 | (Some(algs), _) => algs, |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
391 | (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
|
392 | (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
|
393 | }; |
0 | 394 | |
395 | // Set up operators | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
396 | let depth = DynamicDepth(8); |
0 | 397 | let opA = DefaultSG::new(domain, sensor_count, sensor, spread, depth); |
398 | let op𝒟 = DefaultSeminormOp::new(depth, domain, kernel); | |
399 | ||
400 | // Set up random number generator. | |
401 | let mut rng = StdRng::seed_from_u64(noise_seed); | |
402 | ||
403 | // Generate the data and calculate SSNR statistic | |
404 | let b_hat = opA.apply(μ_hat); | |
405 | let noise = DVector::from_distribution(b_hat.len(), &noise_distr, &mut rng); | |
406 | let b = &b_hat + &noise; | |
407 | // Need to wrap calc_ssnr into a function to hide ultra-lame nalgebra::RealField | |
408 | // overloading log10 and conflicting with standard NumTraits one. | |
409 | let stats = ExperimentStats::new(&b, &noise); | |
410 | ||
411 | // Save experiment configuration and statistics | |
412 | let mkname_e = |t| format!("{prefix}{t}.json", prefix = prefix, t = t); | |
413 | std::fs::create_dir_all(&prefix)?; | |
414 | 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
|
415 | write_json(mkname_e("config"), cli)?; |
0 | 416 | write_json(mkname_e("stats"), &stats)?; |
417 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
418 | plotall(cli, &prefix, &domain, &sensor, &kernel, &spread, |
0 | 419 | &μ_hat, &op𝒟, &opA, &b_hat, &b, kernel_plot_width)?; |
420 | ||
421 | // Run the algorithm(s) | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
422 | for named @ Named { name : alg_name, data : alg } in algorithms.iter() { |
0 | 423 | let this_prefix = format!("{}{}/", prefix, alg_name); |
424 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
425 | let running = || if !cli.quiet { |
0 | 426 | println!("{}\n{}\n{}", |
427 | 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
|
428 | format!("{:?}", iterator_options).bright_black(), |
0 | 429 | format!("{:?}", alg).bright_black()); |
430 | }; | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
431 | let not_implemented = || { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
432 | let msg = format!("Algorithm “{alg_name}” not implemented for \ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
433 | dataterm {dataterm:?} and regularisation {regularisation:?}. \ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
434 | Skipping.").red(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
435 | eprintln!("{}", msg); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
436 | }; |
0 | 437 | // Create Logger and IteratorFactory |
438 | let mut logger = Logger::new(); | |
439 | let findim_data = prepare_optimise_weights(&opA); | |
440 | let inner_config : InnerSettings<F> = Default::default(); | |
441 | let inner_it = inner_config.iterator_options; | |
442 | let logmap = |iter, Timed { cpu_time, data }| { | |
443 | let IterInfo { | |
444 | value, | |
445 | n_spikes, | |
446 | inner_iters, | |
447 | merged, | |
448 | pruned, | |
449 | postprocessing, | |
450 | this_iters, | |
451 | .. | |
452 | } = data; | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
453 | let post_value = match (postprocessing, dataterm, regularisation) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
454 | (Some(mut μ), DataTerm::L2Squared, Regularisation::Radon(α)) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
455 | // Comparison postprocessing is only implemented for the case handled |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
456 | // by the FW variants. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
457 | optimise_weights( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
458 | &mut μ, &opA, &b, α, &findim_data, &inner_config, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
459 | inner_it |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
460 | ); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
461 | dataterm.value_at_residual(opA.apply(&μ) - &b) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
462 | + regularisation.apply(&μ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
463 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
464 | _ => value, |
0 | 465 | }; |
466 | CSVLog { | |
467 | iter, | |
468 | value, | |
469 | post_value, | |
470 | n_spikes, | |
471 | cpu_time : cpu_time.as_secs_f64(), | |
472 | inner_iters, | |
473 | merged, | |
474 | pruned, | |
475 | this_iters | |
476 | } | |
477 | }; | |
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
478 | let iterator = iterator_options.instantiate() |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
479 | .timed() |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
480 | .mapped(logmap) |
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
481 | .into_log(&mut logger); |
0 | 482 | let plotgrid = lingrid(&domain, &[if N==1 { 1000 } else { 100 }; N]); |
483 | ||
484 | // 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
|
485 | let plot_count = if cli.plot >= PlotLevel::Iter { 2000 } else { 0 }; |
0 | 486 | let plotter = SeqPlotter::new(this_prefix, plot_count, plotgrid); |
487 | ||
488 | // Run the algorithm | |
489 | let start = Instant::now(); | |
490 | let start_cpu = ProcessTime::now(); | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
491 | let μ = match alg { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
492 | AlgorithmConfig::FB(ref algconfig) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
493 | match (regularisation, dataterm) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
494 | (Regularisation::NonnegRadon(α), DataTerm::L2Squared) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
495 | running(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
496 | pointsource_fb_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
497 | &opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
498 | iterator, plotter |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
499 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
500 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
501 | (Regularisation::Radon(α), DataTerm::L2Squared) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
502 | running(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
503 | pointsource_fb_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
504 | &opA, &b, RadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
505 | iterator, plotter |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
506 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
507 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
508 | _ => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
509 | not_implemented(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
510 | continue |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
511 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
512 | } |
0 | 513 | }, |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
514 | AlgorithmConfig::PDPS(ref algconfig) => { |
0 | 515 | running(); |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
516 | match (regularisation, dataterm) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
517 | (Regularisation::NonnegRadon(α), DataTerm::L2Squared) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
518 | pointsource_pdps_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
519 | &opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
520 | iterator, plotter, L2Squared |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
521 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
522 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
523 | (Regularisation::Radon(α),DataTerm::L2Squared) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
524 | pointsource_pdps_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
525 | &opA, &b, RadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
526 | iterator, plotter, L2Squared |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
527 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
528 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
529 | (Regularisation::NonnegRadon(α), DataTerm::L1) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
530 | pointsource_pdps_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
531 | &opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
532 | iterator, plotter, L1 |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
533 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
534 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
535 | (Regularisation::Radon(α), DataTerm::L1) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
536 | pointsource_pdps_reg( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
537 | &opA, &b, RadonRegTerm(α), &op𝒟, algconfig, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
538 | iterator, plotter, L1 |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
539 | ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
540 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
541 | } |
0 | 542 | }, |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
543 | AlgorithmConfig::FW(ref algconfig) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
544 | match (regularisation, dataterm) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
545 | (Regularisation::Radon(α), DataTerm::L2Squared) => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
546 | running(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
547 | pointsource_fw(&opA, &b, α, algconfig, iterator, plotter) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
548 | }, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
549 | _ => { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
550 | not_implemented(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
551 | continue |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
552 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
553 | } |
0 | 554 | } |
555 | }; | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
556 | |
0 | 557 | let elapsed = start.elapsed().as_secs_f64(); |
558 | let cpu_time = start_cpu.elapsed().as_secs_f64(); | |
559 | ||
560 | println!("{}", format!("Elapsed {elapsed}s (CPU time {cpu_time}s)… ").yellow()); | |
561 | ||
562 | // Save results | |
563 | println!("{}", "Saving results…".green()); | |
564 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
565 | let mkname = |t| format!("{prefix}{alg_name}_{t}"); |
0 | 566 | |
567 | write_json(mkname("config.json"), &named)?; | |
568 | write_json(mkname("stats.json"), &AlgorithmStats { cpu_time, elapsed })?; | |
569 | μ.write_csv(mkname("reco.txt"))?; | |
570 | logger.write_csv(mkname("log.txt"))?; | |
571 | } | |
572 | ||
573 | Ok(()) | |
574 | } | |
575 | } | |
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
576 | // *** macro end boiler plate *** |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
577 | }} |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
578 | // *** actual code *** |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
579 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
580 | impl_experiment!(ExperimentV2, regularisation, std::convert::identity); |
0 | 581 | |
582 | /// Plot experiment setup | |
583 | #[replace_float_literals(F::cast_from(literal))] | |
584 | 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
|
585 | cli : &CommandLineArgs, |
0 | 586 | prefix : &String, |
587 | domain : &Cube<F, N>, | |
588 | sensor : &Sensor, | |
589 | kernel : &Kernel, | |
590 | spread : &Spread, | |
591 | μ_hat : &DiscreteMeasure<Loc<F, N>, F>, | |
592 | op𝒟 : &𝒟, | |
593 | opA : &A, | |
594 | b_hat : &A::Observable, | |
595 | b : &A::Observable, | |
596 | kernel_plot_width : F, | |
597 | ) -> DynError | |
598 | where F : Float + ToNalgebraRealField, | |
599 | Sensor : RealMapping<F, N> + Support<F, N> + Clone, | |
600 | Spread : RealMapping<F, N> + Support<F, N> + Clone, | |
601 | Kernel : RealMapping<F, N> + Support<F, N>, | |
602 | Convolution<Sensor, Spread> : RealMapping<F, N> + Support<F, N>, | |
603 | 𝒟 : DiscreteMeasureOp<Loc<F, N>, F>, | |
604 | 𝒟::Codomain : RealMapping<F, N>, | |
605 | A : ForwardModel<Loc<F, N>, F>, | |
606 | A::PreadjointCodomain : RealMapping<F, N> + Bounded<F>, | |
607 | PlotLookup : Plotting<N>, | |
608 | Cube<F, N> : SetOrd { | |
609 | ||
9
21b0e537ac0e
Command line parameter passing simplifications and make `-o` required.
Tuomo Valkonen <tuomov@iki.fi>
parents:
0
diff
changeset
|
610 | if cli.plot < PlotLevel::Data { |
0 | 611 | return Ok(()) |
612 | } | |
613 | ||
614 | let base = Convolution(sensor.clone(), spread.clone()); | |
615 | ||
616 | let resolution = if N==1 { 100 } else { 40 }; | |
617 | let pfx = |n| format!("{}{}", prefix, n); | |
618 | let plotgrid = lingrid(&[[-kernel_plot_width, kernel_plot_width]; N].into(), &[resolution; N]); | |
619 | ||
620 | PlotLookup::plot_into_file(sensor, plotgrid, pfx("sensor"), "sensor".to_string()); | |
621 | PlotLookup::plot_into_file(kernel, plotgrid, pfx("kernel"), "kernel".to_string()); | |
622 | PlotLookup::plot_into_file(spread, plotgrid, pfx("spread"), "spread".to_string()); | |
623 | PlotLookup::plot_into_file(&base, plotgrid, pfx("base_sensor"), "base_sensor".to_string()); | |
624 | ||
625 | let plotgrid2 = lingrid(&domain, &[resolution; N]); | |
626 | ||
627 | let ω_hat = op𝒟.apply(μ_hat); | |
628 | let noise = opA.preadjoint().apply(opA.apply(μ_hat) - b); | |
629 | PlotLookup::plot_into_file(&ω_hat, plotgrid2, pfx("omega_hat"), "ω̂".to_string()); | |
630 | PlotLookup::plot_into_file(&noise, plotgrid2, pfx("omega_noise"), | |
631 | "noise Aᵀ(Aμ̂ - b)".to_string()); | |
632 | ||
633 | let preadj_b = opA.preadjoint().apply(b); | |
634 | let preadj_b_hat = opA.preadjoint().apply(b_hat); | |
635 | //let bounds = preadj_b.bounds().common(&preadj_b_hat.bounds()); | |
636 | PlotLookup::plot_into_file_spikes( | |
637 | "Aᵀb".to_string(), &preadj_b, | |
638 | "Aᵀb̂".to_string(), Some(&preadj_b_hat), | |
639 | plotgrid2, None, &μ_hat, | |
640 | pfx("omega_b") | |
641 | ); | |
642 | ||
643 | // Save true solution and observables | |
644 | let pfx = |n| format!("{}{}", prefix, n); | |
645 | μ_hat.write_csv(pfx("orig.txt"))?; | |
646 | opA.write_observable(&b_hat, pfx("b_hat"))?; | |
647 | opA.write_observable(&b, pfx("b_noisy")) | |
648 | } | |
649 | ||
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
650 | // |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
651 | // Deprecated interface |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
652 | // |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
653 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
654 | /// Struct for experiment configurations |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
655 | #[derive(Debug, Clone, Serialize)] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
656 | pub struct Experiment<F, NoiseDistr, S, K, P, const N : usize> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
657 | where F : Float, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
658 | [usize; N] : Serialize, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
659 | NoiseDistr : Distribution<F>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
660 | S : Sensor<F, N>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
661 | P : Spread<F, N>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
662 | K : SimpleConvolutionKernel<F, N>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
663 | { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
664 | /// Domain $Ω$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
665 | pub domain : Cube<F, N>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
666 | /// Number of sensors along each dimension |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
667 | pub sensor_count : [usize; N], |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
668 | /// Noise distribution |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
669 | pub noise_distr : NoiseDistr, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
670 | /// Seed for random noise generation (for repeatable experiments) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
671 | pub noise_seed : u64, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
672 | /// Sensor $θ$; $θ * ψ$ forms the forward operator $𝒜$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
673 | pub sensor : S, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
674 | /// Spread $ψ$; $θ * ψ$ forms the forward operator $𝒜$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
675 | pub spread : P, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
676 | /// Kernel $ρ$ of $𝒟$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
677 | pub kernel : K, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
678 | /// True point sources |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
679 | pub μ_hat : DiscreteMeasure<Loc<F, N>, F>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
680 | /// Regularisation parameter |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
681 | #[deprecated(note = "Use [`ExperimentV2`], which replaces `α` by more generic `regularisation`")] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
682 | pub α : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
683 | /// For plotting : how wide should the kernels be plotted |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
684 | pub kernel_plot_width : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
685 | /// Data term |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
686 | pub dataterm : DataTerm, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
687 | /// A map of default configurations for algorithms |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
688 | #[serde(skip)] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
689 | pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
690 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
691 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
23
diff
changeset
|
692 | impl_experiment!(Experiment, α, Regularisation::NonnegRadon); |