--- a/src/experiments.rs Mon Jun 29 00:20:42 2026 +0300
+++ b/src/experiments.rs Sun Oct 04 13:02:26 2026 -0500
@@ -9,20 +9,23 @@
use alg_tools::direct_product::Pair;
use alg_tools::error::{DynError, DynResult};
use alg_tools::loc::Loc;
+use alg_tools::sets::Cube;
use anyhow::bail;
use clap::Parser;
use log::debug;
+use pointsource_algs::domain::{AdmissibleDomain, FullSpace, UniformRandom};
use pointsource_algs::measures::DiscreteMeasure;
use pointsource_algs::prox_penalty::RadonSquared;
use pointsource_algs::regularisation::{NonnegRadonRegTerm, RadonRegTerm, Regularisation};
use pointsource_algs::run::{
- run_fb, run_fb_pair, AlgorithmConfig, DefaultAlgorithm, Named, RunError::NotImplemented,
- RunnableExperiment, RunnableExperimentExtras,
+ run_fb_dom, run_fb_pair_dom, DefaultAlgorithm, RunError::NotImplemented, RunnableExperiment,
+ RunnableExperimentExtras,
};
use pointsource_algs::{AlgorithmOverrides, CommandLineArgs, ExperimentSetup};
use pyo3::ffi::c_str;
use pyo3::prelude::*;
use pyo3::types::{PyDict, PyFunction};
+use rand::prelude::{SeedableRng, StdRng};
use serde::{Deserialize, Serialize};
use serde_json;
use serde_with::skip_serializing_none;
@@ -264,12 +267,20 @@
/// Plotter
#[pyo3(set, get)]
plot_factory: Option<Py<PyAny>>,
+
+ /// Domain
+ #[pyo3(set, get)]
+ domain: Option<Py<PyAny>>,
+
+ /// Randon number generator seed
+ #[pyo3(set, get)]
+ seed: Option<u64>,
}
#[pymethods]
impl Problem {
#[new]
- #[pyo3(signature = (dataterm, regterm, auxterm = None, auxinit = None, μinit = None, plot_factory=None))]
+ #[pyo3(signature = (dataterm, regterm, auxterm = None, auxinit = None, μinit = None, plot_factory=None, domain = None, seed = None))]
fn new(
dataterm: Py<PyAny>,
regterm: RegTermPy,
@@ -277,8 +288,10 @@
auxinit: Option<Py<PyAny>>,
μinit: Option<DiscreteMeasure<Loc<2, f64>, f64>>,
plot_factory: Option<Py<PyAny>>,
+ domain: Option<Py<PyAny>>,
+ seed: Option<u64>,
) -> Self {
- Self { dataterm, regterm, auxterm, auxinit, μinit, plot_factory }
+ Self { dataterm, regterm, auxterm, auxinit, μinit, plot_factory, seed, domain }
}
#[getter]
@@ -287,6 +300,180 @@
}
}
+impl PythonExperiment {
+ fn do_runall_phase2<D: AdmissibleDomain<Loc<2, f64>>>(
+ &self,
+ cli: &CommandLineArgs,
+ cli_algorithm_overrides: &AlgorithmOverrides<f64>,
+ prefix: String,
+ problem: PyRef<Problem>,
+ domain: &mut D,
+ py: Python,
+ ) -> DynError {
+ // Set up algorithms
+ let algorithms = self.collect_algs(
+ cli,
+ cli_algorithm_overrides,
+ |alg| self.algorithm_overrides.get(&alg).cloned(),
+ Some(&problem.regterm.into()),
+ &[DefaultAlgorithm::RadonFB, DefaultAlgorithm::RadonSlidingFB],
+ )?;
+
+ let save_extra = |_, ()| Ok(());
+ let μ0 = problem.μinit.clone();
+
+ match (problem.regterm, &problem.auxterm, &problem.auxinit) {
+ (regterm, None, None) => {
+ debug!("… Extracting data term.");
+ let dataterm: exp_f64_2::DataTerm<'_> = problem.dataterm.extract(py)?;
+ debug!("… Extracting plotter.");
+ let plot_factory: exp_f64_2::PlotFactory<'_> =
+ if let Some(ref p) = problem.plot_factory {
+ p.extract(py)?
+ } else {
+ PythonPlotFactory::dummy()
+ };
+ let make_plotter = |prefix| plot_factory.produce(prefix);
+
+ debug!("Running.");
+ self.do_runall(
+ &prefix,
+ cli,
+ algorithms,
+ make_plotter,
+ save_extra,
+ μ0,
+ |p| match regterm {
+ RegTermPy::Radon(α) => run_fb_dom(
+ &dataterm,
+ &RadonRegTerm(α),
+ &RadonSquared,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|μ| (μ, ())),
+ RegTermPy::NonnegRadon(α) => run_fb_dom(
+ &dataterm,
+ &NonnegRadonRegTerm(α),
+ &RadonSquared,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|μ| (μ, ())),
+ },
+ )
+ }
+ (regterm, Some(auxterm), Some(z_)) => {
+ debug!("… Extracting auxiliary variable.");
+ let z0: PyResult<exp_f64_2_auxvar::AuxVar<'_>> = z_.extract(py);
+
+ match z0 {
+ Ok(z) => {
+ debug!("… Extracting data term.");
+ let dataterm: exp_f64_2_auxvar::DataTerm<'_> =
+ problem.dataterm.extract(py)?;
+ debug!("… Extracting plotter.");
+ let plot_factory: exp_f64_2::PlotFactory<'_> =
+ if let Some(ref p) = problem.plot_factory {
+ p.extract(py)?
+ } else {
+ PythonPlotFactory::dummy()
+ };
+ let make_plotter = |prefix| plot_factory.produce(prefix);
+ debug!("… Extracting auxiliary term.");
+ let auxterm: exp_f64_2_auxvar::AuxTerm<'_> = auxterm.extract(py)?;
+
+ debug!("Running.");
+ self.do_runall(
+ &prefix,
+ cli,
+ algorithms,
+ make_plotter,
+ save_extra,
+ (μ0, z),
+ |p| match regterm {
+ RegTermPy::Radon(α) => run_fb_pair_dom(
+ &dataterm,
+ &RadonRegTerm(α),
+ &RadonSquared,
+ &auxterm,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|Pair(μ, _)| (μ, ())),
+ RegTermPy::NonnegRadon(α) => run_fb_pair_dom(
+ &dataterm,
+ &NonnegRadonRegTerm(α),
+ &RadonSquared,
+ &auxterm,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|Pair(μ, _)| (μ, ())),
+ },
+ )
+ }
+ Err(_) => {
+ let z: exp_f64_2_auxvar_scalar::AuxVar<'_> = z_.extract(py)?;
+ debug!("… Extracting data term.");
+ let dataterm: exp_f64_2_auxvar_scalar::DataTerm<'_> =
+ problem.dataterm.extract(py)?;
+ debug!("… Extracting plotter.");
+ let plot_factory: exp_f64_2::PlotFactory<'_> =
+ if let Some(ref p) = problem.plot_factory {
+ p.extract(py)?
+ } else {
+ PythonPlotFactory::dummy()
+ };
+ let make_plotter = |prefix| plot_factory.produce(prefix);
+ debug!("… Extracting auxiliary term.");
+ let auxterm: exp_f64_2_auxvar_scalar::AuxTerm<'_> = auxterm.extract(py)?;
+
+ debug!("Running.");
+ self.do_runall(
+ &prefix,
+ cli,
+ algorithms,
+ make_plotter,
+ save_extra,
+ (μ0, z),
+ |p| match regterm {
+ RegTermPy::Radon(α) => run_fb_pair_dom(
+ &dataterm,
+ &RadonRegTerm(α),
+ &RadonSquared,
+ &auxterm,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|Pair(μ, _)| (μ, ())),
+ RegTermPy::NonnegRadon(α) => run_fb_pair_dom(
+ &dataterm,
+ &NonnegRadonRegTerm(α),
+ &RadonSquared,
+ &auxterm,
+ p,
+ domain,
+ |_| Err(NotImplemented.into()),
+ )
+ .map(|Pair(μ, _)| (μ, ())),
+ },
+ )
+ }
+ }
+ }
+ (_, _, _) => {
+ bail!("Errors in problem {} setup: auxiliary term or auxiliary variable initialisation given without the other", self.name);
+ }
+ }
+ }
+}
+
impl RunnableExperiment<f64> for PythonExperiment {
fn name(&self) -> &str {
&self.name
@@ -295,11 +482,8 @@
fn runall(
&self,
cli: &CommandLineArgs,
- algs: Option<Vec<Named<AlgorithmConfig<f64>>>>,
+ cli_algorithm_overrides: &AlgorithmOverrides<f64>,
) -> DynError {
- // Set up algorithms
- let algorithms = algs.unwrap_or_else(|| vec![DefaultAlgorithm::FB.get_named()]);
-
debug!(
"Entered PythonExperiment::runall for experimen {}. Attaching Python.",
self.name
@@ -310,6 +494,7 @@
let prefix = self.start(cli)?;
+ // Get the problem description
let problem: PyRef<Problem> = process_error(
"PythonExperiment::runall",
py,
@@ -318,158 +503,28 @@
.and_then(|r| r.extract(py).map_err(PyErr::from)),
)?;
- let save_extra = |_, ()| Ok(());
- let μ0 = problem.μinit.clone();
-
- match (problem.regterm, &problem.auxterm, &problem.auxinit) {
- (regterm, None, None) => {
- debug!("… Extracting data term.");
- let dataterm: exp_f64_2::DataTerm<'_> = problem.dataterm.extract(py)?;
- debug!("… Extracting plotter.");
- let plot_factory: exp_f64_2::PlotFactory<'_> =
- if let Some(ref p) = problem.plot_factory {
- p.extract(py)?
- } else {
- PythonPlotFactory::dummy()
- };
- let make_plotter = |prefix| plot_factory.produce(prefix);
-
- debug!("Running.");
- self.do_runall(
- &prefix,
- cli,
- algorithms,
- make_plotter,
- save_extra,
- μ0,
- |p| match regterm {
- RegTermPy::Radon(α) => {
- run_fb(&dataterm, &RadonRegTerm(α), &RadonSquared, p, |_| {
- Err(NotImplemented.into())
- })
- .map(|μ| (μ, ()))
- }
- RegTermPy::NonnegRadon(α) => {
- run_fb(&dataterm, &NonnegRadonRegTerm(α), &RadonSquared, p, |_| {
- Err(NotImplemented.into())
- })
- .map(|μ| (μ, ()))
- }
- },
- )
- }
- (regterm, Some(auxterm), Some(z_)) => {
- debug!("… Extracting auxiliary variable.");
- let z0: PyResult<exp_f64_2_auxvar::AuxVar<'_>> = z_.extract(py);
- match z0 {
- Ok(z) => {
- debug!("… Extracting data term.");
- let dataterm: exp_f64_2_auxvar::DataTerm<'_> =
- problem.dataterm.extract(py)?;
- debug!("… Extracting plotter.");
- let plot_factory: exp_f64_2::PlotFactory<'_> =
- if let Some(ref p) = problem.plot_factory {
- p.extract(py)?
- } else {
- PythonPlotFactory::dummy()
- };
- let make_plotter = |prefix| plot_factory.produce(prefix);
- debug!("… Extracting auxiliary term.");
- let auxterm: exp_f64_2_auxvar::AuxTerm<'_> = auxterm.extract(py)?;
-
- debug!("Running.");
- self.do_runall(
- &prefix,
- cli,
- algorithms,
- make_plotter,
- save_extra,
- (μ0, z),
- |p| match regterm {
- RegTermPy::Radon(α) => run_fb_pair(
- &dataterm,
- &RadonRegTerm(α),
- &RadonSquared,
- &auxterm,
- p,
- |_| Err(NotImplemented.into()),
- )
- .map(|Pair(μ, _)| (μ, ())),
- RegTermPy::NonnegRadon(α) => run_fb_pair(
- &dataterm,
- &NonnegRadonRegTerm(α),
- &RadonSquared,
- &auxterm,
- p,
- |_| Err(NotImplemented.into()),
- )
- .map(|Pair(μ, _)| (μ, ())),
- },
- )
- }
- Err(_) => {
- let z: exp_f64_2_auxvar_scalar::AuxVar<'_> = z_.extract(py)?;
- debug!("… Extracting data term.");
- let dataterm: exp_f64_2_auxvar_scalar::DataTerm<'_> =
- problem.dataterm.extract(py)?;
- debug!("… Extracting plotter.");
- let plot_factory: exp_f64_2::PlotFactory<'_> =
- if let Some(ref p) = problem.plot_factory {
- p.extract(py)?
- } else {
- PythonPlotFactory::dummy()
- };
- let make_plotter = |prefix| plot_factory.produce(prefix);
- debug!("… Extracting auxiliary term.");
- let auxterm: exp_f64_2_auxvar_scalar::AuxTerm<'_> =
- auxterm.extract(py)?;
-
- debug!("Running.");
- self.do_runall(
- &prefix,
- cli,
- algorithms,
- make_plotter,
- save_extra,
- (μ0, z),
- |p| match regterm {
- RegTermPy::Radon(α) => run_fb_pair(
- &dataterm,
- &RadonRegTerm(α),
- &RadonSquared,
- &auxterm,
- p,
- |_| Err(NotImplemented.into()),
- )
- .map(|Pair(μ, _)| (μ, ())),
- RegTermPy::NonnegRadon(α) => run_fb_pair(
- &dataterm,
- &NonnegRadonRegTerm(α),
- &RadonSquared,
- &auxterm,
- p,
- |_| Err(NotImplemented.into()),
- )
- .map(|Pair(μ, _)| (μ, ())),
- },
- )
- }
- }
- }
- (_, _, _) => {
- bail!("Errors in problem {} setup: auxiliary term or auxiliary variable initialisation given without the other", self.name);
- }
+ let rng = StdRng::seed_from_u64(problem.seed.unwrap_or(42));
+ if let Some(ref d) = problem.domain {
+ let mut domain: UniformRandom<Cube<2, f64>, _> =
+ UniformRandom(d.extract::<[[f64; 2]; 2]>(py)?.into(), rng);
+ self.do_runall_phase2(
+ cli,
+ cli_algorithm_overrides,
+ prefix,
+ problem,
+ &mut domain,
+ py,
+ )
+ } else {
+ self.do_runall_phase2(
+ cli,
+ cli_algorithm_overrides,
+ prefix,
+ problem,
+ &mut FullSpace,
+ py,
+ )
}
- })?;
-
- // Should run the experiment here, going through all algorithms.
- Ok(())
- }
-
- fn algorithm_overrides(&self, alg: DefaultAlgorithm) -> AlgorithmOverrides<f64> {
- self.algorithm_overrides
- .get(&alg)
- .cloned()
- .unwrap_or(Default::default())
+ })
}
}