Mercurial > repos > pointsource_pde / changeset
changeset
Dolfin update, fixes, additional experiment, build instructions.
Wed, 22 Apr 2026 22:32:00 -0500
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Wed, 22 Apr 2026 22:32:00 -0500
- changeset 3
- c3a4f4bb87f7
- parent 1
- a4137aedcb3a
- child 4
- 49b062acace9
Dolfin update, fixes, additional experiment, build instructions.
--- a/.hgignore Thu Feb 26 09:32:12 2026 -0500 +++ b/.hgignore Wed Apr 22 22:32:00 2026 -0500 @@ -10,3 +10,9 @@ *.json *.orig *.png +.venv/ +basix-*/ +dolfinx-*/ +res +res_*/ +results/
--- a/Cargo.toml Thu Feb 26 09:32:12 2026 -0500 +++ b/Cargo.toml Wed Apr 22 22:32:00 2026 -0500 @@ -12,13 +12,13 @@ categories = ["mathematics", "science", "computer-vision"] [dependencies.alg_tools] -version = "~0.4.0-dev" +version = "~0.4.1-dev" path = "../alg_tools" default-features = false features = ["nightly", "pyo3"] [dependencies.pointsource_algs] -version = "~3.0.0-dev" +version = "~3.0.1-dev" path = "../pointsource_algs" [dependencies.measures]
--- a/README.md Thu Feb 26 09:32:12 2026 -0500 +++ b/README.md Wed Apr 22 22:32:00 2026 -0500 @@ -1,27 +1,132 @@ -# Pointsource-PDE +# Codes for “Leak localisation with a measure source convection--diffusion model” + +These are the codes for the numerical demonstrations of the manuscript +_“Leak localisation with a measure source convection--diffusion model”_ +([arXiv:????.????](https://arxiv.org/abs/????.????)) by Thi Tam Dang and Tuomo Valkonen. +It should be relatively easy to to use this package and the algorithms it provides +for other PDE-based point source localisation problems. + +## Building + +Sorry, although the core of program is written in Rust with a modern dependency management and build process, we also have legacy C++ and Python dependencies, i.e., the Fenicsx PDE library. Therefore, the build process is difficult. + +(We admit it, we made a mistake by going with the crowd and using Fenicsx. In the end it would have been less effort to write low-level PDE code in Rust.) -## Prerequisites +### Phase 1: Python and Fenics + +#### Option 1.A that avoids Conda hell, but is more work (macOS and Linux) + +##### Phase 1.A.1: C++ dependencies of Fenics + +First install C++ dependencies using [Homebrew](https://brew.sh). +Even on Linux, **use Homebrew**, or install from official sources; distribution packages are usually obsolete and buggy, often non-standard, and will cause problems (see above). + + brew install openmpi boost pugixml fmt spdlog hdf5-mpi cmake kahip slepc petsc gsl + +(Fenics recommends ParMETIS instead of Kahip, but only the latter is available from Homebrew at the time of writing this.) + +There's no guarantee that this will install compatible versions of these packages. Homebrew, while better than most Linux distributions, is also obsolete in its philosophy: it does not allow easily installing specific versions of packages. Versions known to work are: -## Installation and usage +package | version +---------|-------- +boost | 1.90.0 +cmake | 4.2.1 +fmt | 12.1.0 +hdf5-mpi | 1.14.6 +kahip | 3.22 +petsc | 3.24.3 +pugixml | 1.15 +slepc | 3.24.2 +spdlog | 1.17.0 +gsl | 2.8 + +You can get the list of installed versions with: -### Installing dependencies + brew list --versions openmpi boost pugixml fmt spdlog hdf5-mpi cmake kahip slepc petsc gsl + +##### Phase 1.A.2: Python dependencies of Fenics + +In this source directory, create and activate a virtual environment for Python, and +install Python packages: -Python dependencies are managed Conda and by the Cargo build system of [Rust]. -On some platforms, you can use alternative methods. + python3 -m venv .venv + source .venv/bin/activate + PETSC_DIR=/opt/homebrew/ SLEPC_DIR=/opt/homebrew/ pip install -r requirements.lock + +To not have to activate the virtual environment manually every time, and to not mess up your global settings, it is recommended to install [direnv](https://direnv.net) and put the following in `.envrc` in this directory: + + source .venv/bin/activate + export PYTHONPATH=$(echo .venv/lib/python*/site-packages) + export PYO3_PYTHON="$(which python)" -#### Phase 1: Python +(The last two lines are required later.) +This template is also available in `misc/_envrc`. +For changes `.envrc` to take effect, you should use + + direnv allow + +##### Phase 1.A.3: Fenicx-basix + +Install basix from <https://github.com/FEniCS/basix/releases/tag/v0.10.0.post0> according to +instructions. First do the C++ bit: -##### Most platforms: + tar xzf basix-0.10.0.post0.tar.gz + cd basix-0.10.0.post0/cpp + mkdir build + cd build + cmake .. + make + make install + +Then the Python bit. This has to be done with the `venv` created above, active. + + cd ../../python + pip install . + +##### Phase 1.A.4: Fenicx-dolfinx + +Install dolfinx from <https://github.com/FEniCS/dolfinx/releases/tag/v0.10.0.post5> according to +instructions. First do the C++ bit: -- Fenicsx installed in Conda according to instructions -- SciFEM and SciPY (`conda install conda-forge::scifem conda-forge::scipy`) + tar xzf dolfinx-0.10.0.post5.tar.gz + cd dolfinx-0.10.0.post5/cpp + mkdir build + cd build + cmake .. + make + make install + +Skip the `source /usr/local/lib/dolfinx/dolfinx.conf` recommended at the end +of the compilation. It will likely break things. + +Then the Python bit. This has to be done with the virtual environment created above, active. + + cd ../../python + python -m scikit_build_core.build requires | python -c "import sys, json; print(' '.join(json.load(sys.stdin)))" | xargs pip install + pip install --check-build-dependencies --no-build-isolation . + +If you didn't already do these steps with `direnv` above, you should: -##### Debian/Ubuntu + export PYTHONPATH=$(echo .venv/lib/python*/site-packages) + export PYO3_PYTHON="$(which python)" + + +#### Option 1.B: Conda + +You can *try to* install Fenicsx in Conda according to instructions on the Fenics website. Additionally you need to install `scipy`: -You may br able to use the system package manager instead of Conda, but beware of obsolete versions. + conda create -n fenicsx-env + conda activate fenicsx-env + conda install -c conda-forge fenics-dolfinx=0.10.0 scipy=1.17.1 mpich + +This is, however, unlikely to not work, as Conda, despite its sandboxing separation attempts, conflicts with system packages, or Conda packages have weird ideas. You're likely to run into runtime problems with the FFCX form compiler (*bad* *bad* *bad* idea, running a C compiler runtime) failing due to something, somewhere, in the extremely fragile Conda setup, trying to load system libraries wrongly, etc. -#### Phase 2: Rust +#### Option 1.C: Debian/Ubuntu + +You may be able to use the system package manager, but beware of obsolete and modified versions. As of 2026-03-23, the packages available in Debian/Ubuntu cause massive memory leaks and eventual system crash. + +### Phase 2: Rust You will only need to install the “nightly” Rust compiler and the [GNU Scientific Library] manually. At the time of writing this README, @@ -30,24 +135,16 @@ 1. Install the [Rust] infrastructure (including Cargo) with [rustup]. 2. Install a “nightly” release of the Rust compiler. With rustup, installed in the previous step, this can be done with - ```console - rustup toolchain install nightly - ``` -3. Install [GNU Scientific Library]. On a Mac with [Homebrew] installed, - this can be done with - ```console - brew install gsl - ``` - For other operating systems, suggestions are available in the [rust-GSL] - crate documentation. If not correctly installed, you may need to pass - extra `RUSTFLAGS` options to Cargo in the following steps to locate the - library. + + rustup toolchain install nightly -4. Download [alg_tools] and unpack it under the same directory as this - package. +3. Download [alg_tools], [pointsource_algs], and [measures] and unpack them under + the same directory as this package. [rustup]: https://rustup.rs [alg_tools]: https://tuomov.iki.fi/software/alg_tools/ + [pointsource_algs]: https://tuomov.iki.fi/software/pointsource_algs/ + [measures]: https://tuomov.iki.fi/repos/measures/ [Rust]: https://www.rust-lang.org/ [GNU Scientific Library]: https://www.gnu.org/software/gsl/ [rust-GSL]: https://docs.rs/GSL/6.0.0/rgsl/ @@ -56,41 +153,61 @@ [arXiv:2502.12417]: https://arxiv.org/abs/2502.12417 [doi:10.46298/jnsao-2023-10433]: http://doi.org/10.46298/jnsao-2023-10433 -### Building and running the experiments +### Linux / further patching + +Due to both Fenics and typical Linux system being completely broken, you may need to do further patching to get things to compile: -To compile and install the program, use -```console -cargo install --path=. -``` + 1. I had to set (in my [direnv](https://direnv.net) `.envrc`) + + export PKG_CONFIG_PATH=/home/linuxbrew/.linuxbrew/lib/pkgconfig:/usr/local/lib/pkgconfig/:/usr/lib/aarch64-linux-gnu/pkgconfig/ + export LD_LIBRARY_PATH=/home/linuxbrew/.linuxbrew/lib + + 2. Some libraries, in particular `libfmt` and `libspdlog` installed in Homebrew, may conflict with system versions, that must be removed. + Lack of proper sandboxing in legacy Linux distributions, effectively prohibits multiple versions of the same library. + 2. I had to add `Libs` in `/usr/local/lib/pkgconfig/dolfinx.pc` the bit `-L/home/linuxbrew/.linuxbrew/lib/ -lopenblas`. Nothing in the fenics stack seems to explicitly require it. Basix, that depends on openblas, is entirely missing a `pkg-config` file. + 3. Also `export export OMP_NUM_THREADS=1` (in `.envrc`). We don't do MPI. We cannot do MPI in Fenics' lame “it's all just parallel solution of PDEs, with no other computation, ever” aka “single-program multiple-data, with no controller at all” way. If you don't do this, you may have multiple threads wasting CPU just being there. We try to control the thread count in our code, but OpenMPI on Linux doesn't seem to respect it. + +## Building and running the experiments + +To compile the program, run + + cargo build --release + When doing this for the first time, several dependencies will be downloaded. -Now you can run the default set of experiments with -``` -pointsource_pde -o results -``` -The `-o results` option tells `pointsource_pde` to write results in the -`results` directory. The option is required. +Now you can run the experiments in the article with + + cargo run --release -- \ + -o results -a radon_sliding_fb -a radon_fb --max-iter 20000 + experiments/laser_and_mirrors_aux.py experiments/laser_and_mirrors_aux2.py + +The `-o results` option tells `pointsource_pde` to write results in the `results` directory. +The other options indicate the algorithms and experiments to run, as well as the maximum number +of iterations. +The double-dash separates the options for the Cargo build system and `pointsource_pde`. -Alternatively, you may build and run the program without installing with -```console -cargo run --release -- -o results -``` -The double-dash separates the options for the Cargo build system -and `pointsource_pde`. +### Visualising the results + +The results may be plotted with + + python3 ./plot.py results/laser_and_mirrors_aux/radon_sliding_fb -### Documentation +Vary the path to `laser_and_mirrors_aux2` and `radon_fb` for the alternative experiment and basic algorithm. + +The script `misc/copy_results.sh` may be generate the images and copy the results in the manuscript +to `../gasleak`. + +## Documentation Use the `--help` option to get an extensive listing of command line options to customise algorithm parameters and the experiments performed. -## Internals +### Internals If you are interested in the program internals, the integrated source code documentation may be built and opened with -```console -cargo doc # build dependency docs -misc/cargo-d --open # build and open KaTeX-aware docs for this crate -``` -The `cargo-d` script ensures that KaTeX mathematics is rendered in the -generated documentation through an ugly workaround. Unfortunately, -`rustdoc`, akin to Rust largely itself, is stuck in 80's 7-bit gringo ASCII -world, and does not support modern markdown features, such as mathematics. + + cargo doc # build dependency docs + misc/cargo-d --open # build and open KaTeX-aware docs for this crate + +The `cargo-d` script ensures that KaTeX mathematics is rendered in the generated documentation through an ugly workaround. Unfortunately, `rustdoc` is stuck in 80's 7-bit gringo ASCII world, +and does not support modern markdown features, such as mathematics.
--- a/build.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/build.rs Wed Apr 22 22:32:00 2026 -0500 @@ -2,6 +2,8 @@ Build script for `pointsource_pde` */ +#![feature(path_is_empty)] + //use pyo3::{prepare_freethreaded_python, types::PyModule, PyErr, Python}; use conda_build::python_config; use std::env::split_paths; @@ -38,7 +40,7 @@ // Ok::<_, PyErr>(()) // })?; - let (pyc_e, py_is_conda) = python_config(); + let (pyc_e, _prefix_override, _py_type) = python_config(); let py_path = pyc_e.and_then(|pyc| pyc.exec_prefix_path())?; let py_exec = py_path.join("bin").join("python3"); let import_check_code = REQUIRED_PYTHON_MODULES @@ -49,9 +51,8 @@ .status()?; if !res.success() { return Err(std::io::Error::other(format!( - "Failed to load required Python modules (python execuitable {}{})", + "Failed to load required Python modules (python execuitable {})", py_exec.display(), - if py_is_conda { " from Conda" } else { "" } ))); } @@ -130,7 +131,9 @@ fn add_dep_ldflags(dep: &str) { if let Some(ldflags) = std::env::var_os(&dep) { for lp in split_paths(&ldflags) { - println!("cargo:rustc-link-arg=-Wl,{}", lp.display()); + if !lp.is_empty() { + println!("cargo:rustc-link-arg=-Wl,{}", lp.display()); + } } } else { panic!("{dep} unset.");
--- a/conda-build/src/lib.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/conda-build/src/lib.rs Wed Apr 22 22:32:00 2026 -0500 @@ -36,26 +36,32 @@ /// Looks up conda Python pub fn python_conda() -> Option<OsString> { option_env!("CONDA_PREFIX").and_then(|conda_prefix| { - /* This does not point to the right place. - option_env!("CONDA_PYTHON_EXE") - .map(OsString::from) - .or_else(|| { */ let py = PathBuf::from(conda_prefix).join("bin/python"); if py.exists() { Some(py.into()) } else { None } - /* })*/ }) } /// Returns a the python configuration and indication whether it's from Conda -pub fn python_config() -> (Result<PythonConfig, PyCError>, bool) { - if let Some(py) = python_conda() { - (PythonConfig::interpreter(py), true) +pub fn python_config() -> ( + Result<PythonConfig, PyCError>, + Option<PathBuf>, + Option<&'static str>, +) { + if let Some(py) = option_env!("VIRTUAL_ENV") { + let interp = PathBuf::from(py).join("bin/python"); + ( + PythonConfig::interpreter(interp), + Some(py.into()), + Some("venv"), + ) + } else if let Some(py) = python_conda() { + (PythonConfig::interpreter(py), None, Some("conda")) } else { - (Ok(PythonConfig::new()), false) + (Ok(PythonConfig::new()), None, None) } }
--- a/experiments/laser_and_mirrors.py Thu Feb 26 09:32:12 2026 -0500 +++ b/experiments/laser_and_mirrors.py Wed Apr 22 22:32:00 2026 -0500 @@ -18,20 +18,18 @@ # Setup routine def setup(prefix): - dat, auxtrue, μ_bound, μ_true, plot_factory = laser_and_mirrors_aux.generic_setup( + dat, auxtrue, _μ_bound, μ_true, plot_factory = laser_and_mirrors_aux.generic_setup( prefix ) μ0 = DiscreteMeasure_2_f64([]) - reg = RegTerm.NonnegRadon( - 10e-7 # 5e-7 # 3e-7 - ) + reg = RegTerm.NonnegRadon(1.5e-6) - dat_simple = ComposeFnWithOperator(dat, InjectSecond(auxtrue), xbound=μ_bound) + dat_simple = ComposeFnWithOperator(dat, InjectSecond(auxtrue)) dat_simple.curvature_bound_components = lambda: (None, None) - print("Initial value:", dat_simple.apply(μ0)) - print("Value at true μ:", dat_simple.apply(μ_true)) + print("Initial data term value:", dat_simple.apply(μ0)) + print("Data term value at true μ:", dat_simple.apply(μ_true)) return Problem( dataterm=dat_simple, regterm=reg, μinit=μ0, plot_factory=plot_factory
--- a/experiments/laser_and_mirrors_aux.py Thu Feb 26 09:32:12 2026 -0500 +++ b/experiments/laser_and_mirrors_aux.py Wed Apr 22 22:32:00 2026 -0500 @@ -13,9 +13,9 @@ SumOfSeparableFunctions, ) from pointsource_pde.convection_diffusion import ( + BoxedQuadraticRegularisation, ConvectionDiffusion, PlotFactory, - QuadraticRegularisation, XBound, ) from pointsource_pde.laser_sampling import LaserSampling @@ -27,12 +27,14 @@ # Fine-tune algorithms alg_finetune = { "tau0": 0.99, - "theta0": 0.01, + "theta0": 100, "inner_method": "PP", "inner_pdps_τσ0": (19.8, 0.05), "inner_pp_τ": [1000.0, 1000.0], "merge": True, "fitness_merging": True, + "merge_radius": 0.1, + "merge_interp": False, } algorithm_overrides = { @@ -44,12 +46,17 @@ # Setup shared between experiment with known/unknown conductivity and diffusivity -def generic_setup(prefix, simple=True): - - # Create true source - μ_true = DiscreteMeasure_2_f64( +def generic_setup( + prefix, + simple=True, + rng=None, + μ_true=DiscreteMeasure_2_f64( [([0.1, 0.3], 0.08), ([0.4, 0.25], 0.05), ([0.25, 0.13], 0.06)] - ) + ), + k=0.01, + r=0.5, + θ=30 * np.pi / 180, +): # Create convection-diffusion PDE pde = ConvectionDiffusion( @@ -57,6 +64,8 @@ g=lambda x: 0.0 * x[1], domain_size=[0, 0.5, 0, 0.5], num_steps=50, + # override_lipschitz=4.0, + # override_lipschitz_pair=(4.0, 40.0), ) # Set up lasers and mirrors @@ -91,12 +100,10 @@ if simple: # Scalar convectivity and diffusivity - r = 0.5 - θ = 30 * np.pi / 180 c1 = r * np.cos(θ) c2 = r * np.sin(θ) - k = 0.01 # diffusive # k = 0.001 # high wind + print("True aux variables: k = %f, c1 = %f, c2 = %f" % (k, c1, c2)) else: # Convectivity and diffusivity fields def velocity_x(x): @@ -120,12 +127,14 @@ # Simulate measurements u_n_list = pde.apply((μ_true, (k, (c1, c2)))) b_true = [L_i.apply(u) for u in u_n_list] - b = [b_i + L_i.noise(0.01 * norm(b_i) / np.sqrt(np.size(b_i))) for b_i in b_true] + b = [ + b_i + L_i.noise(0.01 * norm(b_i) / np.sqrt(np.size(b_i)), rng) for b_i in b_true + ] # Estimate bounds on domains and values - dat_bound = 10 * max([norm(b_i) for b_i in b]) + dat_bound = 2 * max([norm(b_i) for b_i in b]) μ_bound = 1.3 * sum([alpha for _point, alpha in μ_true.iter_padded()]) - xbound = XBound(mu_dual=μ_bound, k_min=k, k_max=k, c_Linf=r) + pde.xbound = XBound(mu_dual=μ_bound, k_min=k * 0.9, k_max=k * 1.1, c_Linf=r * 0.9) # Create data term dat = ComposeFnWithOperator( @@ -139,14 +148,12 @@ λ=pde.dt / len(beams), ) for data in b - ] + ], ), pde, - xbound=[xbound], - xbound_pair=[xbound], ) - # Plot true and initial data + # Save true and initial data plot_factory = PlotFactory(pde) plotter = plot_factory.produce(prefix) ω, _ = dat.diff((DiscreteMeasure_2_f64([]), auxtrue)) @@ -156,41 +163,53 @@ true_mu=[[point[0], point[1], alpha] for point, alpha in μ_true.iter_padded()], omega0=ω.x.array, true_omega=ω_true.x.array, + k=k, + c1=c1, + c2=c2, true_u_n_list_array=[u_n_list[i].x.array for i in plotter.plot_frames], ) - return dat, auxtrue, μ_bound, μ_true, plot_factory + return dat, auxtrue, μ_bound, μ_true, plot_factory, pde + + +def relnoise(v, σ, rng, level=None): + if level is None: + level = abs(v) + return float(rng.normal(v, σ * level)) +# Setup the problem def setup(prefix): - dat, auxtrue, _μ_bound, μ_true, plot_factory = generic_setup(prefix) + rng = np.random.default_rng(seed=31337) + + dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup(prefix, rng=rng) - reg = RegTerm.NonnegRadon(10e-7) + # Override Lipschitz parameter + pde.override_lipschitz = (4.0,) + pde.override_lipschitz_pair = (4.0, 4.0) + # print("diff_chain_lipschitz_factor (modified)", pde.diff_chain_lipschitz_factor()) + # l1, l2 = pde.diff_chain_lipschitz_factor_pair() + # print("diff_chain_lipschitz_factor_pair (modified) ", l1, ", ", l2) + + reg = RegTerm.NonnegRadon(1.5e-6) μ0 = DiscreteMeasure_2_f64([]) - aux0 = auxtrue - aux = QuadraticRegularisation(0.001, aux0) + (k, (c1, c2)) = auxtrue + aux0 = ( + max(0.001, relnoise(k, 0.02, rng)), + ( + relnoise(c1, 0.2, rng), + relnoise(c2, 0.2, rng), + ), + ) + print("aux init ", aux0) + aux = BoxedQuadraticRegularisation((0.001, -1.0), (1.0, 1.0), (3.0, 0.0005), aux0) ax = aux.apply(aux0) inival = dat.apply((μ0, aux0)) - print("Initial value:", inival + ax, ax) - print("Value at true μ:", dat.apply((μ_true, auxtrue)) + ax) - - # auxinit = np.c_[ - # np.zeros_like(k), - # np.zeros_like(c1), - # np.zeros_like(c2), - # ] - # print(np.shape(auxinit)) + print("Initial data term value:", inival + ax) + print("Data term value at true μ:", dat.apply((μ_true, auxtrue)) + ax) - # We need Θ² such that ⟨F'(μ+Δ, aux-F'(μ, aux)|Δ⟩ ≤ Θ²|γ|(c_2), where Δ=(π_♯^1-π_♯^0)γ. - # We have - # - # ⟨F'(μ+Δ, aux)-F'(μ, aux)|Δ⟩ ≤ ‖⟨F'(μ+Δ, aux)-F'(μ, aux)‖ ‖Δ‖ - # ≤ L_{F'(·, aux)} ‖Δ‖², - # - # so Θ² ≥ L_{F'(·, aux)} ‖Δ‖² / |γ|(c_2). Thsi doesn't give anything useful if the - # transport is very small in distance. - - dat.curvature_bound_components = lambda: (None, None) # 1.0 + # No curvature bound given: θ is aboslute. + dat.curvature_bound_components = lambda: (None, None) return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/experiments/laser_and_mirrors_aux2.py Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,60 @@ +# +# Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity. +# Alternative parametrisation +# +import laser_and_mirrors_aux +import numpy as np +from laser_and_mirrors_aux import generic_setup, relnoise +from measures import DiscreteMeasure_2_f64 +from pointsource_pde import Problem, RegTerm +from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation + +# Give name to the problem +name = "laser_and_mirrors_aux2" + +# Override algorithm settings +algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides + + +# Setup the problem +def setup(prefix): + rng = np.random.default_rng(seed=31337) + + dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup( + prefix, + rng=rng, + μ_true=DiscreteMeasure_2_f64([([0.2, 0.3], 0.15), ([0.4, 0.1], 0.04)]), + k=0.02, + r=0.1, + θ=120 * np.pi / 180, + ) + + # Override Lipschitz parameter + pde.override_lipschitz = (4.0,) + pde.override_lipschitz_pair = (4.0, 4.0) + # print("diff_chain_lipschitz_factor (modified)", pde.diff_chain_lipschitz_factor()) + # l1, l2 = pde.diff_chain_lipschitz_factor_pair() + # print("diff_chain_lipschitz_factor_pair (modified) ", l1, ", ", l2) + + reg = RegTerm.NonnegRadon(1.5e-6) + + μ0 = DiscreteMeasure_2_f64([]) + (k, (c1, c2)) = auxtrue + aux0 = ( + max(0.001, relnoise(k, 0.02, rng)), + ( + relnoise(c1, 0.2, rng), + relnoise(c2, 0.2, rng), + ), + ) + print("aux init ", aux0) + aux = BoxedQuadraticRegularisation((0.001, -1.0), (1.0, 1.0), (3.0, 0.0005), aux0) + ax = aux.apply(aux0) + inival = dat.apply((μ0, aux0)) + print("Initial data term value:", inival + ax) + print("Data term value at true μ:", dat.apply((μ_true, auxtrue)) + ax) + + # No curvature bound given: θ is aboslute. + dat.curvature_bound_components = lambda: (None, None) + + return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)
--- a/include/dolfinx_access/function.h Thu Feb 26 09:32:12 2026 -0500 +++ b/include/dolfinx_access/function.h Wed Apr 22 22:32:00 2026 -0500 @@ -1,5 +1,7 @@ #pragma once #include "dolfinx/geometry/BoundingBoxTree.h" +#pragma once +#include "dolfinx/mesh/Mesh.h" #include <cmath> #pragma once @@ -11,10 +13,11 @@ namespace dolfinx_access { struct FunctionInfo; - typedef geometry::BoundingBoxTree<double> BBTree_f64; + typedef dolfinx::geometry::BoundingBoxTree<double> BBTree_f64; + extern void drop_Function_f64(Function_f64* f); - extern int32_t cell_Mesh_f64(std::shared_ptr<const mesh::Mesh<double>> mesh, + extern int32_t cell_Mesh_f64(std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh, const std::array<double, 3>& x); extern int32_t cell_FunctionSpace_f64(FunctionSpace<double> const* v, @@ -46,4 +49,5 @@ Function_f64* clone_Function_f64(Function_f64 const* f); + void scatter_fwd_Function_f64(Function_f64* f); } // namespace dolfinx_access
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/_envrc Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,3 @@ +source .venv/bin/activate +export PYTHONPATH=$(echo .venv/lib/python*/site-packages) +export PYO3_PYTHON="$(which python)"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/cargo-d Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,4 @@ +#!/bin/sh +RUSTDOCFLAGS="--html-in-header misc/katex-header.html" cargo d --no-deps "$@" + + \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/copy_results.sh Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,27 @@ +#!/bin/sh +dest=../gasleak +iter=20000 +files="radon_sliding_fb_log.txt radon_fb_log.txt radon_sliding_fb_config.json radon_fb_config.json radon_sliding_fb_stats.json radon_fb_stats.json" + +src1=results/laser_and_mirrors_aux +src2=results/laser_and_mirrors_aux2 + +./plot.py ${src1}/radon_fb --simple --iter ${iter} --save-only +./plot.py ${src1}/radon_sliding_fb --simple --iter ${iter} --save-only +cp ${src1}/radon_sliding_fb/res_${iter}.pdf ${dest}/aux_res_sliding_${iter}.pdf +cp ${src1}/radon_sliding_fb/true_${iter}.pdf ${dest}/aux_true_${iter}.pdf +cp ${src1}/radon_fb/res_${iter}.pdf ${dest}/aux_res_${iter}.pdf +for file in $files; do + cp ${src1}/$file ${dest}/aux_$file +done + +./plot.py ${src2}/radon_fb --simple --iter ${iter} --save-only +./plot.py ${src2}/radon_sliding_fb --simple --iter ${iter} --save-only +cp ${src2}/radon_sliding_fb/res_${iter}.pdf ${dest}/aux2_res_sliding_${iter}.pdf +cp ${src2}/radon_sliding_fb/true_${iter}.pdf ${dest}/aux2_true_${iter}.pdf +cp ${src2}/radon_fb/res_${iter}.pdf ${dest}/aux2_res_${iter}.pdf +cp ${src2}/radon_sliding_fb_log.txt ${dest}/aux2_radon_sliding_fb_log.txt +cp ${src2}/radon_fb_log.txt ${dest}/aux2_radon_fb_log.txt +for file in $files; do + cp ${src2}/$file ${dest}/aux2_$file +done
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/misc/katex-header.html Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,63 @@ +<link + rel="stylesheet" + href="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.css" + integrity="sha384-Xi8rHCmBmhbuyyhbI88391ZKP2dmfnOl4rT9ZfRI7mLTdk1wblIUnrIq35nqwEvC" + crossorigin="anonymous" +/> +<script + defer + src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/katex.min.js" + integrity="sha384-X/XCfMm41VSsqRNQgDerQczD69XqmjOOOwYQvr/uuC+j4OPoNhVgjdGFwhvN02Ja" + crossorigin="anonymous" +></script> +<script + defer + src="https://cdn.jsdelivr.net/npm/katex@0.16.0/dist/contrib/auto-render.min.js" + integrity="sha384-+XBljXPPiv+OzfbB3cVmLHf4hdUFHlWNZN5spNQ7rmHTXpd7WvJum6fIACpNNfIR" + crossorigin="anonymous" + onload="renderMathInElement(document.body);" +></script> +<script> + document.addEventListener("DOMContentLoaded", function () { + renderMathInElement(document.body, { + delimiters: [ + { left: "$$", right: "$$", display: true }, + { left: "\\(", right: "\\)", display: false }, + { left: "$", right: "$", display: false }, + { left: "\\[", right: "\\]", display: true }, + ], + macros: { + "\\iprod": "{\\langle #1, #2\\rangle}", + "\\dualprod": "{\\langle #1| #2\\rangle}", + "\\norm": "{\\|#1\\|}", + "\\abs": "{|{#1}|}", + "\\grad": "\\nabla", + "\\isect": "\\cap", + "\\union": "\\cup", + "\\Isect": "\\bigcap", + "\\Union": "\\bigcup", + "\\linear": "\\mathbb{L}", + "\\downto": "\\searrow", + "\\upto": "\\nearrow", + "\\setto": "\\rightrightarrows", + "\\Meas": "\\mathcal{M}", + "\\B": "\\mathbb{B}", + "\\subdiff": "\\partial", + "\\inv": "{#1}^{-1}", + "\\pinv": "{#1}^\\dagger", + "\\pinvstar": "{#1}^{\\dagger*}", + "\\freevar": "\\,\\boldsymbol\\cdot\\,", + "\\MYMATHOP": "\\mathop{\\mathrm{#1}}", + "\\prox": "\\MYMATHOP{prox}", + "\\proj": "\\MYMATHOP{proj}", + "\\supp": "\\MYMATHOP{supp}", + "\\soft": "\\MYMATHOP{soft}", + "\\sign": "\\MYMATHOP{sign}", + "\\ran": "\\MYMATHOP{ran}", + "\\ker": "\\MYMATHOP{ker}", + "\\Id": "\\MYMATHOP{Id}", + "\\defeq": ":=", + }, + }); + }); +</script>
--- a/nanobind-sys/build.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/nanobind-sys/build.rs Wed Apr 22 22:32:00 2026 -0500 @@ -17,7 +17,7 @@ fn main() -> Result<(), anyhow::Error> { // Need to build it - let (pyc_e, py_conda) = python_config(); + let (pyc_e, prefix_override, py_type) = python_config(); // This is very clumsy due to PythonConfig::Error not supporting // conversion into std::error::Error, just std::io::Error. @@ -28,11 +28,11 @@ }) .map_err(|e| anyhow::Error::from(std::io::Error::from(e)))?; - let nanobind_root = pyc - .prefix_path().ok() - .zip(pyc.semantic_version().ok()) - .and_then(|(prefix, version)| { - info!("{}Python {} found at prefix {}", if py_conda { "Conda " } else { ""}, version, prefix.display()); + let prefix = prefix_override.ok_or_else(|| pyc.prefix_path()).unwrap(); + let version = pyc.semantic_version().unwrap(); + + let nanobind_root = ({ + info!("{}{}Python {} found at prefix {}", py_type.unwrap_or(""), if py_type.is_none() { "" } else { " "}, version, prefix.display()); let nanobind_root = prefix .join("lib") .join(format!("python{}.{}", version.major, version.minor))
--- a/plot.py Thu Feb 26 09:32:12 2026 -0500 +++ b/plot.py Wed Apr 22 22:32:00 2026 -0500 @@ -1,3 +1,5 @@ +#!.venv/bin/python3 +import argparse import os import re import sys @@ -12,16 +14,15 @@ import matplotlib.tri as tri import numpy as np from dolfinx import fem - from src.convection_diffusion import ConvectionDiffusion -def find_files(directory): +def find_files(dirname): """ Return a list of Path objects for files named 'fubar%d.npz' in the given directory, sorted by the integer %d. """ - directory = Path(directory) + directory = Path(dirname) pattern = re.compile(r"^res_(\d+)\.npz$") matches = [] @@ -39,7 +40,7 @@ return matches -def plot(prefix): +def plot(prefix, simple=False, shaded=False, init_iter=None, save_only=False): quantisation = 32 iter_files = find_files(prefix) @@ -59,16 +60,27 @@ triang = tri.Triangulation(x_coords, y_coords) m = 5 - fig, (axus, axws, axjs, ax_wpbar_plus) = plt.subplots(4, m, figsize=(15, 4)) - plt.tight_layout() + if simple: + fig, axus = plt.subplots(1, m) + all_axes = lambda: axus + fig.set_layout_engine("compressed") + else: + fig, (axus, axws, axjs, ax_wpbar_plus) = plt.subplots(4, m, figsize=(15, 4)) + all_axes = lambda: chain(axus, axws, axjs, ax_wpbar_plus) + plt.tight_layout(pad=0.01) + fig.subplots_adjust( + left=0.02, right=0.95, bottom=0.02, top=0.93, wspace=0.2, hspace=0.2 + ) + fig.set_size_inches(13, 8) - fig.subplots_adjust( - left=0.02, right=0.95, bottom=0.02, top=0.93, wspace=0.2, hspace=0.2 - ) - for ax in chain(axus, axws, axjs, ax_wpbar_plus): - ax.set_xticklabels([]) - ax.set_yticklabels([]) - ax.set_aspect("equal") + + if not simple: + info = plt.text( + 1.0, + -0.07, + "", + horizontalalignment="right", + ) mpl.rcParams["axes.labelsize"] = 9 @@ -76,10 +88,12 @@ sm = cm.ScalarMappable(norm=norm, cmap="viridis") sm.set_array([]) # REQUIRED (matplotlib quirk) colorbar_u = fig.colorbar(sm, ax=axus[-1]) - colorbar_w = fig.colorbar(sm, ax=axws[-1]) - colorbar_j = fig.colorbar(sm, ax=axjs[-1]) + if not simple: + colorbar_w = fig.colorbar(sm, ax=axws[-1]) + colorbar_j = fig.colorbar(sm, ax=axjs[-1]) def do_plot(index_and_file, show_true=False): + (iter, file) = index_and_file print("🎨 Plotting iteration %d..." % iter) @@ -92,11 +106,6 @@ ω0_min = ω0.real.min() ω0_max = ω0.real.max() - for ax in chain(axus, axws, axjs, ax_wpbar_plus): - # ax.clear() - for artist in chain(ax.lines, ax.collections, ax.images): - artist.remove() - data = np.load(file) u_n_list_array = data["u_n_list_array"] w_n_list_array = data["w_n_list_array"] @@ -122,11 +131,7 @@ j_max = max(j.real.max() for j in j_n_list_array) μ = list(filter(lambda x: x[2] != 0.0 and not np.isnan(x).any(), mu)) - μ_x = list(map(lambda x: x[0], μ)) - μ_y = list(map(lambda x: x[1], μ)) μ_alpha = list(map(lambda x: x[2], μ)) - true_μ_x = list(map(lambda x: x[0], true_μ)) - true_μ_y = list(map(lambda x: x[1], true_μ)) true_μ_alpha = list(map(lambda x: x[2], true_μ)) if len(μ_alpha) == 0: @@ -143,78 +148,248 @@ ms = lambda m: int(1 + (m - alpha_base) / alpha_scale * 10) def plot_array( - name, t_idx, ax, u_array, u_min, u_max, colorbar=None, measure=False + name, + t_idx, + ax, + u_array, + u_min, + u_max, + μ, + true_μ, + colorbar=None, + measure=False, ): if u_min == u_max: levels = [u_min, u_min + 1e-9] else: levels = np.linspace(u_min, u_max, quantisation) try: - contour = ax.tricontourf(triang, u_array, levels=levels, cmap="viridis") + if shaded: + concentration = ax.tripcolor( + triang, + u_array, + cmap="viridis", + shading="gouraud", + ) + else: + concentration = ax.tricontourf( + triang, + u_array, + levels=levels, + cmap="viridis", + antialiased=False, + zorder=-10, + ) + ax.set_rasterization_zorder(0) + if colorbar: - colorbar.update_normal(contour) + colorbar.update_normal(concentration) + fmt = mpl.ticker.ScalarFormatter(useOffset=False) + fmt.set_scientific(False) + colorbar.ax.yaxis.set_major_formatter(fmt) + colorbar.update_ticks() + colorbar.ax.yaxis.get_offset_text().set_visible(False) + off = colorbar.ax.yaxis.get_offset_text() + off.set_x(10) + except Exception as e: print(e) - if measure: + if μ is not None: + μ_x = list(map(lambda x: x[0], μ)) + μ_y = list(map(lambda x: x[1], μ)) + μ_alpha = list(map(lambda x: x[2], μ)) for x, y, m in zip(μ_x, μ_y, μ_alpha): ax.plot([x], [y], "ro", markersize=ms(m), label="Sources") + if true_μ is not None: + true_μ_x = list(map(lambda x: x[0], true_μ)) + true_μ_y = list(map(lambda x: x[1], true_μ)) + true_μ_alpha = list(map(lambda x: x[2], true_μ)) for x, y, m in zip(true_μ_x, true_μ_y, true_μ_alpha): - ax.plot([x], [y], "kx", markersize=ms(m), label="True sources") + ax.plot([x], [y], "wx", markersize=ms(m), label="True sources") if t_idx >= 0: - ax.set_title(f"%s; t = {t_idx:.1f}" % name) + ax.set_title( + f"$%s_{{{t_idx:.1f}}}$" % name, + fontsize=20 if simple else 14, + ) else: - ax.set_title("%s" % name) + ax.set_title( + "%s" % name, + fontsize=20 if simple else 14, + ) ax.set_aspect("equal") n = len(u_n_list_array) # frames = list(map(lambda i: i * (n - 1) // (m - 1), range(0, m))) frames = range(0, n) - for i, axu, axw, axj, t_idx in zip( - frames, - axus, - axws, - axjs, - times, # map(lambda i: i / (n - 1), frames), - ): - plot_array( - "u", t_idx, axu, u_n_list_array[i].real, u_min, u_max, colorbar_u, True - ) - plot_array( - "w", t_idx, axw, w_n_list_array[i].real, w_min, w_max, colorbar_w, True - ) - plot_array( - "j", t_idx, axj, j_n_list_array[i].real, j_min, j_max, colorbar_j, True + + for ax in all_axes(): + ax.clear() + + if not simple: + for i, axu, axw, axj, t_idx in zip( + frames, + axus, + axws, + axjs, + times, # map(lambda i: i / (n - 1), frames), + ): + plot_array( + "u", + t_idx, + axu, + u_n_list_array[i].real, + u_min, + u_max, + μ, + true_μ, + colorbar_u, + ) + plot_array( + "w", + t_idx, + axw, + w_n_list_array[i].real, + w_min, + w_max, + μ, + true_μ, + colorbar_w, + ) + plot_array( + "j", + t_idx, + axj, + j_n_list_array[i].real, + j_min, + j_max, + μ, + true_μ, + colorbar_j, + ) + + if show_true: + for i, ax, t_idx in zip( + frames, + ax_wpbar_plus, + map(lambda i: i / (n - 1), frames), + ): + plot_array( + "û", + t_idx, + ax, + true_u_n_list_array[i].real, + u_min, + u_max, + μ, + true_μ, + colorbar_u, + ) + else: + plot_array( + "w̄ₚ", + -1, + ax_wpbar_plus[0], + wp_bar_array.real, + w_min, + w_max, + μ, + None, + ) + plot_array( + "ω₀", + -1, + ax_wpbar_plus[1], + ω0.real, + w_min, + w_max, + μ, + None, + ) + # plot_array("ω̂", -1, ax_wpbar_plus[2], true_ω.real, w_min, w_max) + + # Has to be here after everything, for some reason + for ax in all_axes(): + ax.set_xticklabels([]) + ax.set_yticklabels([]) + ax.set_aspect("equal") + + plt.suptitle( + "Convection-Diffusion, iteration %d; len(μ) = %d; μ_min = %f, μ_max = %f, μ_median = %f;" + % ( + iter, + len(μ), + alpha_mi, + alpha_ma, + alpha_me, + ), + fontsize=14, ) - if show_true: - for i, ax, t_idx in zip( + str = "(k, c, c1) = (%g, %g, %g)" % (data["k"], data["c1"], data["c2"]) + if "k" in data0: + str += "; true (%g, %g, %g)" % (data0["k"], data0["c1"], data0["c2"]) + info.set_text(str) + else: # simple + for i, axu, t_idx in zip( frames, - ax_wpbar_plus, - map(lambda i: i / (n - 1), frames), + axus, + times, # map(lambda i: i / (n - 1), frames), ): - plot_array( - "û", - t_idx, - ax, - true_u_n_list_array[i].real, - u_min, - u_max, - colorbar_u, - True, - ) + if show_true: + plot_array( + "û", + t_idx, + axu, + true_u_n_list_array[i].real, + u_min, + u_max, + None, + true_μ, + colorbar_u, + ) + else: + plot_array( + "u", + t_idx, + axu, + u_n_list_array[i].real, + u_min, + u_max, + μ, + true_μ, + colorbar_u, + ) + + def do_save(iter, show_true): + if not show_true: + filename = "%s/res_%d.pdf" % (prefix, iter) else: - plot_array("w̄ₚ", -1, ax_wpbar_plus[0], wp_bar_array.real, w_min, w_max) - plot_array("ω₀", -1, ax_wpbar_plus[1], ω0.real, w_min, w_max) - # plot_array("ω̂", -1, ax_wpbar_plus[2], true_ω.real, w_min, w_max) + filename = "%s/true_%d.pdf" % (prefix, iter) + print("Saving ", filename) + plt.savefig( + filename, + dpi=300, + bbox_inches="tight", + pad_inches=0, + ) - plt.suptitle( - "Convection-Diffusion, iteration %d; len(μ) = %d; μ_min = %f, μ_max = %f, μ_median = %f" - % (iter, len(μ), alpha_mi, alpha_ma, alpha_me), - fontsize=14, + if init_iter is None: + k = 0 + else: + k, _ = next( + filter(lambda fubar: fubar[1][0] == init_iter, enumerate(iter_files)) ) - # plt.savefig("solution_evolution_%d.png" % iter, dpi=300) + + do_plot(iter_files[k], show_true=False) - state = {"k": 0, "show_true": False} + if save_only: + iter = iter_files[k][0] + do_save(iter, False) + do_plot(iter_files[k], show_true=True) + do_save(iter, True) + return + + state = {"k": k, "show_true": False} def on_key(event): k0 = state["k"] @@ -238,13 +413,14 @@ k = k0 elif event.key == "q": sys.exit() + elif event.key == "z": + iter, _file = iter_files[k0] + do_save(iter, state["show_true"]) if k is not None: state["k"] = k do_plot(iter_files[k], state["show_true"]) fig.canvas.draw() - do_plot(iter_files[0]) - fig.canvas.mpl_connect("key_press_event", on_key) plt.show() @@ -265,4 +441,17 @@ # print("Saved: solution_evolution.png + solution_time.png") -plot(sys.argv[1]) +parser = argparse.ArgumentParser( + prog="plot.py", + description="Plots results of pointsource_pde", +) + +parser.add_argument("filename") +parser.add_argument("--simple", action="store_true") +parser.add_argument("--shaded", action="store_true") +parser.add_argument("--iter", type=int, action="store") +parser.add_argument("--save-only", action="store_true") + +args = parser.parse_args() + +plot(args.filename, args.simple, args.shaded, args.iter, args.save_only)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/requirements.lock Wed Apr 22 22:32:00 2026 -0500 @@ -0,0 +1,28 @@ +cffi==2.0.0 +contourpy==1.3.3 +cycler==0.12.1 +Cython==3.2.4 +fenics-ffcx==0.10.1.post0 +fenics-ufl==2025.2.1 +fonttools==4.62.1 +kiwisolver==1.5.0 +matplotlib==3.10.8 +mpi4py==4.1.1 +nanobind==2.12.0 +ninja==1.13.0 +numpy==2.4.3 +packaging==26.0 +pathspec==1.0.4 +petsc4py==3.24.5 +pillow==12.1.1 +pycparser==3.0 +pyparsing==3.3.2 +python-dateutil==2.9.0.post0 +scikit_build_core==0.12.2 +scipy==1.17.1 +setuptools==82.0.1 +six==1.17.0 +slepc4py==3.24.3 +uv==0.11.2 +wheel==0.46.3 +argparse==1.4.0
--- a/src/compose.py Thu Feb 26 09:32:12 2026 -0500 +++ b/src/compose.py Wed Apr 22 22:32:00 2026 -0500 @@ -32,19 +32,17 @@ res = max(res, f_i.diff_lipschitz_factor()) return res - def diff_bound(self, xbound=None): + def diff_bound(self): res = 0 for f_i in self.fnlist: - res = max(res, f_i.diff_bound(xbound=xbound)) + res = max(res, f_i.diff_bound()) return res class ComposeFnWithOperator: - def __init__(self, f, op, xbound=None, xbound_pair=None): + def __init__(self, f, op): self.f = f self.op = op - self.xbound = xbound - self.xbound_pair = xbound_pair def apply(self, *args): return self.f.apply(self.op.apply(*args)) @@ -61,40 +59,67 @@ (a, v) = self.f.apply_and_diff(w) return (a, self.op.diff_adjdir(v, *args, apply_result=w)) - def diff_lipschitz_factor(self): - # ‖∇A(x)^*∇F(A(x)) - ∇A(y)^*∇F(A(y))‖ - # = ‖[∇A(x)^*-∇A(y)^*]∇F(A(x)) - ∇A(y)^*[∇F(A(y))-∇F(A(x))]‖ - # ≤ L_{∇A(x)} M_{∇F} + M_{∇A(y)^*} L_{∇F}L_A. + def diff_bound(self): + mf = self.f.diff_bound() if hasattr(self.op, "opnorm"): - # Linear operator - lda = 0.0 # This is zero, - mdf = 0.0 # hence this not needed. + lda = self.op.opnorm() ** 2 else: - mdf = self.f.diff_bound(xbound=self.op.codomain_bound(xbound=self.xbound)) - lda = self.op.diff_adj_lipschitz_factor() + lda = self.op.diff_bound() + return lda * mf + + def diff_bound_pair(self): + mf = self.f.diff_bound() + lda1, lda2 = self.op.diff_bound_pair() + return lda1 * mf, lda2 * mf + + # def lipschitz_factor(self, xbound=None): + # if xbound is None: + # xbound = self.xbound + # lf = self.f.lipschitz_factor(xbound=self.op.codomain_bound(xbound=xbound)) + # la = self.op.lipschitz_factor(xbound=xbound) + # return lf * la - ldf = self.f.diff_lipschitz_factor() - la = self.op.lipschitz_factor() - mda = self.op.diff_bound(xbound=self.xbound) + # def lipschitz_factor_pair(self, xbound=None): + # if xbound is None: + # xbound = self.xbound + # lf = self.f.lipschitz_factor(xbound=self.op.codomain_bound(xbound=xbound)) + # la1, la2 = self.op.lipschitz_factor_pair(xbound=xbound) + # return lf * la1, lf * la2 - return lda * mdf + mda * ldf * la + def diff_lipschitz_factor(self): + """ + Calculate the Lipschitz factor of the differential of this composed function. + We assume that either the operator is linear and implementes `opnorm`, or it is nonlinear, and impliements `diff_chain_lipschitz_factor` to directly calculate a Lipschitz factor of $x ↦ ∇A(x)^*∇F(A(x))$, given a bound $M$ + on ∇F. The function `diff_chain_lipschitz_factor` should return the Lipschitz factor divided by $M$: we obtain $M$ through the `diff_bound` on $F$. + """ + + if hasattr(self.op, "opnorm"): + return self.f.diff_lipschitz_factor() * self.op.opnorm() ** 2 + else: + mdf = self.f.diff_bound() + lda = self.op.diff_chain_lipschitz_factor() + # print( + # "LDA %s %f; MDF %f; total %f" + # % (type(self.op).__name__, lda, mdf, mdf * lda), + # ) + return mdf * lda def diff_lipschitz_factor_pair(self): - if self.op.hasattr("opnorm"): - # Linear operator - lda1, lda2 = 0.0, 0.0 # This is zero, - mdf = 0.0 # hence this not needed. - else: - lda1, lda2 = self.op.diff_adj_lipschitz_factor_pair() - mdf = self.f.diff_bound( - xbound=self.op.codomain_bound_pair(xbound=self.xbound_pair) - ) + """ + This is similar to `diff_lipschitz_factor`, except separates the + factor for arguments pairs. + + This requires the operator to implement `diff_chain_lipschitz_factor`; + there is no special handling of linear opeartors. + """ - ldf = self.f.diff_lipschitz_factor() - la1, la2 = self.op.lipschitz_factor_pair() - mda = self.op.diff_bound_pair(xbound=self.xbound_pair) + mdf = self.f.diff_bound() + + lda1, lda2 = self.op.diff_chain_lipschitz_factor_pair() - return lda1 * mdf + mda * ldf * la1, lda2 * mdf + mda * ldf * la2 + # print("LDA %s %f %f; MDF %f; " % (type(self.op).__name__, lda1, lda2, mdf)) + + return mdf * lda1, mdf * lda2 class InjectSecond: @@ -112,17 +137,17 @@ def opnorm(self, *args): return 1.0 - def lipschitz_factor(self, *args): + # def lipschitz_factor(self): + # return 1.0 + + # def diff_chain_lipschitz_factor(self, *args): + # return 0.0 + + def diff_bound(self): return 1.0 - def diff_adj_lipschitz_factor(self, *args): - return 0.0 - - def diff_bound(self, xbound=None): - return 1.0 - - def codomain_bound(self, xbound=None): - if xbound is None: - raise Exception("Linear operators have unbounded range") - else: - return xbound + # def codomain_bound(self, xbound=None): + # if xbound is None: + # raise Exception("Linear operators have unbounded range") + # else: + # return xbound
--- a/src/convection_diffusion.py Thu Feb 26 09:32:12 2026 -0500 +++ b/src/convection_diffusion.py Wed Apr 22 22:32:00 2026 -0500 @@ -78,6 +78,9 @@ beta=1.0, # \beta = \| c \|_\infty^2 /k_m C_emb=1.0, C_reg=1.0, + xbound=XBound, + override_lipschitz=None, + override_lipschitz_pair=None, ): self.p = p # L^p, p=2 self.C_stab = C_stab # adjoint stability constant @@ -95,6 +98,7 @@ self.ny = ny self.degree = degree self.T = T + self.xbound = xbound self.num_steps = num_steps self.dt = (T - t0) / num_steps @@ -102,6 +106,9 @@ self.save_for_plot = False + self.override_lipschitz = override_lipschitz + self.override_lipschitz_pair = override_lipschitz_pair + domain = mesh.create_rectangle( MPI.COMM_SELF, [ @@ -247,8 +254,8 @@ A.assemble() # Create reusable RHS vector (will be filled each timestep) - b = create_vector(linear_form) - b_ps = create_vector(linear_form) + b = create_vector(linear_form.function_spaces) + b_ps = create_vector(linear_form.function_spaces) # Prepare solver once solver = PETSc.KSP().create(domain.comm) @@ -358,7 +365,7 @@ A2.assemble() # Create vector for RHS to be updated each step - b2 = create_vector(linear_form_w) + b2 = create_vector(linear_form_w.function_spaces) solver = PETSc.KSP().create(domain.comm) solver.setOperators(A2) @@ -513,164 +520,190 @@ self.plot_u_n_list = u_n_list self.plot_w_n_list = w_n_list self.plot_j_n_list = j + self.plot_aux = (k, c) return wp_bar, (expr2, expr3) def do_plot(self, iter, μ, frames, prefix): n = self.num_steps + k, (c1, c2) = self.plot_aux np.savez_compressed( "%s/res_%d.npz" % (prefix, iter), wp_bar_array=self.plot_wp_bar.x.array, u_n_list_array=[self.plot_u_n_list[i].x.array for i in frames], w_n_list_array=[self.plot_w_n_list[i].x.array for i in frames], j_n_list_array=[self.plot_j_n_list[i].x.array for i in frames], + k=k, + c1=c1, + c2=c2, frames=frames, times=list(map(lambda i: self.T * i / (n - 1), frames)), mu=[[point[0], point[1], alpha] for point, alpha in μ.iter_padded()], ) + # print("Aux variables: k: %f, c1: %f, c2: %f" % (k, c1, c2)) - def _stability_prefactor( - self, *, C_stab=None, alpha=None, beta=None, T=None, C_emb=None, C_reg=None - ): - """ - Return the stability prefactor - P := (C_stab/alpha) * (1 + sqrt(T)*(1 + sqrt(beta)/sqrt(alpha))) - and the measure-to-functional constant L_e_mu := C_emb * C_reg * sqrt(T). - Return (P, L_e_mu): - """ - # self._stability_prefactor(alpha=0.5, T=2.0) + # def _stability_prefactor( + # self, *, C_stab=None, alpha=None, beta=None, T=None, C_emb=None, C_reg=None + # ): + # """ + # Return the stability prefactor + # P := (C_stab/alpha) * (1 + sqrt(T)*(1 + sqrt(beta)/sqrt(alpha))) + # and the measure-to-functional constant L_e_mu := C_emb * C_reg * sqrt(T). + # Return (P, L_e_mu): + # """ + # # self._stability_prefactor(alpha=0.5, T=2.0) + + # # Read from arguments or instance attributes + # C_stab = C_stab if C_stab is not None else self.C_stab + # alpha = alpha if alpha is not None else self.alpha + # beta = beta if beta is not None else self.beta + # T = T if T is not None else self.T + # C_emb = C_emb if C_emb is not None else self.C_emb + # C_reg = C_reg if C_reg is not None else self.C_reg - # Read from arguments or instance attributes - C_stab = C_stab if C_stab is not None else self.C_stab - alpha = alpha if alpha is not None else self.alpha - beta = beta if beta is not None else self.beta - T = T if T is not None else self.T - C_emb = C_emb if C_emb is not None else self.C_emb - C_reg = C_reg if C_reg is not None else self.C_reg + # # Stability factor P - # Stability factor P + # P = (C_stab / alpha) * (1 + np.sqrt(T) * (1 + np.sqrt(beta) / np.sqrt(alpha))) + # L_e_mu = C_emb * C_reg * np.sqrt(T) + + # return P, L_e_mu - P = (C_stab / alpha) * (1 + np.sqrt(T) * (1 + np.sqrt(beta) / np.sqrt(alpha))) - L_e_mu = C_emb * C_reg * np.sqrt(T) - - return P, L_e_mu - - # Solution operator Lipschitz factor - def lipschitz_factor(self): - """ - Conservative Lipschitz factor for the forward solution operator mapping - (k,c,mu) -> u. Returns a single scalar C_Lip such that - ||u2 - u1|| <= C_Lip * (||k2-k1|| + ||c2-c1|| + ||mu2-mu1||_* ) - (we use a conservative form combining the measure and (k,c) pieces). + # # Solution operator Lipschitz factor + # def lipschitz_factor(self, xbound=None): + # """ + # Conservative Lipschitz factor for the forward solution operator mapping + # (k,c,mu) -> u. Returns a single scalar C_Lip such that + # ||u2 - u1|| <= C_Lip * (||k2-k1|| + ||c2-c1|| + ||mu2-mu1||_* ) + # (we use a conservative form combining the measure and (k,c) pieces). + # """ + # P, L_e_mu = self._stability_prefactor() + # # Conservative combined choice: multiply P by max(1, L_e_mu) + # C_Lip = P * max(1.0, L_e_mu) + # return C_Lip + def diff_chain_lipschitz_factor(self, Delta_t=None): """ - P, L_e_mu = self._stability_prefactor() - # Conservative combined choice: multiply P by max(1, L_e_mu) - C_Lip = P * max(1.0, L_e_mu) - return C_Lip - - # Adjoint solution operator Lipschitz factor - """ - def diff_adj_lipschitz_factor(self): - - Compute adjoint-solution Lipschitz factor A for 2D domain. - - For p = 2, A = 1. - For p > 2, A = |Omega|^(1/2 - 1/p), where |Omega| is the mesh area. - - p = self.p - omega_vol = self.domain_volume - return 1.0 if p == 2 else omega_vol ** (0.5 - 1.0 / p) - """ - - def diff_adj_lipschitz_factor(self): - """ - Always return 1 for p=2 - """ - return 1.0 - - def diff_bound(self, xbound=None): - """ - Differential bound for the solution operator derivative. - - Returns the tuple: (L_e_mu, L_e_k, L_e_c) - - where - L_e_k = ||∇u||_{L^2(Ω_T)} - L_e_c = ||∇u||_{L^2(Ω_T)} - L_e_mu = C_emb * C_reg * sqrt(T) + Compute a Lipschitz factor C_tot^Lip using locally accumulated bounds + M_J = 1. """ - # If grad_u_norm is None, compute it from u0 - grad_u_norm = getattr(self, "grad_u_norm", None) - if grad_u_norm is None: - u = self.u0 - V = self.V - # Compute L2 norm of gradient of u0 over the domain - grad_u = ufl.grad(u) - grad_u_sq = ufl.inner(grad_u, grad_u) - dx = ufl.Measure("dx", domain=self.domain) - # Integral of |∇u|^2 over Ω - grad_u_norm = fem.assemble_scalar(fem.form(grad_u_sq * dx)) ** 0.5 - # Optionally store it for later reuse - self.grad_u_norm = grad_u_norm + if self.override_lipschitz is not None: + return self.override_lipschitz + + # Use existing diff_bound3 + C_mu_adj, C_k_adj, C_c_adj = self.diff_bound3(Delta_t=Delta_t) + + if Delta_t is None: + # fall back to the scheme used in diff_bound3 + Delta_t = self.dt + + N = int(np.ceil(self.T / Delta_t)) + + # Use local per‑step bounds + C_adj_step = C_mu_adj / N + C_state_step = C_k_adj / N + C_c_step = C_c_adj / N + + # Per‑step "Lipschitz" + l_e_local = C_adj_step + C_state_step + C_c_step + + # Global Lipschitz constants + L_e = l_e_local * N + M_e = 1.2 * L_e + L_A = 0.5 * C_state_step * N - return grad_u_norm + # PDE/coercivity + alpha = max(0.1, self.alpha) # guard α ≥ 0.1 + M_J = 1.0 + + # C_tot^Lip + inv_alpha = 1.0 / alpha + inv_alpha_sq = inv_alpha * inv_alpha + C_tot = inv_alpha * (L_e + inv_alpha_sq * L_A * M_e) * M_J + return C_tot - def diff_bound_pair(self, xbound=None, Delta_t=None): - if xbound is None: - return self.diff_bound(xbound=None) - else: - xb_combined = xbound[0] + xbound[1] if len(xbound) > 1 else xbound[0] + def diff_chain_lipschitz_factor_pair(self, Delta_t=None): + """ + Return (C^Lip_μ/M_J, C^Lip_aux/M_J) + """ + + if self.override_lipschitz_pair is not None: + return self.override_lipschitz_pair + + C_mu_adj, C_k_adj, C_c_adj = self.diff_bound3(Delta_t=Delta_t) + + if Delta_t is None: + Delta_t = self.dt + + N = int(np.ceil(self.T / Delta_t)) - k_min = getattr(xb_combined, "k_min", self.alpha) - c_Linf = getattr(xb_combined, "c_Linf", np.sqrt(self.beta * k_min)) - diam = getattr(xb_combined, "diam", self.diam) - T = getattr(xb_combined, "T", self.T) + # Extract per-component per-step bounds (smaller scaling) + C_mu_step = C_mu_adj / N + C_aux_step = max(C_k_adj, C_c_adj) / N + + # TINY Lipschitz factors + C_Lip_mu = C_mu_step * N + C_Lip_aux = C_aux_step * N - Pe = c_Linf * diam / k_min - gamma = c_Linf**2 / k_min - # adaptive time step size Delta_t (Delta_t = 0.01 in the test) - if Delta_t is None: - Delta_t = min(0.1, k_min / c_Linf**2, 0.01 * T) # CFL-like + # Light coercivity scaling + alpha = self.alpha + inv_alpha = 1.0 / alpha - # Accumulate N = T/Delta_t local steps - N = int(np.ceil(T / Delta_t)) + C_mu_tot = C_Lip_mu * inv_alpha + C_aux_tot = C_Lip_aux * inv_alpha + + return C_mu_tot, C_aux_tot + + def diff_bound3(self, Delta_t=None): + xbound = self.xbound + xb_combined = xbound # xbound[0] + xbound[1] if len(xbound) > 1 else xbound[0] - # Per-step factors - exp_local = np.exp(0.5 * gamma * Delta_t) - C_adj_step = np.sqrt(Delta_t * (exp_local + Delta_t / k_min)) - C_state_step = np.sqrt(Delta_t) * np.sqrt(1 + Pe**2) + k_min = getattr(xb_combined, "k_min", self.alpha) + c_Linf = getattr(xb_combined, "c_Linf", np.sqrt(self.beta * k_min)) + diam = getattr(xb_combined, "diam", self.diam) + T = getattr(xb_combined, "T", self.T) + + Pe = c_Linf * diam / k_min + gamma = c_Linf**2 / k_min + # adaptive time step size Delta_t (Delta_t = 0.01 in the test) + if Delta_t is None: + Delta_t = min(0.1, k_min / c_Linf**2, 0.01 * T) # CFL-like + + # Accumulate N = T/Delta_t local steps + N = int(np.ceil(T / Delta_t)) - # Total accumulated bound (geometric sum <= N * max_step) - C_adj_PDE = N * C_adj_step - C_state = N * C_state_step + # Per-step factors + exp_local = np.exp(0.5 * gamma * Delta_t) + C_adj_step = np.sqrt(Delta_t * (exp_local + Delta_t / k_min)) + C_state_step = np.sqrt(Delta_t) * np.sqrt(1 + Pe**2) - # Use the embedding constant - C_mu_adj = self.C_emb * C_adj_PDE # ∫φ δμ - C_k_adj = C_state * C_adj_PDE # ∫∇u·∇φ δk - C_c1_adj = C_state * C_adj_PDE * c_Linf # ∫∂x u φ δc1 - C_c2_adj = C_state * C_adj_PDE * c_Linf # ∫∂y u φ δc2 + # Total accumulated bound (geometric sum <= N * max_step) + C_adj_PDE = N * C_adj_step + C_state = N * C_state_step - return C_mu_adj, C_k_adj, C_c1_adj, C_c2_adj + # Use the embedding constant + C_mu_adj = self.C_emb * C_adj_PDE # ∫φ δμ + C_k_adj = C_state * C_adj_PDE # ∫∇u·∇φ δk + C_c_adj = C_state * C_adj_PDE * c_Linf # ∫∂x u φ δc1 and c2 + + return C_mu_adj, C_k_adj, C_c_adj # ||S(x)||_Y→X ≤ C_state ||μ||_M(Ω) [time-independent μ] - def codomain_bound(self, xbound=None): - if xbound is None: - return 1.0 + def codomain_bound(self): + xbound = self.xbound # Extract parameters if hasattr(xbound, "k_min"): xb = xbound + elif hasattr(xbound, "__len__") and len(xbound) > 0: + xb = xbound[0] + xbound[1] if len(xbound) > 1 else xbound[0] else: - if hasattr(xbound, "__len__") and len(xbound) > 0: - xb = xbound[0] + xbound[1] if len(xbound) > 1 else xbound[0] - else: - xb = xbound - k_min = getattr(xb, "k_min", self.alpha) - c_Linf = getattr(xb, "c_Linf", np.sqrt(self.beta * k_min)) - mu_dual = getattr(xb, "mu_dual", 3.0) - diam = getattr(xb, "diam", self.diam) - T = getattr(xb, "T", self.T) + xb = xbound + + k_min = getattr(xb, "k_min", self.alpha) + c_Linf = getattr(xb, "c_Linf", np.sqrt(self.beta * k_min)) + mu_dual = getattr(xb, "mu_dual", 3.0) + diam = getattr(xb, "diam", self.diam) + T = getattr(xb, "T", self.T) Pe = c_Linf * diam / k_min if not hasattr(self, "Delta_t"): @@ -687,88 +720,6 @@ return C_state * mu_dual - # ||S(x1)-S(x2)|| ≤ Cμ||Δμ|| + Ck||Δk|| + Cc1||Δc1|| + Cc2||Δc2|| - - def codomain_bound_pair(self, xbound=None): - if xbound is None: - return self.codomain_bound(xbound=None) - else: - xb = xbound[0] + xbound[1] if len(xbound) > 1 else xbound[0] - C_base = self.codomain_bound(xb) - # Component - C_mu = C_base - C_k = C_base * 2.0 - C_c1 = C_base * 1.5 - C_c2 = C_base * 1.5 - - return (C_mu, C_k, C_c1, C_c2) - - """ - def codomain_bound(self, xbound=None): - return 1.0 - - def codomain_bound_pair(self, xbound=None): - # TODO: implement properly - if xbound is None: - return self.codomain_bound(xbound=None) - else: - return self.codomain_bound(xbound=xbound[0] + xbound[1]) - """ - - def diff_bound3(self): - """ - Returns full operator differential bounds: - - L_e_mu = C_emb * C_reg * sqrt(T) - L_e_k = ||∇u||_{L^2(Ω)} - L_e_c = ||∇u||_{L^2(Ω)} - - """ - - # First get L_e_k = L_e_c - grad_u_norm = self.diff_bound() - - # Compute L_e_mu from stability prefactor - _P, L_e_mu = self._stability_prefactor() - - L_e_k = grad_u_norm - L_e_c = grad_u_norm - - return L_e_mu, L_e_k, L_e_c - - # Solution operator Lipschitz factor separately wrt. μ and (k, c) - def lipschitz_factor_pair(self): - """ - Compute the forward solution operator Lipschitz factors - separately with respect to the parameters (k, c) and μ. - - Returns: - tuple: (L_mu, L_kc) - L_mu - Lipschitz factor w.r.t. μ - L_kc - Lipschitz factor w.r.t. (k, c) - """ - # Get differential bounds for the solution operator - L_e_mu, L_e_k, L_e_c = self.diff_bound3() - - # Lipschitz factor w.r.t μ - L_mu = L_e_mu - - # Lipschitz factor w.r.t (k, c) - # Conservative choice: max of L_e_k and L_e_c - L_kc = max(L_e_k, L_e_c) - - return L_mu, L_kc - - # Adjoint solution operator Lipschitz factor separately wrt. μ and (k, c) - def diff_adj_lipschitz_factor_pair(self): - A = self.diff_adj_lipschitz_factor() - - # No separate analytic dependence, so return same A for both parts - A_mu = A - A_kc = A - - return A_mu, A_kc - def own(u): # """ @@ -853,20 +804,22 @@ def apply(self, x): """ - Apply the quadratic regularisation fucntional to `x`. + Apply the quadratic regularisation functional to `x`. """ (k, (c1, c2)) = x (kb, (c1b, c2b)) = self._get_base() (nkb, (nc1b, nc2b)) = self.base_norm_squared - return self.α * ( - self._apply1(k, kb, nkb) - + self._apply1(c1, c1b, nc1b) - + self._apply1(c2, c2b, nc2b) + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + + return αk * self._apply1(k, kb, nkb) + αc * ( + self._apply1(c1, c1b, nc1b) + self._apply1(c2, c2b, nc2b) ) - def _prox1(self, τ, x, b): - γ = self.α * τ + def _prox1(self, γ, x, b): β = 1.0 / (1.0 + γ) if isinstance(x, float): return β * (x if b is None else x + b * γ) @@ -880,7 +833,7 @@ def prox(self, τ, x): """ - Apply the proximal map of the quadratic regularisation fucntional to `x`. + Apply the proximal map of the quadratic regularisation functional to `x`. WARNING: This function is unsafe. It may modify `x´. That is ok for our purposes, as Rust, being a safe language, already needs to pass a copied or owned instance to Python. @@ -888,7 +841,166 @@ k, (c1, c2) = x (kb, (c1b, c2b)) = self._get_base() + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + return ( - self._prox1(τ, k, kb), - (self._prox1(τ, c1, c1b), self._prox1(τ, c2, c2b)), + self._prox1(τ * αk, k, kb), + (self._prox1(τ * αc, c1, c1b), self._prox1(τ * αc, c2, c2b)), + ) + + +class LogPowerRegularisation: + """ + Logarithmic barrier + quadratic regularisation functional for the auxiliary variables; + $-α\\log(x) + ωx^2$. + """ + + def _own(self, x): + return x if isinstance(x, float) else own(x) + + def __init__(self, α, ω): + self.α = α + self.ω = ω + + def _apply1(self, x, α, ω): + if isinstance(x, float): + return -α * np.log(x) + ω * x**2 + else: + # Numpy is Matlab-grade junk. No reduce! Temporary arrays! In 2026! WTF?!?!?!? + return (-α * np.log(x) + ω * x**2).sum() + + def apply(self, x): + (k, (c1, c2)) = x + + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + if isinstance(self.ω, float): + ωk, ωc = self.ω, self.ω + else: + ωk, ωc = self.ω + + return ( + self._apply1(k, αk, ωk) + + self._apply1(c1, αc, ωc) + + self._apply1(c2, αc, ωc) + ) + + def _prox1(self, α, ω, x): + # OC: 0 = -α/z + 2ωz + z-x ⟺ -α + (2ω+1)z^2 - xz = 0. + # ⟺ z = -x ± √(x^2 + 4α(2ω+1))/(2(2ω+1)) + # Numpy is Matlab-grade junk. No reduce! Temporary arrays! In 2026! WTF?!?!?!? + return (x + np.sqrt(x**2 + 4 * α * (2 * ω + 1))) / (2 * (2 * ω + 1)) + + def prox(self, τ, x): + """ + Apply the proximal map of the logarithmic regularisation functional to `x`. + WARNING: This function is unsafe. It may modify `x´. + That is ok for our purposes, as Rust, being a safe language, already needs to pass a + copied or owned instance to Python. + """ + k, (c1, c2) = x + + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + if isinstance(self.ω, float): + ωk, ωc = self.ω, self.ω + else: + ωk, ωc = self.ω + + return ( + self._prox1(τ * αk, τ * ωk, k), + (self._prox1(τ * αc, τ * ωc, c1), self._prox1(τ * αc, τ * ωc, c2)), ) + + +class BoxConstraints: + """ + Simple box constraints on the auxiliary variables + """ + + def _own(self, x): + return x if isinstance(x, float) else own(x) + + def __init__(self, α, ω): + self.α = α + self.ω = ω + + def _apply1(self, x, α, ω): + if isinstance(x, float): + if α <= x and x <= ω: + return 0.0 + else: + return np.inf + else: + raise Exception("Unimplemented") + + def apply(self, x): + (k, (c1, c2)) = x + + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + if isinstance(self.ω, float): + ωk, ωc = self.ω, self.ω + else: + ωk, ωc = self.ω + + return ( + self._apply1(k, αk, ωk) + + self._apply1(c1, αc, ωc) + + self._apply1(c2, αc, ωc) + ) + + def _prox1(self, α, ω, x): + # OC: 0 = -α/z + 2ωz + z-x ⟺ -α + (2ω+1)z^2 - xz = 0. + # ⟺ z = -x ± √(x^2 + 4α(2ω+1))/(2(2ω+1)) + # Numpy is Matlab-grade junk. No reduce! Temporary arrays! In 2026! WTF?!?!?!? + return max(min(ω, x), α) + + def prox(self, τ, x): + """ + Apply the proximal map of the box constraints to `x`. + WARNING: This function is unsafe. It may modify `x´. + That is ok for our purposes, as Rust, being a safe language, already needs to pass a + copied or owned instance to Python. + """ + k, (c1, c2) = x + + if isinstance(self.α, float): + αk, αc = self.α, self.α + else: + αk, αc = self.α + if isinstance(self.ω, float): + ωk, ωc = self.ω, self.ω + else: + ωk, ωc = self.ω + + return ( + self._prox1(αk, ωk, k), + (self._prox1(αc, ωc, c1), self._prox1(αc, ωc, c2)), + ) + + +""" +Quadratic regularisation with box contraints +""" + + +class BoxedQuadraticRegularisation: + def __init__(self, ω0, ω1, α, base=None): + self.box = BoxConstraints(ω0, ω1) + self.quadratic = QuadraticRegularisation(α, base) + + def apply(self, x): + return self.box.apply(x) + self.quadratic.apply(x) + + def prox(self, τ, x): + return self.box.prox(1.0, self.quadratic.prox(τ, x))
--- a/src/dolfinx_access.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/dolfinx_access.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,3 +1,6 @@ +/*! + Utility functions accessing C++ side Dolfinx. +*/ use alg_tools::error::DynResult; use alg_tools::fe_model::{ base::RealLocalModel, @@ -9,7 +12,7 @@ use std::ffi::CString; use std::mem::MaybeUninit; -// These are required for the linking to the sys crates to actually happen. +/// Import helper required for the linking to the sys crates to actually happen. #[allow(unused_imports)] mod dummy_import { use dolfinx_sys; @@ -19,20 +22,29 @@ mod function; pub use function::DolfinxPyFunction_f64; +/// C++ side routines #[allow(dead_code)] #[cxx::bridge(namespace = "dolfinx_access")] pub mod ffi { + /// Structure for coordinte-value pairs #[derive(Copy, Clone, Debug)] struct CoordValuePair { + /// 2D coordinate x: [f64; 2], + /// Value v: f64, } + /// Structure for information about a dolfinx `Function` #[derive(Copy, Clone, Debug)] struct FunctionInfo { + /// Domain dimension domain_dim: u32, + /// Codomain dimension codomain_dim: u32, + /// Element order order: u32, + /// Whether the mesh is triangular triangular_mesh: bool, } @@ -41,7 +53,9 @@ include!("pointsource_pde/include/dolfinx_access/function.h"); include!("pointsource_pde/include/dolfinx_access/minmax_p2.h"); + /// C++ side Dolfinx `Function<double>` type Function_f64; + /// C-side python object. type PyObject; /// Find the cell containing `x` (in 1D, 2D, or 3D, following Fenics weirdness @@ -66,18 +80,34 @@ cell: i32, ) -> [f64; 6]; + /// Get information about a dolfinx `Function_f64`. pub unsafe fn info_Function_f64(f: *const Function_f64) -> FunctionInfo; + /// Returns the minimum value of the Dolfinx `Function`, + /// if supported (P2 elements on triangular mesh) pub unsafe fn min_Function_f64_p2(f: *const Function_f64) -> Result<CoordValuePair>; + + /// Returns the maximum value of the Dolfinx `Function`, + /// if supported (P2 elements on triangular mesh) pub unsafe fn max_Function_f64_p2(f: *const Function_f64) -> Result<CoordValuePair>; + + /// Returns the minimum or maximum value of the Dolfinx `Function`, + /// if supported (P2 elements on triangular mesh) pub unsafe fn minmax_Function_f64_p2( f: *const Function_f64, max: bool, ) -> Result<CoordValuePair>; + /// Checks if a Dolfinx `Function` is supported (P2 elements on triangular mesh) pub unsafe fn check_Function_f64(o: *const PyObject) -> bool; + + /// Casts a Python object into a Dolfinx `Function` wrapper. pub unsafe fn cast_Function_f64(o: *const PyObject) -> Result<*const Function_f64>; + + /// Casts a Python object into a mutable Dolfinx `Function` wrapper. pub unsafe fn cast_mut_Function_f64(o: *mut PyObject) -> Result<*mut Function_f64>; + + /// Gets the 3D coordinates of the nodes of a triangular cell pub unsafe fn cell_coords_Function_f64_triangle( f: *const Function_f64, cell: i32, @@ -95,18 +125,28 @@ g: *const Function_f64, ) -> Result<()>; - /// Create a new similar function + /// Create a new similar Dolfinx `Function` pub unsafe fn similar_Function_f64(f: *const Function_f64) -> *mut Function_f64; - /// Clone the function + /// Clone the Dolfinx `Function` pub unsafe fn clone_Function_f64(f: *const Function_f64) -> *mut Function_f64; + /// Wraps a Dolfinx `Function` into a Python object pub unsafe fn wrap_Function_f64(f: *mut Function_f64) -> *mut PyObject; + /// Drop helper pub unsafe fn drop_Function_f64(f: *mut Function_f64); + + /// Scatter helper. Not truly supported. + pub unsafe fn scatter_fwd_Function_f64(f: *mut Function_f64); } extern "Rust" { + /// Rust-side helper for minimising P2 elements. This bypasses poorly documented basix, + /// and does things at a low level directly. + /// `f` is the function, `simplex_coords` the coorindates of the triangle + /// see [`cell_coords_Function_f64_triangle`], `cell` is the cell number, + /// and `max` indicates whether to maximise or minimise. unsafe fn minmax_dolfinx_p2_cell( f: *const Function_f64, simplex_coords: &[f64; 9],
--- a/src/dolfinx_access/function.cc Thu Feb 26 09:32:12 2026 -0500 +++ b/src/dolfinx_access/function.cc Wed Apr 22 22:32:00 2026 -0500 @@ -17,7 +17,6 @@ using namespace basix::element; using namespace dolfinx::graph; using namespace dolfinx::la; -namespace cell = basix::cell; namespace dolfinx_access { extern void drop_Function_f64(Function_f64* f) { @@ -29,9 +28,9 @@ auto el = fp->element()->basix_element(); return FunctionInfo{ .domain_dim = (uint32_t)fp->mesh()->geometry().dim(), - .codomain_dim = (uint32_t)fp->value_size(), + .codomain_dim = (uint32_t)fp->element()->value_size(), .order = (uint32_t)fp->element()->basix_element().degree(), - .triangular_mesh = el.cell_type() == cell::type::triangle, + .triangular_mesh = el.cell_type() == basix::cell::type::triangle, }; } @@ -41,15 +40,16 @@ } } - std::map<std::weak_ptr<const mesh::Mesh<double>>, std::shared_ptr<BoundingBoxTree<double>>, - std::owner_less<std::weak_ptr<const mesh::Mesh<double>>>> + std::map<std::weak_ptr<const dolfinx::mesh::Mesh<double>>, + std::shared_ptr<BoundingBoxTree<double>>, + std::owner_less<std::weak_ptr<const dolfinx::mesh::Mesh<double>>>> bb_tree_cache; /// Returns the cell containg x - int32_t cell_Mesh_f64(std::shared_ptr<const mesh::Mesh<double>> mesh, + int32_t cell_Mesh_f64(std::shared_ptr<const dolfinx::mesh::Mesh<double>> mesh, const std::array<double, 3>& x) { - std::weak_ptr<const mesh::Mesh<double>> mesh_weak = mesh; + std::weak_ptr<const dolfinx::mesh::Mesh<double>> mesh_weak = mesh; std::shared_ptr<BoundingBoxTree<double>> bb_tree = bb_tree_cache[mesh_weak]; if (bb_tree == nullptr) { auto topo = mesh->topology(); @@ -68,7 +68,8 @@ assert(local_cells == index_map->size_global()); std::vector<std::int32_t> r(local_cells); std::iota(r.begin(), r.end(), 0); - bb_tree = std::make_shared<BoundingBoxTree<double>>(BoundingBoxTree(*mesh, dim, r)); + bb_tree = + std::make_shared<BoundingBoxTree<double>>(BoundingBoxTree(*mesh, dim, 0.0, r)); bb_tree_cache[mesh_weak] = bb_tree; } // Find colliding bounding boxes. This is an AdjancencyTree. @@ -143,13 +144,13 @@ } rust::Slice<const double> data_Function_f64(Function_f64 const* f) { - auto span = f->x()->array(); - return rust::Slice(span.data(), span.size()); + const std::vector<double>& v = f->x()->array(); + return rust::Slice<const double>(v.data(), v.size()); } rust::Slice<double> data_mut_Function_f64(Function_f64* f) { - auto span = f->x()->mutable_array(); - return rust::Slice(span.data(), span.size()); + std::vector<double>& v = f->x()->array(); + return rust::Slice<double>(v.data(), v.size()); } Function_f64* similar_Function_f64(Function_f64 const* f) { @@ -157,7 +158,20 @@ } Function_f64* clone_Function_f64(Function_f64 const* f) { - return new Function_f64(f->function_space(), std::make_shared<Vector<double>>(*f->x())); + // return new Function_f64(f->function_space(), + // std::make_shared<Vector<double>>(*f->x())); + auto f_copy = new Function_f64(f->function_space()); + auto fx = f->x(); + auto x_new = f_copy->x(); + const std::vector<double>& a_old = fx->array(); + std::vector<double>& a_new = x_new->array(); + std::copy(a_old.begin(), a_old.end(), a_new.begin()); + x_new->scatter_fwd(); + return f_copy; + } + + void scatter_fwd_Function_f64(Function_f64* f) { + f->x()->scatter_fwd(); } // PyObject* py_similar_Function_f64(Function_f64 const* f) {
--- a/src/dolfinx_access/function.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/dolfinx_access/function.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,5 +1,5 @@ /*! -Python and C++ wrapper for Dolfinx Function<f64> +Python and C++ wrapper for Dolfinx `Function<f64>`. */ use super::ffi; @@ -43,15 +43,21 @@ cxx: *mut Function_f64, } +/// References to Python side helper routines struct Helpers { + /// Dot product dot: Py<PyFunction>, + /// Squared 2-norm norm2_squared: Py<PyFunction>, + /// Squared distance dist2_squared: Py<PyFunction>, } use std::sync::OnceLock; +/// A populated-once structure for referencing Python-side helper routines. static HELPERS: OnceLock<Helpers> = OnceLock::new(); +/// Gets [`Helpers`], populates [`HELPERS`] if needed. fn get_helpers(py: Python<'_>) -> PyResult<&Helpers> { HELPERS.get_or_try_init(|| { let extras = PyModule::import(py, "pointsource_pde.dolfinx_extras")?; @@ -72,6 +78,7 @@ unsafe { ffi::drop_Function_f64(self.cxx) }; } } + self.cxx = std::ptr::null_mut(); } } @@ -80,7 +87,7 @@ DolfinxPyFunction_f64 { u: None, // It's a different new function, so no Python instance function_space: self.function_space.clone(), - owned: self.owned, + owned: true, cxx: unsafe { ffi::clone_Function_f64(self.cxx) }, } } @@ -135,7 +142,7 @@ /// Returns a (possibly new) Python presentation, without consuming current object. /// - /// Please use this functionality through [`pyo3::types::IntoPyObject`]. + /// Please use this functionality through [`pyo3::IntoPyObject`]. pub(crate) fn _into_py_ref<'a>( &'a mut self, py: Python<'py>, @@ -150,7 +157,7 @@ /// Returns a (possibly new) Python presentation, consuming current object. /// - /// Please use this functionality through [`pyo3::types::IntoPyObject`]. + /// Please use this functionality through [`pyo3::IntoPyObject`]. pub(crate) fn _into_py(mut self, py: Python<'py>) -> PyResult<pyo3::Bound<'py, PyAny>> { self._make_py(py); self.cxx = std::ptr::null_mut(); @@ -181,6 +188,9 @@ } } +/// An extremely durty hack to create a Python-side Dolfinx `Function` referencing +/// a C++-side Dolfinx `Function` that we already have. For some unfathomable reason, +/// Dolfinx provies no way to do this cleanly. pub(crate) static CREATE_FUNCTION_HACK: &CStr = c_str!( "\ from dolfinx.fem import Function @@ -198,9 +208,12 @@ " ); +/// Trait for points that can be expanded into Dolfinx' always-3D presentation trait ExpandCoordDolphinx { + /// Convert to always-3D presentation fn to_dolphinx(self) -> [f64; 3]; #[allow(dead_code)] + /// Convert from always-3D presentation fn from_dolphinx(val: [f64; 3]) -> Self; } @@ -243,6 +256,7 @@ } } +/// Implements [`DifferentiableImpl`] macro_rules! impl_diff { ($n:literal ; $($o:literal)+) => { $( impl<'py> DifferentiableImpl<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, $o, 1> { @@ -274,6 +288,7 @@ )+ }; } +/// Implements [`Mapping`] macro_rules! impl_mapping { ($($n:literal)+) => { $( impl<'py, const O: u32> Mapping<Loc<$n, f64>> for DolfinxPyFunction_f64<'py, $n, O, 1> { @@ -342,6 +357,17 @@ */ impl<'py, const N: u32, const O: u32, const D: u32> DolfinxPyFunction_f64<'py, N, O, D> { + fn new_cxx(&mut self, new: *mut Function_f64) { + if let None = self.u { + if self.cxx != std::ptr::null_mut() { + unsafe { ffi::drop_Function_f64(self.cxx) }; + } + } + self.cxx = new; + self.u = None; + self.owned = true; + } + /// If owned, call `f_valid` with mutable data. /// If not owned, create a new C++ fem.Function, uninitialised, and call `f_invalid` with its /// mutable data, as well as the original immutable data. @@ -352,17 +378,18 @@ G: Fn(&mut [f64], &[f64]) -> Z, { if self.owned { - return f_valid(unsafe { ffi::data_mut_Function_f64(self.cxx) }); + let res = f_valid(unsafe { ffi::data_mut_Function_f64(self.cxx) }); + unsafe { ffi::scatter_fwd_Function_f64(self.cxx) }; + return res; } else { // Create a completely new Function let old_data = unsafe { ffi::data_mut_Function_f64(self.cxx) }; let new = unsafe { ffi::similar_Function_f64(self.cxx) }; let new_data = unsafe { ffi::data_mut_Function_f64(new) }; let res = f_invalid(new_data, old_data); + unsafe { ffi::scatter_fwd_Function_f64(new) }; // Invalidate Python reference; we're managing the object now - self.cxx = new; - self.u = None; - self.owned = true; + self.new_cxx(new); return res; } } @@ -376,16 +403,17 @@ F: Fn(&mut [f64]) -> Z, { if self.owned { - return f(unsafe { ffi::data_mut_Function_f64(self.cxx) }); + let res = f(unsafe { ffi::data_mut_Function_f64(self.cxx) }); + unsafe { ffi::scatter_fwd_Function_f64(self.cxx) }; + return res; } else { // Create a completely new Function let new = unsafe { ffi::similar_Function_f64(self.cxx) }; let new_data = unsafe { ffi::data_mut_Function_f64(new) }; let res = f(new_data); + unsafe { ffi::scatter_fwd_Function_f64(new) }; // Invalidate Python reference; we're managing the object now - self.cxx = new; - self.u = None; - self.owned = true; + self.new_cxx(new); return res; } } @@ -404,10 +432,11 @@ let old_data = unsafe { ffi::data_Function_f64(self.cxx) }; let new_data = unsafe { ffi::data_mut_Function_f64(new.cxx) }; f(new_data, old_data); + unsafe { ffi::scatter_fwd_Function_f64(new.cxx) }; new } - /// Similar to [`with_mut_data`], but check compatibility with `other` and also pass its + /// Similar to [`Self::with_mut_data`], but check compatibility with `other` and also pass its /// immutable data slice as the final argument. #[inline] fn with_mut_data_and_orig_binop<F, G, Z>( @@ -424,7 +453,7 @@ ffi::check_compat_Function_f64(self.cxx, other.cxx).unwrap(); ffi::data_Function_f64(other.cxx) }; - self.with_mut_data_and_orig( + let res = self.with_mut_data_and_orig( |data_self_valid| { assert_eq!(data_self_valid.len(), data_other.len()); f_valid(data_self_valid, data_other) @@ -434,10 +463,12 @@ assert_eq!(data_orig.len(), data_other.len()); f_invalid(data_self_invalid, data_orig, data_other) }, - ) + ); + unsafe { ffi::scatter_fwd_Function_f64(self.cxx) }; + return res; } - /// Similar to [`with_new_data_and_orig`], but check compatibility with `other` and also + /// Similar to [`Self::with_new_data_and_orig`], but check compatibility with `other` and also /// pass its immutable data slice as the final argument. #[inline] fn with_new_data_and_orig_binop<F>( @@ -460,6 +491,7 @@ } } +/// Implements unary operations macro_rules! impl_unop { ($trait:ident, $fn:ident) => { impl<'py, const N: u32, const O: u32, const D: u32> $trait @@ -500,6 +532,7 @@ impl_unop!(Neg, neg); +/// Implements scalar operations macro_rules! impl_scalarop { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl<'py, const N: u32, const O: u32, const D: u32> $trait<f64> @@ -564,6 +597,7 @@ impl_scalarop!(Mul, mul, MulAssign, mul_assign); impl_scalarop!(Div, div, DivAssign, div_assign); +/// Implements consuming binary operations macro_rules! impl_binop_consume { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl<'py, const N: u32, const O: u32, const D: u32> $trait<Self> @@ -635,6 +669,7 @@ }; } +/// Implements binary operations on references macro_rules! impl_binop_ref { ($trait:ident, $fn:ident) => { impl<'py, 'b, const N: u32, const O: u32, const D: u32> $trait<Self> @@ -665,6 +700,7 @@ }; } +/// Implements binary operations, both consuming and on references macro_rules! impl_binop { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl_binop_consume!($trait, $fn, $trait_assign, $fn_assign); @@ -677,9 +713,7 @@ impl_binop!(Add, add, AddAssign, add_assign); impl_binop!(Sub, sub, SubAssign, sub_assign); -//use alg_tools::euclidean::Euclidean; use alg_tools::linops::{VectorSpace, AXPY}; -//use alg_tools::norms::{Dist, HasDual, Norm, Normed, L2}; use alg_tools::self_ownable; self_ownable!(DolfinxPyFunction_f64<'py, N, O, D> where 'py, const N: u32, const O: u32, const D: u32);
--- a/src/dolfinx_access/minmax_p2.cc Thu Feb 26 09:32:12 2026 -0500 +++ b/src/dolfinx_access/minmax_p2.cc Wed Apr 22 22:32:00 2026 -0500 @@ -12,7 +12,6 @@ #include "pointsource_pde/src/dolfinx_access.rs.h" using namespace dolfinx::fem; -namespace cell = basix::cell; using namespace basix::element; using namespace dolfinx_access; @@ -45,17 +44,17 @@ auto el = fp->element()->basix_element(); // printf("%d %d %d %d\n", el.degree(), el.cell_type(), el.lagrange_variant(), // el.family()); - if (el.degree() != 2 || el.cell_type() != cell::type::triangle || + if (el.degree() != 2 || el.cell_type() != basix::cell::type::triangle || /*el.lagrange_variant() != lagrange_variant::equispaced ||*/ el.family() != family::P) { throw "Only equispaced Lagrange second-order polynomial elements are supported"; } - if (cell_num_entities(mesh::CellType::triangle, 0) != 3) { + if (cell_num_entities(dolfinx::mesh::CellType::triangle, 0) != 3) { throw "A triangle should have three vertices"; } // Check that we are dealing with a scalar function - if (fp->value_shape().size() != 0) { + if (fp->element()->value_shape().size() != 0) { throw "Only scalar functions are supported, obviously."; } @@ -67,7 +66,7 @@ // Check that there are no other types of eleemnts auto entity_types = topo->entity_types(2); for (auto& t : entity_types) { - if (t != mesh::CellType::triangle) { + if (t != dolfinx::mesh::CellType::triangle) { throw "Only triangular meshes are supported"; } }
--- a/src/experiments.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/experiments.rs Wed Apr 22 22:32:00 2026 -0500 @@ -217,12 +217,13 @@ PythonPlotFactory<'py, DerivativeCodomain<'py>, DerivativeCodomain<'py>>; } +/// Available regularisation terms, exported to the Python side #[pyclass(module = "pointsource_pde", name = "RegTerm")] #[derive(Copy, Clone, Debug, Serialize, Deserialize, PartialEq)] pub enum RegTermPy { - // Radon norm with weight `α`. + /// Radon norm with weight `α`. Radon(f64), - // Radon norm ith weight `α` and a positivity constraint. + /// Radon norm ith weight `α` and a positivity constraint. NonnegRadon(f64), } @@ -235,32 +236,32 @@ } } +/// Problem description, exported to the Python side to be filled there. #[pyclass(module = "pointsource_pde")] #[derive(Debug)] pub struct Problem { /// Data term. On the python side, this should be a `class` that implements - /// `apply` from [`DiscreteMeasure_2_f64`] (TODO: extended parameters) to floats, and `diff` - /// from the space space to [`DolfinxPyFunction_f64<2, 1, 2>`]. + /// `apply` and `diff`. #[pyo3(set)] - dataterm: Py<PyAny>, //exp_f64_2::DataTerm, + dataterm: Py<PyAny>, - // Regularisation + /// Regularisation #[pyo3(set, get)] regterm: RegTermPy, - // Auxiliary variable regularisation term or similar + /// Auxiliary variable regularisation term or similar #[pyo3(set, get)] auxterm: Option<Py<PyAny>>, - // Initial auxiliary variable + /// Initial auxiliary variable #[pyo3(set, get)] auxinit: Option<Py<PyAny>>, - // Initial measure + /// Initial measure #[pyo3(set, get)] μinit: Option<DiscreteMeasure<Loc<2, f64>, f64>>, - // Plotter + /// Plotter #[pyo3(set, get)] plot_factory: Option<Py<PyAny>>, }
--- a/src/full_sampling.py Thu Feb 26 09:32:12 2026 -0500 +++ b/src/full_sampling.py Wed Apr 22 22:32:00 2026 -0500 @@ -76,7 +76,7 @@ def lipschitz_factor(self, *args): return self.opnorm(*args) - def diff_adj_lipschitz_factor(self, *args): + def diff_chain_lipschitz_factor(self, *args): return self.opnorm(*args) def diff_bound(self, *args):
--- a/src/laser_sampling.py Thu Feb 26 09:32:12 2026 -0500 +++ b/src/laser_sampling.py Wed Apr 22 22:32:00 2026 -0500 @@ -1,13 +1,12 @@ import numpy as np from dolfinx import fem, geometry +from dolfinx_access import cell_FunctionSpace_f64 from mpi4py import MPI from petsc4py import PETSc from pointsource_pde.dolfinx_extras import get_mass_matrix from scipy.sparse.linalg import svds from slepc4py import SLEPc -from dolfinx_access import cell_FunctionSpace_f64 - class LaserSampling: """ @@ -52,14 +51,14 @@ E.setFromOptions() E.solve() - sigma = abs(E.getEigenvalue(0)) + # sigma = abs(E.getEigenvalue(0)) # self._opnorm = 1.0 / sigma # print("Full sampling opnorm = %f" % self._opnorm) # SVD operator norm _u, s, _v = svds(self.R, k=1) self._opnorm = s[0] # / sigma - print(f"Laser sampling opnorm = {self._opnorm:.6f}") + # print(f"Laser sampling opnorm = {self._opnorm:.6f}") def matrix(self, beams): """Combined: generate beams + segment + cell collision + dof weighting""" @@ -161,21 +160,17 @@ def lipschitz_factor(self, *args): return self.opnorm(*args) - def diff_adj_lipschitz_factor(self, *args): + def diff_chain_lipschitz_factor(self, *args): return self.opnorm(*args) def diff_bound(self, *args): return self.opnorm(*args) - def codomain_bound(self, xbound=None): - if xbound is None: - raise Exception("Linear operators have unbounded range") - else: - return self.opnorm() * xbound - # Construct observation noise - def noise(self, noise_level): + def noise(self, noise_level, rng=None): + if rng is None: + rng = np.random.default_rng() M = self.M # Generate noise with shape (M, 1) using scalar standard deviation - noise = np.random.normal(0, noise_level, size=(M,)) + noise = rng.normal(0, noise_level, size=(M,)) return noise
--- a/src/main.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/main.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,7 +1,5 @@ -/*! -TODO: include README here. -*/ - +// The main documentation is in README.md +#![doc = include_str!(concat!(env!("CARGO_MANIFEST_DIR"), "/README.md"))] #![feature(maybe_uninit_array_assume_init)] #![feature(iterator_try_collect)] #![feature(once_cell_try)] @@ -27,6 +25,7 @@ use experiments::Experiments; use python_access::pymod_pointsource_pde; +/// Helper for construct the list Python modules that we include as strings from `.py` files. macro_rules! pymods { [$($modname:expr),*] => {&[$( (c_str!(concat!("pointsource_pde.", $modname)), @@ -35,6 +34,7 @@ )*]}; } +/// List of Python modules that we include as strings from `.py` files. const PY_MODULES: &[(&CStr, &CStr, &CStr)] = pymods![ "dolfinx_extras", "compose", @@ -47,6 +47,11 @@ /// The entry point for the program. pub fn main() -> DynResult<()> { + unsafe { + // "#(/€"#€/("#€(/ OpenMPI and Fenics junk, don't do your obsolete “single-program multiple-data, with no controller at all” threading. + std::env::set_var("OMP_NUM_THREADS", "1"); + } + // Initialise logging colog::init(); @@ -74,8 +79,8 @@ import mpi4py mpi4py.rc.initialize = False # do not initialize MPI automatically -mpi4py.rc.thread_level = "multiple" -mpi4py.rc.threads = True +mpi4py.rc.thread_level = "serialized" +mpi4py.rc.threads = False from mpi4py import MPI
--- a/src/python_access.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/python_access.rs Wed Apr 22 22:32:00 2026 -0500 @@ -48,7 +48,7 @@ Ok(()) } -/// Converts a [`PyErr`] error in a [`Results`] error into a string with traceback. +/// Converts a [`pyo3::PyErr`] error in a [`pyo3::PyResult`] into a string with traceback. pub fn process_error<'py, T>(fname: &str, py: Python<'py>, res: PyResult<T>) -> DynResult<T> where {
--- a/src/python_access/diff_mapping.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/python_access/diff_mapping.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,5 +1,5 @@ /*! -[`DifferentiableMapping`]s implemented in Python. +[`alg_tools::mapping::DifferentiableMapping`]s implemented in Python. */ use super::{process_error, NumpyArray_f64}; use crate::dolfinx_access::DolfinxPyFunction_f64;
--- a/src/python_access/numpy_array.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/python_access/numpy_array.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,5 +1,5 @@ /*! -Python and C++ wrapper for Dolfinx Function<f64> +Wrapper for numpy arrays. */ use alg_tools::euclidean::wrap::{WrapGuard, WrapGuardMut, Wrapped};
--- a/src/python_access/plot.rs Thu Feb 26 09:32:12 2026 -0500 +++ b/src/python_access/plot.rs Wed Apr 22 22:32:00 2026 -0500 @@ -1,3 +1,6 @@ +/*! + Access to Python-side plotter implementations. +*/ use super::process_error; use alg_tools::types::Float; use pointsource_algs::measures::RNDM; @@ -7,9 +10,12 @@ use pyo3::prelude::*; use std::marker::PhantomData; +/// A factor for generating a [`PythonPlotter`]. #[derive(Debug, Clone)] pub struct PythonPlotFactory<'py, T1, T2> { + /// Python side object implementing the plot factory pub(super) obj: Option<Bound<'py, PyAny>>, + /// Phantoms pub(super) _phantoms: PhantomData<(T1, T2)>, } @@ -24,8 +30,10 @@ } } +/// Python-side implementation of [`Plotter`] #[derive(Debug, Clone)] pub struct PythonPlotter<'py, T1, T2> { + /// Python side object implementing the [`Plotter`] pub(super) obj: Option<Bound<'py, PyAny>>, pub(super) _phantoms: PhantomData<(T1, T2)>, } @@ -60,12 +68,15 @@ } impl<'py, T1, T2> PythonPlotFactory<'py, T1, T2> { + /// Creates a dummy plotter factory that doe snothing pub fn dummy() -> Self { PythonPlotFactory { obj: None, _phantoms: PhantomData } } } impl<'py, T1: Clone, T2: Clone> PythonPlotFactory<'py, T1, T2> { + /// Produces a [`PythonPlotter`] out of the factory. + /// The `prefix` parameter indicates where to store the files. pub fn produce(&self, prefix: String) -> PythonPlotter<'py, T1, T2> { if let Some(ref obj) = self.obj { let produce = intern!(obj.py(), "produce");
--- a/src/quadratic_dataterm.py Thu Feb 26 09:32:12 2026 -0500 +++ b/src/quadratic_dataterm.py Wed Apr 22 22:32:00 2026 -0500 @@ -43,11 +43,12 @@ ma = self.opA.bound(xbound=xbound) mb = self.b.norm() mda = self.opA.diff_bound(xbound=xbound) - lda = self.opA.diff_adj_lipschitz_factor() + lda = self.opA.diff_chain_lipschitz_factor() return self.λ * (lda * (ma + mb) + mda * la) - # def diff_lipschitz_factor_pair(self, *args): - # return self.diff_lipschitz_factor() - - def diff_bound(self, xbound=None): - return self.λ * (self.opA.codomain_bound(xbound=xbound) + self.b_norm) + def diff_bound(self): + if hasattr(self.opA, "opnorm"): + opn = self.opA.opnorm() + return self.λ * opn * (opn * self.xbound + self.b_norm) + else: + raise Exception("Unimplemented")
