pointsource_algs: comparison src/run.rs

-:9869fa1e0ccd
+:d29d1fcf5423
 use alg_tools::tabledump::TableDump;
 use alg_tools::sets::Cube;
 use alg_tools::mapping::RealMapping;
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::euclidean::Euclidean;
-use alg_tools::norms::{Norm, L1};
+use alg_tools::norms::L1;
 use alg_tools::lingrid::lingrid;
 use alg_tools::sets::SetOrd;
 use crate::kernels::*;
 use crate::types::*;
 use crate::measures::*;
 use crate::measures::merging::SpikeMerging;
 use crate::forward_model::*;
 use crate::fb::{
 FBConfig,
-pointsource_fb,
+pointsource_fb_reg,
-FBMetaAlgorithm, FBGenericConfig,
+FBMetaAlgorithm,
+FBGenericConfig,
 };
 use crate::pdps::{
 PDPSConfig,
 L2Squared,
-pointsource_pdps,
+pointsource_pdps_reg,
 };
 use crate::frank_wolfe::{
 FWConfig,
 FWVariant,
 pointsource_fw,
 use crate::subproblem::InnerSettings;
 use crate::seminorms::*;
 use crate::plot::*;
 use crate::{AlgorithmOverrides, CommandLineArgs};
 use crate::tolerance::Tolerance;
+use crate::regularisation::{Regularisation, RadonRegTerm, NonnegRadonRegTerm};
 /// Available algorithms and their configurations
 #[derive(Copy, Clone, Debug, Serialize, Deserialize)]
 pub enum AlgorithmConfig<F : Float> {
 FB(FBConfig<F>),
 fn write_json<T : Serialize>(filename : String, data : &T) -> DynError {
 serde_json::to_writer_pretty(std::fs::File::create(filename)?, data)?;
 Ok(())
 }
+/// Struct for experiment configurations
+#[derive(Debug, Clone, Serialize)]
+pub struct ExperimentV2<F, NoiseDistr, S, K, P, const N : usize>
+where F : Float,
+[usize; N] : Serialize,
+NoiseDistr : Distribution<F>,
+S : Sensor<F, N>,
+P : Spread<F, N>,
+K : SimpleConvolutionKernel<F, N>,
+{
+/// Domain $Ω$.
+pub domain : Cube<F, N>,
+/// Number of sensors along each dimension
+pub sensor_count : [usize; N],
+/// Noise distribution
+pub noise_distr : NoiseDistr,
+/// Seed for random noise generation (for repeatable experiments)
+pub noise_seed : u64,
+/// Sensor $θ$; $θ * ψ$ forms the forward operator $𝒜$.
+pub sensor : S,
+/// Spread $ψ$; $θ * ψ$ forms the forward operator $𝒜$.
+pub spread : P,
+/// Kernel $ρ$ of $𝒟$.
+pub kernel : K,
+/// True point sources
+pub μ_hat : DiscreteMeasure<Loc<F, N>, F>,
+/// Regularisation term and parameter
+pub regularisation : Regularisation<F>,
+/// For plotting : how wide should the kernels be plotted
+pub kernel_plot_width : F,
+/// Data term
+pub dataterm : DataTerm,
+/// A map of default configurations for algorithms
+#[serde(skip)]
+pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>,
+}
+/// Trait for runnable experiments
+pub trait RunnableExperiment<F : ClapFloat> {
+/// Run all algorithms provided, or default algorithms if none provided, on the experiment.
+fn runall(&self, cli : &CommandLineArgs,
+algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError;
+/// Return algorithm default config
+fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>)
+-> Named<AlgorithmConfig<F>>;
+}
+// *** macro boilerplate ***
+macro_rules! impl_experiment {
+($type:ident, $reg_field:ident, $reg_convert:path) => {
+// *** macro ***
+impl<F, NoiseDistr, S, K, P, const N : usize> RunnableExperiment<F> for
+Named<$type<F, NoiseDistr, S, K, P, N>>
+where F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>,
+[usize; N] : Serialize,
+S : Sensor<F, N> + Copy + Serialize + std::fmt::Debug,
+P : Spread<F, N> + Copy + Serialize + std::fmt::Debug,
+Convolution<S, P>: Spread<F, N> + Bounded<F> + LocalAnalysis<F, Bounds<F>, N> + Copy,
+AutoConvolution<P> : BoundedBy<F, K>,
+K : SimpleConvolutionKernel<F, N> + LocalAnalysis<F, Bounds<F>, N>
++ Copy + Serialize + std::fmt::Debug,
+Cube<F, N>: P2Minimise<Loc<F, N>, F> + SetOrd,
+PlotLookup : Plotting<N>,
+DefaultBT<F, N> : SensorGridBT<F, S, P, N, Depth=DynamicDepth> + BTSearch<F, N>,
+BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
+DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
+NoiseDistr : Distribution<F> + Serialize + std::fmt::Debug {
+fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>)
+-> Named<AlgorithmConfig<F>> {
+alg.to_named(
+self.data
+.algorithm_defaults
+.get(&alg)
+.map_or_else(|| alg.default_config(),
+|config| config.clone())
+.cli_override(cli)
+)
+}
+fn runall(&self, cli : &CommandLineArgs,
+algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError {
+// Get experiment configuration
+let &Named {
+name : ref experiment_name,
+data : $type {
+domain, sensor_count, ref noise_distr, sensor, spread, kernel,
+ref μ_hat, /*regularisation,*/ kernel_plot_width, dataterm, noise_seed,
+..
+}
+} = self;
+#[allow(deprecated)]
+let regularisation = $reg_convert(self.data.$reg_field);
+println!("{}\n{}",
+format!("Performing experiment {}…", experiment_name).cyan(),
+format!("{:?}", &self.data).bright_black());
+// Set up output directory
+let prefix = format!("{}/{}/", cli.outdir, self.name);
+// Set up algorithms
+let iterator_options = AlgIteratorOptions{
+max_iter : cli.max_iter,
+verbose_iter : cli.verbose_iter
+.map_or(Verbose::Logarithmic(10),
+|n| Verbose::Every(n)),
+quiet : cli.quiet,
+};
+let algorithms = match (algs, self.data.dataterm) {
+(Some(algs), _) => algs,
+(None, DataTerm::L2Squared) => vec![DefaultAlgorithm::FB.get_named()],
+(None, DataTerm::L1) => vec![DefaultAlgorithm::PDPS.get_named()],
+};
+// Set up operators
+let depth = DynamicDepth(8);
+let opA = DefaultSG::new(domain, sensor_count, sensor, spread, depth);
+let op𝒟 = DefaultSeminormOp::new(depth, domain, kernel);
+// Set up random number generator.
+let mut rng = StdRng::seed_from_u64(noise_seed);
+// Generate the data and calculate SSNR statistic
+let b_hat = opA.apply(μ_hat);
+let noise = DVector::from_distribution(b_hat.len(), &noise_distr, &mut rng);
+let b = &b_hat + &noise;
+// Need to wrap calc_ssnr into a function to hide ultra-lame nalgebra::RealField
+// overloading log10 and conflicting with standard NumTraits one.
+let stats = ExperimentStats::new(&b, &noise);
+// Save experiment configuration and statistics
+let mkname_e = |t| format!("{prefix}{t}.json", prefix = prefix, t = t);
+std::fs::create_dir_all(&prefix)?;
+write_json(mkname_e("experiment"), self)?;
+write_json(mkname_e("config"), cli)?;
+write_json(mkname_e("stats"), &stats)?;
+plotall(cli, &prefix, &domain, &sensor, &kernel, &spread,
+&μ_hat, &op𝒟, &opA, &b_hat, &b, kernel_plot_width)?;
+// Run the algorithm(s)
+for named @ Named { name : alg_name, data : alg } in algorithms.iter() {
+let this_prefix = format!("{}{}/", prefix, alg_name);
+let running = || if !cli.quiet {
+println!("{}\n{}\n{}",
+format!("Running {} on experiment {}…", alg_name, experiment_name).cyan(),
+format!("{:?}", iterator_options).bright_black(),
+format!("{:?}", alg).bright_black());
+};
+let not_implemented = || {
+let msg = format!("Algorithm “{alg_name}” not implemented for \
+dataterm {dataterm:?} and regularisation {regularisation:?}. \
+Skipping.").red();
+eprintln!("{}", msg);
+};
+// Create Logger and IteratorFactory
+let mut logger = Logger::new();
+let findim_data = prepare_optimise_weights(&opA);
+let inner_config : InnerSettings<F> = Default::default();
+let inner_it = inner_config.iterator_options;
+let logmap = |iter, Timed { cpu_time, data }| {
+let IterInfo {
+value,
+n_spikes,
+inner_iters,
+merged,
+pruned,
+postprocessing,
+this_iters,
+..
+} = data;
+let post_value = match (postprocessing, dataterm, regularisation) {
+(Some(mut μ), DataTerm::L2Squared, Regularisation::Radon(α)) => {
+// Comparison postprocessing is only implemented for the case handled
+// by the FW variants.
+optimise_weights(
+&mut μ, &opA, &b, α, &findim_data, &inner_config,
+inner_it
+);
+dataterm.value_at_residual(opA.apply(&μ) - &b)
++ regularisation.apply(&μ)
+},
+_ => value,
+};
+CSVLog {
+iter,
+value,
+post_value,
+n_spikes,
+cpu_time : cpu_time.as_secs_f64(),
+inner_iters,
+merged,
+pruned,
+this_iters
+}
+};
+let iterator = iterator_options.instantiate()
+.timed()
+.mapped(logmap)
+.into_log(&mut logger);
+let plotgrid = lingrid(&domain, &[if N==1 { 1000 } else { 100 }; N]);
+// Create plotter and directory if needed.
+let plot_count = if cli.plot >= PlotLevel::Iter { 2000 } else { 0 };
+let plotter = SeqPlotter::new(this_prefix, plot_count, plotgrid);
+// Run the algorithm
+let start = Instant::now();
+let start_cpu = ProcessTime::now();
+let μ = match alg {
+AlgorithmConfig::FB(ref algconfig) => {
+match (regularisation, dataterm) {
+(Regularisation::NonnegRadon(α), DataTerm::L2Squared) => {
+running();
+pointsource_fb_reg(
+&opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter
+)
+},
+(Regularisation::Radon(α), DataTerm::L2Squared) => {
+running();
+pointsource_fb_reg(
+&opA, &b, RadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter
+)
+},
+_ => {
+not_implemented();
+continue
+}
+}
+},
+AlgorithmConfig::PDPS(ref algconfig) => {
+running();
+match (regularisation, dataterm) {
+(Regularisation::NonnegRadon(α), DataTerm::L2Squared) => {
+pointsource_pdps_reg(
+&opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter, L2Squared
+)
+},
+(Regularisation::Radon(α),DataTerm::L2Squared) => {
+pointsource_pdps_reg(
+&opA, &b, RadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter, L2Squared
+)
+},
+(Regularisation::NonnegRadon(α), DataTerm::L1) => {
+pointsource_pdps_reg(
+&opA, &b, NonnegRadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter, L1
+)
+},
+(Regularisation::Radon(α), DataTerm::L1) => {
+pointsource_pdps_reg(
+&opA, &b, RadonRegTerm(α), &op𝒟, algconfig,
+iterator, plotter, L1
+)
+},
+}
+},
+AlgorithmConfig::FW(ref algconfig) => {
+match (regularisation, dataterm) {
+(Regularisation::Radon(α), DataTerm::L2Squared) => {
+running();
+pointsource_fw(&opA, &b, α, algconfig, iterator, plotter)
+},
+_ => {
+not_implemented();
+continue
+}
+}
+}
+};
+let elapsed = start.elapsed().as_secs_f64();
+let cpu_time = start_cpu.elapsed().as_secs_f64();
+println!("{}", format!("Elapsed {elapsed}s (CPU time {cpu_time}s)… ").yellow());
+// Save results
+println!("{}", "Saving results…".green());
+let mkname = |t| format!("{prefix}{alg_name}_{t}");
+write_json(mkname("config.json"), &named)?;
+write_json(mkname("stats.json"), &AlgorithmStats { cpu_time, elapsed })?;
+μ.write_csv(mkname("reco.txt"))?;
+logger.write_csv(mkname("log.txt"))?;
+}
+Ok(())
+}
+}
+// *** macro end boiler plate ***
+}}
+// *** actual code ***
+impl_experiment!(ExperimentV2, regularisation, std::convert::identity);
+/// Plot experiment setup
+#[replace_float_literals(F::cast_from(literal))]
+fn plotall<F, Sensor, Kernel, Spread, 𝒟, A, const N : usize>(
+cli : &CommandLineArgs,
+prefix : &String,
+domain : &Cube<F, N>,
+sensor : &Sensor,
+kernel : &Kernel,
+spread : &Spread,
+μ_hat : &DiscreteMeasure<Loc<F, N>, F>,
+op𝒟 : &𝒟,
+opA : &A,
+b_hat : &A::Observable,
+b : &A::Observable,
+kernel_plot_width : F,
+) -> DynError
+where F : Float + ToNalgebraRealField,
+Sensor : RealMapping<F, N> + Support<F, N> + Clone,
+Spread : RealMapping<F, N> + Support<F, N> + Clone,
+Kernel : RealMapping<F, N> + Support<F, N>,
+Convolution<Sensor, Spread> : RealMapping<F, N> + Support<F, N>,
+𝒟 : DiscreteMeasureOp<Loc<F, N>, F>,
+𝒟::Codomain : RealMapping<F, N>,
+A : ForwardModel<Loc<F, N>, F>,
+A::PreadjointCodomain : RealMapping<F, N> + Bounded<F>,
+PlotLookup : Plotting<N>,
+Cube<F, N> : SetOrd {
+if cli.plot < PlotLevel::Data {
+return Ok(())
+}
+let base = Convolution(sensor.clone(), spread.clone());
+let resolution = if N==1 { 100 } else { 40 };
+let pfx = |n| format!("{}{}", prefix, n);
+let plotgrid = lingrid(&[[-kernel_plot_width, kernel_plot_width]; N].into(), &[resolution; N]);
+PlotLookup::plot_into_file(sensor, plotgrid, pfx("sensor"), "sensor".to_string());
+PlotLookup::plot_into_file(kernel, plotgrid, pfx("kernel"), "kernel".to_string());
+PlotLookup::plot_into_file(spread, plotgrid, pfx("spread"), "spread".to_string());
+PlotLookup::plot_into_file(&base, plotgrid, pfx("base_sensor"), "base_sensor".to_string());
+let plotgrid2 = lingrid(&domain, &[resolution; N]);
+let ω_hat = op𝒟.apply(μ_hat);
+let noise =  opA.preadjoint().apply(opA.apply(μ_hat) - b);
+PlotLookup::plot_into_file(&ω_hat, plotgrid2, pfx("omega_hat"), "ω̂".to_string());
+PlotLookup::plot_into_file(&noise, plotgrid2, pfx("omega_noise"),
+"noise Aᵀ(Aμ̂ - b)".to_string());
+let preadj_b =  opA.preadjoint().apply(b);
+let preadj_b_hat =  opA.preadjoint().apply(b_hat);
+//let bounds = preadj_b.bounds().common(&preadj_b_hat.bounds());
+PlotLookup::plot_into_file_spikes(
+"Aᵀb".to_string(), &preadj_b,
+"Aᵀb̂".to_string(), Some(&preadj_b_hat),
+plotgrid2, None, &μ_hat,
+pfx("omega_b")
+);
+// Save true solution and observables
+let pfx = |n| format!("{}{}", prefix, n);
+μ_hat.write_csv(pfx("orig.txt"))?;
+opA.write_observable(&b_hat, pfx("b_hat"))?;
+opA.write_observable(&b, pfx("b_noisy"))
+}
+//
+// Deprecated interface
+//
 /// Struct for experiment configurations
 #[derive(Debug, Clone, Serialize)]
 pub struct Experiment<F, NoiseDistr, S, K, P, const N : usize>
 where F : Float,
 /// Kernel $ρ$ of $𝒟$.
 pub kernel : K,
 /// True point sources
 pub μ_hat : DiscreteMeasure<Loc<F, N>, F>,
 /// Regularisation parameter
+#[deprecated(note = "Use [`ExperimentV2`], which replaces `α` by more generic `regularisation`")]
 pub α : F,
 /// For plotting : how wide should the kernels be plotted
 pub kernel_plot_width : F,
 /// Data term
 pub dataterm : DataTerm,
 /// A map of default configurations for algorithms
 #[serde(skip)]
 pub algorithm_defaults : HashMap<DefaultAlgorithm, AlgorithmConfig<F>>,
 }
-/// Trait for runnable experiments
+impl_experiment!(Experiment, α, Regularisation::NonnegRadon);
-pub trait RunnableExperiment<F : ClapFloat> {
-/// Run all algorithms provided, or default algorithms if none provided, on the experiment.
-fn runall(&self, cli : &CommandLineArgs,
-algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError;
-/// Return algorithm default config
-fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>)
--> Named<AlgorithmConfig<F>>;
-}
-impl<F, NoiseDistr, S, K, P, const N : usize> RunnableExperiment<F> for
-Named<Experiment<F, NoiseDistr, S, K, P, N>>
-where F : ClapFloat + nalgebra::RealField + ToNalgebraRealField<MixedType=F>,
-[usize; N] : Serialize,
-S : Sensor<F, N> + Copy + Serialize + std::fmt::Debug,
-P : Spread<F, N> + Copy + Serialize + std::fmt::Debug,
-Convolution<S, P>: Spread<F, N> + Bounded<F> + LocalAnalysis<F, Bounds<F>, N> + Copy,
-AutoConvolution<P> : BoundedBy<F, K>,
-K : SimpleConvolutionKernel<F, N> + LocalAnalysis<F, Bounds<F>, N>
-+ Copy + Serialize + std::fmt::Debug,
-Cube<F, N>: P2Minimise<Loc<F, N>, F> + SetOrd,
-PlotLookup : Plotting<N>,
-DefaultBT<F, N> : SensorGridBT<F, S, P, N, Depth=DynamicDepth> + BTSearch<F, N>,
-BTNodeLookup: BTNode<F, usize, Bounds<F>, N>,
-DiscreteMeasure<Loc<F, N>, F> : SpikeMerging<F>,
-NoiseDistr : Distribution<F> + Serialize + std::fmt::Debug {
-fn algorithm_defaults(&self, alg : DefaultAlgorithm, cli : &AlgorithmOverrides<F>)
--> Named<AlgorithmConfig<F>> {
-alg.to_named(
-self.data
-.algorithm_defaults
-.get(&alg)
-.map_or_else(|| alg.default_config(),
-|config| config.clone())
-.cli_override(cli)
-)
-}
-fn runall(&self, cli : &CommandLineArgs,
-algs : Option<Vec<Named<AlgorithmConfig<F>>>>) -> DynError {
-// Get experiment configuration
-let &Named {
-name : ref experiment_name,
-data : Experiment {
-domain, sensor_count, ref noise_distr, sensor, spread, kernel,
-ref μ_hat, α, kernel_plot_width, dataterm, noise_seed,
-..
-}
-} = self;
-println!("{}\n{}",
-format!("Performing experiment {}…", experiment_name).cyan(),
-format!("{:?}", &self.data).bright_black());
-// Set up output directory
-let prefix = format!("{}/{}/", cli.outdir, self.name);
-// Set up algorithms
-let iterator_options = AlgIteratorOptions{
-max_iter : cli.max_iter,
-verbose_iter : cli.verbose_iter
-.map_or(Verbose::Logarithmic(10),
-|n| Verbose::Every(n)),
-quiet : cli.quiet,
-};
-let algorithms = match (algs, self.data.dataterm) {
-(Some(algs), _) => algs,
-(None, DataTerm::L2Squared) => vec![DefaultAlgorithm::FB.get_named()],
-(None, DataTerm::L1) => vec![DefaultAlgorithm::PDPS.get_named()],
-};
-// Set up operators
-let depth = DynamicDepth(8);
-let opA = DefaultSG::new(domain, sensor_count, sensor, spread, depth);
-let op𝒟 = DefaultSeminormOp::new(depth, domain, kernel);
-// Set up random number generator.
-let mut rng = StdRng::seed_from_u64(noise_seed);
-// Generate the data and calculate SSNR statistic
-let b_hat = opA.apply(μ_hat);
-let noise = DVector::from_distribution(b_hat.len(), &noise_distr, &mut rng);
-let b = &b_hat + &noise;
-// Need to wrap calc_ssnr into a function to hide ultra-lame nalgebra::RealField
-// overloading log10 and conflicting with standard NumTraits one.
-let stats = ExperimentStats::new(&b, &noise);
-// Save experiment configuration and statistics
-let mkname_e = |t| format!("{prefix}{t}.json", prefix = prefix, t = t);
-std::fs::create_dir_all(&prefix)?;
-write_json(mkname_e("experiment"), self)?;
-write_json(mkname_e("config"), cli)?;
-write_json(mkname_e("stats"), &stats)?;
-plotall(cli, &prefix, &domain, &sensor, &kernel, &spread,
-&μ_hat, &op𝒟, &opA, &b_hat, &b, kernel_plot_width)?;
-// Run the algorithm(s)
-for named @ Named { name : alg_name, data : alg } in algorithms.iter() {
-let this_prefix = format!("{}{}/", prefix, alg_name);
-let running = || if !cli.quiet {
-println!("{}\n{}\n{}",
-format!("Running {} on experiment {}…", alg_name, experiment_name).cyan(),
-format!("{:?}", iterator_options).bright_black(),
-format!("{:?}", alg).bright_black());
-};
-// Create Logger and IteratorFactory
-let mut logger = Logger::new();
-let findim_data = prepare_optimise_weights(&opA);
-let inner_config : InnerSettings<F> = Default::default();
-let inner_it = inner_config.iterator_options;
-let logmap = |iter, Timed { cpu_time, data }| {
-let IterInfo {
-value,
-n_spikes,
-inner_iters,
-merged,
-pruned,
-postprocessing,
-this_iters,
-..
-} = data;
-let post_value = match postprocessing {
-None => value,
-Some(mut μ) => {
-match dataterm {
-DataTerm::L2Squared => {
-optimise_weights(
-&mut μ, &opA, &b, α, &findim_data, &inner_config,
-inner_it
-);
-dataterm.value_at_residual(opA.apply(&μ) - &b) + α * μ.norm(Radon)
-},
-_ => value,
-}
-}
-};
-CSVLog {
-iter,
-value,
-post_value,
-n_spikes,
-cpu_time : cpu_time.as_secs_f64(),
-inner_iters,
-merged,
-pruned,
-this_iters
-}
-};
-let iterator = iterator_options.instantiate()
-.timed()
-.mapped(logmap)
-.into_log(&mut logger);
-let plotgrid = lingrid(&domain, &[if N==1 { 1000 } else { 100 }; N]);
-// Create plotter and directory if needed.
-let plot_count = if cli.plot >= PlotLevel::Iter { 2000 } else { 0 };
-let plotter = SeqPlotter::new(this_prefix, plot_count, plotgrid);
-// Run the algorithm
-let start = Instant::now();
-let start_cpu = ProcessTime::now();
-let μ : DiscreteMeasure<Loc<F, N>, F> = match (alg, dataterm) {
-(AlgorithmConfig::FB(ref algconfig), DataTerm::L2Squared) => {
-running();
-pointsource_fb(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter)
-},
-(AlgorithmConfig::FW(ref algconfig), DataTerm::L2Squared) => {
-running();
-pointsource_fw(&opA, &b, α, &algconfig, iterator, plotter)
-},
-(AlgorithmConfig::PDPS(ref algconfig), DataTerm::L2Squared) => {
-running();
-pointsource_pdps(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter, L2Squared)
-},
-(AlgorithmConfig::PDPS(ref algconfig), DataTerm::L1) => {
-running();
-pointsource_pdps(&opA, &b, α, &op𝒟, &algconfig, iterator, plotter, L1)
-},
-_ =>  {
-let msg = format!("Algorithm “{alg_name}” not implemented for \
-dataterm {dataterm:?}. Skipping.").red();
-eprintln!("{}", msg);
-continue
-}
-};
-let elapsed = start.elapsed().as_secs_f64();
-let cpu_time = start_cpu.elapsed().as_secs_f64();
-println!("{}", format!("Elapsed {elapsed}s (CPU time {cpu_time}s)… ").yellow());
-// Save results
-println!("{}", "Saving results…".green());
-let mkname = |t| format!("{prefix}{alg_name}_{t}");
-write_json(mkname("config.json"), &named)?;
-write_json(mkname("stats.json"), &AlgorithmStats { cpu_time, elapsed })?;
-μ.write_csv(mkname("reco.txt"))?;
-logger.write_csv(mkname("log.txt"))?;
-}
-Ok(())
-}
-}
-/// Plot experiment setup
-#[replace_float_literals(F::cast_from(literal))]
-fn plotall<F, Sensor, Kernel, Spread, 𝒟, A, const N : usize>(
-cli : &CommandLineArgs,
-prefix : &String,
-domain : &Cube<F, N>,
-sensor : &Sensor,
-kernel : &Kernel,
-spread : &Spread,
-μ_hat : &DiscreteMeasure<Loc<F, N>, F>,
-op𝒟 : &𝒟,
-opA : &A,
-b_hat : &A::Observable,
-b : &A::Observable,
-kernel_plot_width : F,
-) -> DynError
-where F : Float + ToNalgebraRealField,
-Sensor : RealMapping<F, N> + Support<F, N> + Clone,
-Spread : RealMapping<F, N> + Support<F, N> + Clone,
-Kernel : RealMapping<F, N> + Support<F, N>,
-Convolution<Sensor, Spread> : RealMapping<F, N> + Support<F, N>,
-𝒟 : DiscreteMeasureOp<Loc<F, N>, F>,
-𝒟::Codomain : RealMapping<F, N>,
-A : ForwardModel<Loc<F, N>, F>,
-A::PreadjointCodomain : RealMapping<F, N> + Bounded<F>,
-PlotLookup : Plotting<N>,
-Cube<F, N> : SetOrd {
-if cli.plot < PlotLevel::Data {
-return Ok(())
-}
-let base = Convolution(sensor.clone(), spread.clone());
-let resolution = if N==1 { 100 } else { 40 };
-let pfx = |n| format!("{}{}", prefix, n);
-let plotgrid = lingrid(&[[-kernel_plot_width, kernel_plot_width]; N].into(), &[resolution; N]);
-PlotLookup::plot_into_file(sensor, plotgrid, pfx("sensor"), "sensor".to_string());
-PlotLookup::plot_into_file(kernel, plotgrid, pfx("kernel"), "kernel".to_string());
-PlotLookup::plot_into_file(spread, plotgrid, pfx("spread"), "spread".to_string());
-PlotLookup::plot_into_file(&base, plotgrid, pfx("base_sensor"), "base_sensor".to_string());
-let plotgrid2 = lingrid(&domain, &[resolution; N]);
-let ω_hat = op𝒟.apply(μ_hat);
-let noise =  opA.preadjoint().apply(opA.apply(μ_hat) - b);
-PlotLookup::plot_into_file(&ω_hat, plotgrid2, pfx("omega_hat"), "ω̂".to_string());
-PlotLookup::plot_into_file(&noise, plotgrid2, pfx("omega_noise"),
-"noise Aᵀ(Aμ̂ - b)".to_string());
-let preadj_b =  opA.preadjoint().apply(b);
-let preadj_b_hat =  opA.preadjoint().apply(b_hat);
-//let bounds = preadj_b.bounds().common(&preadj_b_hat.bounds());
-PlotLookup::plot_into_file_spikes(
-"Aᵀb".to_string(), &preadj_b,
-"Aᵀb̂".to_string(), Some(&preadj_b_hat),
-plotgrid2, None, &μ_hat,
-pfx("omega_b")
-);
-// Save true solution and observables
-let pfx = |n| format!("{}{}", prefix, n);
-μ_hat.write_csv(pfx("orig.txt"))?;
-opA.write_observable(&b_hat, pfx("b_hat"))?;
-opA.write_observable(&b, pfx("b_noisy"))
-}

Mercurial > repos > pointsource_algs / file comparison

comparison: src/run.rs

src/run.rs