# HG changeset patch # User Tuomo Valkonen # Date 1791136946 18000 # Node ID a0d61ac3719e72278739016a8eac275ed2a18441 # Parent c4f92301734935cc28f12175b0a442890d04bd2c Various revisions: * Proper line integrals, * Boundary condition implementation fixes * Higher simulation resolution for experiments * Related parameter adjustments * Updates related to updated pointsource_algs diff -r c4f923017349 -r a0d61ac3719e Cargo.toml --- a/Cargo.toml Mon Jun 29 00:20:42 2026 +0300 +++ b/Cargo.toml Sun Oct 04 13:02:26 2026 -0500 @@ -1,6 +1,6 @@ [package] name = "pointsource_pde" -version = "0.1.1" +version = "0.2.1" edition = "2021" rust-version = "1.85" authors = ["Tuomo Valkonen "] @@ -12,17 +12,17 @@ categories = ["mathematics", "science", "computer-vision"] [dependencies.alg_tools] -version = "~0.4.1-dev" +version = "~0.5.4" path = "../alg_tools" default-features = false features = ["nightly", "pyo3"] [dependencies.pointsource_algs] -version = "~3.0.2-dev" +version = "~3.2.0" path = "../pointsource_algs" [dependencies.measures] -version = "~0.1.0" +version = "~0.3.0" path = "../measures" features = ["pyo3"] @@ -30,7 +30,7 @@ serde = { version = "~1.0", features = ["derive"] } cxx = "~1.0" pyo3 = { version = "0.27.1", features = ["anyhow"] } -anyhow = "~1.0.95" +anyhow = "~1.0" itertools = "~0.14.0" dolfinx-sys = { path = "dolfinx-sys" } nanobind-sys = { path = "nanobind-sys" } @@ -44,14 +44,18 @@ "release_max_level_warn", ] } colog = "1.3.0" +rand = "~0.9.5" [build-dependencies] build-print = "~0.1.1" cc = { version = "~1.2.22", features = ["parallel"] } cxx-build = "~1.0" #pyo3 = "~0.25.0" -regex = "~1.11.0" +regex = "~1.12.0" conda-build = { path = "conda-build" } [profile.release] debug = true + +[lints.cargo] +non_kebab_case_bins = "allow" diff -r c4f923017349 -r a0d61ac3719e README.md --- a/README.md Mon Jun 29 00:20:42 2026 +0300 +++ b/README.md Sun Oct 04 13:02:26 2026 -0500 @@ -1,5 +1,4 @@ - -# Nuerical implementation of “Leak localisation with a measure source convection–diffusion model” +# Numerical implementation of “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”_ @@ -13,7 +12,7 @@ (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.) -### Phase 1: Python and Fenics +### Phase 1: Python and Fenics #### Option 1.A that avoids Conda hell, but is more work (macOS and Linux) @@ -28,18 +27,18 @@ 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: -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 +| 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: @@ -65,7 +64,7 @@ For changes `.envrc` to take effect, you should use direnv allow - + ##### Phase 1.A.3: Fenicx-basix Install basix from according to @@ -111,16 +110,15 @@ 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 can _try to_ install Fenicsx in Conda according to instructions on the Fenics website. Additionally you need to install `scipy`: 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. +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. #### Option 1.C: Debian/Ubuntu @@ -135,37 +133,37 @@ 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 - + rustup toolchain install nightly -3. Download [alg_tools], [pointsource_algs], and [measures] and unpack them 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/ - [Homebrew]: https://brew.sh - [arXiv:2212.02991]: https://arxiv.org/abs/2212.02991 - [arXiv:2502.12417]: https://arxiv.org/abs/2502.12417 - [doi:10.46298/jnsao-2023-10433]: http://doi.org/10.46298/jnsao-2023-10433 +[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/ +[Homebrew]: https://brew.sh +[arXiv:2212.02991]: https://arxiv.org/abs/2212.02991 +[arXiv:2502.12417]: https://arxiv.org/abs/2502.12417 +[doi:10.46298/jnsao-2023-10433]: http://doi.org/10.46298/jnsao-2023-10433 ### 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: - 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. +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. +3. 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. +4. 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 @@ -187,7 +185,7 @@ ### Visualising the results -The results may be plotted with +The results may be plotted with python3 ./plot.py results/laser_and_mirrors_aux/radon_sliding_fb diff -r c4f923017349 -r a0d61ac3719e build.rs --- a/build.rs Mon Jun 29 00:20:42 2026 +0300 +++ b/build.rs Sun Oct 04 13:02:26 2026 -0500 @@ -2,7 +2,6 @@ Build script for `pointsource_pde` */ -#![feature(path_is_empty)] //use pyo3::{prepare_freethreaded_python, types::PyModule, PyErr, Python}; use conda_build::python_config; diff -r c4f923017349 -r a0d61ac3719e conda-build/Cargo.toml --- a/conda-build/Cargo.toml Mon Jun 29 00:20:42 2026 +0300 +++ b/conda-build/Cargo.toml Sun Oct 04 13:02:26 2026 -0500 @@ -10,5 +10,5 @@ crate-type = ["dylib", "rlib"] [dependencies] -anyhow = "~1.0.98" +anyhow = "~1.0" python-config-rs = "0.1.2" diff -r c4f923017349 -r a0d61ac3719e dolfinx-sys/Cargo.toml --- a/dolfinx-sys/Cargo.toml Mon Jun 29 00:20:42 2026 +0300 +++ b/dolfinx-sys/Cargo.toml Sun Oct 04 13:02:26 2026 -0500 @@ -7,7 +7,7 @@ links = "dolfinx" [build-dependencies] -anyhow = "~1.0.98" +anyhow = "~1.0" build-print = "~0.1.1" itertools = "~0.14.0" pkg-config = "~0.3" diff -r c4f923017349 -r a0d61ac3719e experiments/laser_and_mirrors.py --- a/experiments/laser_and_mirrors.py Mon Jun 29 00:20:42 2026 +0300 +++ b/experiments/laser_and_mirrors.py Sun Oct 04 13:02:26 2026 -0500 @@ -2,6 +2,8 @@ # Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity # import laser_and_mirrors_aux +import numpy as np +from laser_and_mirrors_aux import generic_setup, reg_params from measures import DiscreteMeasure_2_f64 from pointsource_pde import Problem, RegTerm from pointsource_pde.compose import ( @@ -16,21 +18,46 @@ algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides -# Setup routine -def setup(prefix): - dat, auxtrue, _μ_bound, μ_true, plot_factory = laser_and_mirrors_aux.generic_setup( - prefix +# Setup the problem +def setup( + prefix, + alpha=reg_params["alpha"], + c_noise=reg_params["c_noise"], + k_noise=reg_params["k_noise"], + noise=reg_params["noise"], + seed=reg_params["seed"], + sim_oversample=reg_params["sim_oversample"], + sim_oversample_time=reg_params["sim_oversample_time"], +): + rng = np.random.default_rng(seed=seed) + + dat, auxtrue, aux0, _μ_bound, μ_true, plot_factory, pde = generic_setup( + prefix, + rng=rng, + c_noise=c_noise, + k_noise=k_noise, + noise=noise, + sim_oversample=sim_oversample, + sim_oversample_time=sim_oversample_time, ) + reg = RegTerm.NonnegRadon(alpha) + μ0 = DiscreteMeasure_2_f64([]) - reg = RegTerm.NonnegRadon(1.5e-6) - dat_simple = ComposeFnWithOperator(dat, InjectSecond(auxtrue)) + dat_simple.curvature_bound_components = lambda: (None, None, None) + dat_simple.override_lipschitz = 1.0 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 + dat_simple, + reg, + None, + None, + μ0, + plot_factory=plot_factory, + domain=pde.domain_dims, ) diff -r c4f923017349 -r a0d61ac3719e experiments/laser_and_mirrors2.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/experiments/laser_and_mirrors2.py Sun Oct 04 13:02:26 2026 -0500 @@ -0,0 +1,67 @@ +# +# Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity +# +import laser_and_mirrors_aux +import numpy as np +from laser_and_mirrors_aux import generic_setup, reg_params +from measures import DiscreteMeasure_2_f64 +from pointsource_pde import Problem, RegTerm +from pointsource_pde.compose import ( + ComposeFnWithOperator, + InjectSecond, +) + +# Give name to the problem +name = "laser_and_mirrors2" + +# Override algorithm settings +algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides + + +# Setup the problem +def setup( + prefix, + alpha=reg_params["alpha"], + c_noise=reg_params["c_noise"], + k_noise=reg_params["k_noise"], + noise=reg_params["noise"], + seed=reg_params["seed"], + sim_oversample=reg_params["sim_oversample"], + sim_oversample_time=reg_params["sim_oversample_time"], +): + rng = np.random.default_rng(seed=seed) + + dat, auxtrue, aux0, _μ_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, + c_noise=c_noise, + k_noise=k_noise, + noise=noise, + sim_oversample=sim_oversample, + sim_oversample_time=sim_oversample_time, + ) + + reg = RegTerm.NonnegRadon(alpha) + + μ0 = DiscreteMeasure_2_f64([]) + + dat_simple = ComposeFnWithOperator(dat, InjectSecond(auxtrue)) + dat_simple.curvature_bound_components = lambda: (None, None, None) + dat_simple.override_lipschitz = 1.0 + + print("Initial data term value:", dat_simple.apply(μ0)) + print("Data term value at true μ:", dat_simple.apply(μ_true)) + + return Problem( + dat_simple, + reg, + None, + None, + μ0, + plot_factory=plot_factory, + domain=pde.domain_dims, + ) diff -r c4f923017349 -r a0d61ac3719e experiments/laser_and_mirrors_aux.py --- a/experiments/laser_and_mirrors_aux.py Mon Jun 29 00:20:42 2026 +0300 +++ b/experiments/laser_and_mirrors_aux.py Sun Oct 04 13:02:26 2026 -0500 @@ -10,6 +10,7 @@ from pointsource_pde import Problem, RegTerm from pointsource_pde.compose import ( ComposeFnWithOperator, + ScaleAuxStepHack, SumOfSeparableFunctions, ) from pointsource_pde.convection_diffusion import ( @@ -21,22 +22,25 @@ from pointsource_pde.laser_sampling import LaserSampling from pointsource_pde.quadratic_dataterm import QuadraticDataTerm +from dolfinx_access import eval_Function_f64 + # Give name to the problem name = "laser_and_mirrors_aux" # Fine-tune algorithms alg_finetune = { - "tau0": 0.99, - "tautheta": 100, - "gradv_lipest_mult": 1.3, - "transport_tolerance": 10, - "inner_method": "PP", - "inner_pdps_τσ0": (19.8, 0.05), - "inner_pp_τ": [1000.0, 1000.0], + "tau0": 0.02, + "sigmap0": 0.0001, + "tautheta": 0.07, + "gradv_lipest_mult": 1, + "transport_tolerance": 1e5, + "inner_method": "Exact", "merge": True, "fitness_merging": True, "merge_radius": 0.1, "merge_interp": False, + "scaling_heuristic": None, + "extra_weight_optimisation_steps": None, } algorithm_overrides = { @@ -46,6 +50,22 @@ "RadonSlidingPDPS": alg_finetune, } +reg_params = { + "alpha": 2e-3, + "alpha_k": 60, + "alpha_c": 0.06, + "c_noise": 0.02, + "k_noise": 0.02, + "seed": 2393, + "true_init": False, + "true_target": False, + "noise": 0.01, + "sim_oversample": 2, + "sim_oversample_time": 2, + "scale_cstep": 300, + "do_non_sliding": True, +} + # Setup shared between experiment with known/unknown conductivity and diffusivity def generic_setup( @@ -58,16 +78,39 @@ k=0.01, r=0.5, θ=30 * np.pi / 180, + c_noise=reg_params["c_noise"], + k_noise=reg_params["k_noise"], + noise=reg_params["noise"], + sim_oversample=reg_params["sim_oversample"], + sim_oversample_time=reg_params["sim_oversample_time"], + scale_cstep=reg_params["scale_cstep"], ): + u0 = lambda x: 0.0 * x[1] + g = lambda x: 0.0 * x[1] + domain_dims = [[0, 0.5], [0, 0.5]] + num_steps = 50 + res = 32 + # Create convection-diffusion PDE pde = ConvectionDiffusion( - u0=lambda x: 0.0 * x[1], - 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), + u0=u0, + g=g, + domain_dims=domain_dims, + num_steps=num_steps, + nx=res, + ny=res, + ) + + # Create convection-diffusion PDE for simulation, to avoid inverse crimes + pde_sim = ConvectionDiffusion( + u0=u0, + g=g, + domain_dims=domain_dims, + # num_steps=(sim_oversample-1) * (num_steps + 1) + num_steps, + num_steps=sim_oversample_time * (num_steps - 1) + 1, + nx=res * sim_oversample, + ny=res * sim_oversample, ) # Set up lasers and mirrors @@ -79,10 +122,10 @@ ] mirrors_per_side = 10 x0, x1, y0, y1 = ( - pde.domain_size[0], - pde.domain_size[1], - pde.domain_size[2], - pde.domain_size[3], + pde.domain_dims[0][0], + pde.domain_dims[0][1], + pde.domain_dims[1][0], + pde.domain_dims[1][1], ) beams = [ (src, np.array(dst)) @@ -97,6 +140,7 @@ ] ] L_i = LaserSampling(pde.V, pde.domain, pde.nx, pde.ny, beams) + L_i_sim = LaserSampling(pde_sim.V, pde_sim.domain, pde_sim.nx, pde_sim.ny, beams) v0 = 1.0 @@ -106,6 +150,9 @@ c2 = r * np.sin(θ) # k = 0.001 # high wind print("True aux variables: k = %f, c1 = %f, c2 = %f" % (k, c1, c2)) + k_sim = k + c1_sim = c1 + c2_sim = c2 else: # Convectivity and diffusivity fields def velocity_x(x): @@ -124,13 +171,33 @@ k0 = 1.0 k.interpolate(lambda x: k0 + 0.2 * np.sin(5 * x[0]) * np.cos(5 * x[1])) # k is + c1_sim = fem.Function(pde_sim.V) + c1_sim.interpolate(velocity_x) + c2_sim = fem.Function(pde_sim.V) + c2_sim.interpolate(velocity_y) + k_sim = fem.Function(pde_sim.V) + k0 = 1.0 + k_sim.interpolate( + lambda x: k0 + 0.2 * np.sin(5 * x[0]) * np.cos(5 * x[1]) + ) # k is + auxtrue = (k, (c1, c2)) # Simulate measurements + u_n_list_sim = pde_sim.apply((μ_true, (k_sim, (c1_sim, c2_sim)))) u_n_list = pde.apply((μ_true, (k, (c1, c2)))) - b_true = [L_i.apply(u) for u in u_n_list] + + b_true = [L_i_sim.apply(u) for u in u_n_list_sim[0::sim_oversample_time]] + + # Add noise + level_k = abs(k) + level_c = r + k_noisy = max(0.001, relnoise(k, k_noise, rng, level=level_k)) + c1_noisy, c2_noisy = relnoise((c1, c2), c_noise, rng, level=level_c) + aux_noisy = (k_noisy, (c1_noisy, c2_noisy)) b = [ - b_i + L_i.noise(0.01 * norm(b_i) / np.sqrt(np.size(b_i)), rng) for b_i in b_true + b_i + L_i_sim.noise(noise * norm(b_i) / np.sqrt(np.size(b_i)), rng) + for b_i in b_true ] # Estimate bounds on domains and values @@ -154,62 +221,118 @@ ), pde, ) + if scale_cstep is not None: + dat = ComposeFnWithOperator(dat, ScaleAuxStepHack((1.0, scale_cstep))) + dat.curvature_bound_components = lambda: (None, None, None) # Save true and initial data plot_factory = PlotFactory(pde) plotter = plot_factory.produce(prefix) ω, _ = dat.diff((DiscreteMeasure_2_f64([]), auxtrue)) - ω_true, _ = dat.diff((μ_true, auxtrue)) + ω_true, (dk, (dc1, dc2)) = dat.diff((μ_true, auxtrue)) np.savez_compressed( "%s/omega_0.npz" % prefix, 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, + k_true=k, + c1_true=c1, + c2_true=c2, + k_noisy=k_noisy, + c1_noisy=c1_noisy, + c2_noisy=c2_noisy, true_u_n_list_array=[u_n_list[i].x.array for i in plotter.plot_frames], ) - return dat, auxtrue, μ_bound, μ_true, plot_factory, pde + delta = 0.001 * k + dk_est = ( + dat.apply((μ_true, (k + delta, (c1, c2)))) - dat.apply((μ_true, auxtrue)) + ) / delta + print("k-diff-test:", dk, ", ", dk_est) + delta = 0.001 * c1 + dc1_est = ( + dat.apply((μ_true, (k, (c1 + delta, c2)))) - dat.apply((μ_true, auxtrue)) + ) / delta + print("c1-diff-test:", dc1, ", ", dc1_est) + delta = 0.001 + μ_perturb = DiscreteMeasure_2_f64( + [([0.1, 0.3], 0.08 + delta), ([0.4, 0.25], 0.05), ([0.25, 0.13], 0.06)] + ) + dcμ_est = (dat.apply((μ_perturb, auxtrue)) - dat.apply((μ_true, auxtrue))) / delta + print( + "μ-diff-test:", + eval_Function_f64(ω_true._cpp_object, [0.1, 0.3, 0.0]), + ", ", + dcμ_est, + ) + + return dat, auxtrue, aux_noisy, μ_bound, μ_true, plot_factory, pde def relnoise(v, σ, rng, level=None): - if level is None: - level = abs(v) - return float(rng.normal(v, σ * level)) + if isinstance(v, float): + return float(rng.normal(v, σ * (abs(v) if level is None else level))) + else: + return map( + float, rng.normal(v, σ * (np.linalg.norm(v) if level is None else level)) + ) # 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) +def setup( + prefix, + alpha=reg_params["alpha"], + alpha_k=reg_params["alpha_k"], + alpha_c=reg_params["alpha_c"], + c_noise=reg_params["c_noise"], + k_noise=reg_params["k_noise"], + noise=reg_params["noise"], + seed=reg_params["seed"], + sim_oversample=reg_params["sim_oversample"], + sim_oversample_time=reg_params["sim_oversample_time"], + true_init=False, + true_target=False, + scale_cstep=reg_params["scale_cstep"], +): + rng = np.random.default_rng(seed=seed) - # 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) + dat, auxtrue, aux0, _μ_bound, μ_true, plot_factory, pde = generic_setup( + prefix, + rng=rng, + c_noise=c_noise, + k_noise=k_noise, + noise=noise, + sim_oversample=sim_oversample, + sim_oversample_time=sim_oversample_time, + scale_cstep=scale_cstep, + ) - reg = RegTerm.NonnegRadon(1.5e-6) + # Override Lipschitz parameter to use direct step length choice + dat.override_lipschitz = 1.0 + dat.override_lipschitz_pair = (0.0, 0.0) + + reg = RegTerm.NonnegRadon(alpha) μ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.005, -1.0), + (1.0, 1.0), + (alpha_k, alpha_c), + auxtrue if true_target else aux0, + (1, scale_cstep) if scale_cstep is not None else None, ) - 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) + print("Initial value:", dat.apply((μ0, aux0)) + aux.apply(aux0)) + print("Value at true μ:", dat.apply((μ_true, auxtrue)) + aux.apply(auxtrue)) - return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory) + return Problem( + dat, + reg, + aux, + auxtrue if true_init else aux0, + μ0, + plot_factory=plot_factory, + domain=pde.domain_dims, + ) diff -r c4f923017349 -r a0d61ac3719e experiments/laser_and_mirrors_aux2.py --- a/experiments/laser_and_mirrors_aux2.py Mon Jun 29 00:20:42 2026 +0300 +++ b/experiments/laser_and_mirrors_aux2.py Sun Oct 04 13:02:26 2026 -0500 @@ -4,7 +4,7 @@ # import laser_and_mirrors_aux import numpy as np -from laser_and_mirrors_aux import generic_setup, relnoise +from laser_and_mirrors_aux import generic_setup, reg_params from measures import DiscreteMeasure_2_f64 from pointsource_pde import Problem, RegTerm from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation @@ -17,41 +17,66 @@ # Setup the problem -def setup(prefix): - rng = np.random.default_rng(seed=31337) +def setup( + prefix, + alpha=reg_params["alpha"], + alpha_k=reg_params["alpha_k"], + alpha_c=reg_params["alpha_c"], + c_noise=reg_params["c_noise"], + k_noise=reg_params["k_noise"], + noise=reg_params["noise"], + seed=reg_params["seed"], + sim_oversample=reg_params["sim_oversample"], + sim_oversample_time=reg_params["sim_oversample_time"], + true_init=False, + true_target=False, + scale_cstep=reg_params["scale_cstep"], +): + rng = np.random.default_rng(seed=seed) - dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup( + dat, auxtrue, aux0, _μ_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, + c_noise=c_noise, + k_noise=k_noise, + noise=noise, + sim_oversample=sim_oversample, + sim_oversample_time=sim_oversample_time, + scale_cstep=scale_cstep, ) # Override Lipschitz parameter - pde.override_lipschitz = (4.0,) - pde.override_lipschitz_pair = (4.0, 4.0) + dat.override_lipschitz = 1.0 + dat.override_lipschitz_pair = (0.0, 0.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) + reg = RegTerm.NonnegRadon(alpha) μ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), + (alpha_k, alpha_c), + auxtrue if true_target else aux0, + (1, scale_cstep) if scale_cstep is not None else None, ) - 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) + print("Initial value:", dat.apply((μ0, aux0)) + aux.apply(aux0)) + print("Value at true μ:", dat.apply((μ_true, auxtrue)) + aux.apply(auxtrue)) - return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory) + return Problem( + dat, + reg, + aux, + auxtrue if true_init else aux0, + μ0, + plot_factory=plot_factory, + domain=pde.domain_dims, + ) diff -r c4f923017349 -r a0d61ac3719e experiments/tests/test_minimisation.py --- a/experiments/tests/test_minimisation.py Mon Jun 29 00:20:42 2026 +0300 +++ b/experiments/tests/test_minimisation.py Sun Oct 04 13:02:26 2026 -0500 @@ -1,9 +1,23 @@ +import numpy as np +from dolfinx import fem, mesh +from mpi4py import MPI + from dolfinx_access import eval_Function_f64, min_Function_f64_p2 # A simple test that the finite-dimensional subproblem is solved correctly -def test_minimisation(pde): - test = fem.Function(pde.V) +def test_minimisation(domain_dims=[[-5, 5], [-5, 5]], degree=2, nx=256, ny=256, n=512): + domain = mesh.create_rectangle( + MPI.COMM_SELF, + [ + np.array([domain_dims[0][0], domain_dims[1][0]]), + np.array([domain_dims[0][1], domain_dims[1][1]]), + ], + [nx, ny], + mesh.CellType.triangle, + ) + V = fem.functionspace(domain, ("Lagrange", degree)) + test = fem.Function(V) # Beale function def beale(x): @@ -13,22 +27,69 @@ + (2.625 - x[0] + x[0] * (x[1] ** 3)) ** 2 ) - test.interpolate(beale) - m = min_Function_f64_p2(test._cpp_object) - print("Beale minimizer:", m.x, m.v) - - def quadratic_testf(x): + def quadratic(x): return x[0] ** 2 - 5 * x[1] ** 2 + x[0] * x[1] - x[0] + x[1] - testf = quadratic_testf - test.interpolate(testf) - n = 256 - dist = 0.0 - for x in np.linspace(pde.domain_size[0], pde.domain_size[1], n): - for y in np.linspace(pde.domain_size[2], pde.domain_size[3], n): - p = [x, y, 0.0] - v = eval_Function_f64(test._cpp_object, p) - v0 = testf(p) - d = np.linalg.norm(v - v0) - dist = max(d, dist) - print("Quadratic function evaluation max discrepancy:", dist) + def sphere(x): + return x[0] ** 2 + x[1] ** 2 + + def booth(x): + return (x[0] + 2 * x[1] - 7) ** 2 + (2 * x[0] + x[1] - 5) ** 2 + + def himmelblau(x): + return (x[0] ** 2 + x[1] - 11) ** 2 + (x[0] + x[1] ** 2 - 7) ** 2 + + for testf, name in [ + (beale, "Beale"), + (quadratic, "Quadratic"), + (sphere, "Sphere"), + (booth, "Booth"), + (himmelblau, "Himmelblau"), + ]: + test.interpolate(testf) + m = min_Function_f64_p2(test._cpp_object) + print( + name, + "minimizer:", + m.v, + "(verification ", + eval_Function_f64(test._cpp_object, [m.x[0], m.x[1], 0]), + " original ", + testf(m.x), + ") @ ", + m.x, + ) + + dist = 0.0 + m = float("inf") + m0 = float("inf") + mp = None + for x in np.linspace(domain_dims[0][0], domain_dims[0][1], n): + for y in np.linspace(domain_dims[1][0], domain_dims[1][1], n): + p = [x, y, 0.0] + v = eval_Function_f64(test._cpp_object, p) + v0 = testf(p) + if v0 < m0: + m = v + m0 = v0 + mp = p + d = np.linalg.norm(v - v0) + dist = max(d, dist) + print(name, "evaluation max discrepancy:", dist) + print( + name, + "rough scan minimiser ", + m, + "(model", + m0, + ") @ [", + float(mp[0]), + ", ", + float(mp[1]), + "]", + ) + print() + + +test_minimisation() +exit(0) diff -r c4f923017349 -r a0d61ac3719e extract_aux.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/extract_aux.py Sun Oct 04 13:02:26 2026 -0500 @@ -0,0 +1,95 @@ +#!.venv/bin/python3 +import argparse +import csv +from pathlib import Path + +import numpy as np + +from plot import find_files + + +def extract_aux(prefix): + + iter_files = find_files(prefix) + + init_filename = "%s/init_aux.csv" % Path(prefix).parent + filename = "%s/res_aux.csv" % prefix + + # Do true and noisy values + + initdata = [ + [ + "k_true", + "c1_true", + "c2_true", + "k_noisy", + "c1_noisy", + "c2_noisy", + ] + ] + + data0 = np.load("%s/omega_0.npz" % Path(prefix).parent) + k_true = data0["k_true"] + c1_true = data0["c1_true"] + c2_true = data0["c2_true"] + initdata.append( + [ + k_true, + c1_true, + c2_true, + data0["k_noisy"], + data0["c1_noisy"], + data0["c2_noisy"], + ] + ) + with open(init_filename, mode="w", newline="") as file: + # Create a csv.writer object + writer = csv.writer(file, delimiter="\t") + # Write data to the CSV file + writer.writerows(initdata) + + # Do reconstructed values + + auxdata = [ + [ + "iter", + "k", + "c1", + "c2", + "k_dist", + "c_dist", + ] + ] + + for iter, file in iter_files: + d = np.load(file) + k = d["k"] + c1 = d["c1"] + c2 = d["c2"] + k_dist = np.abs(k - k_true) + c_dist = np.sqrt((c1 - c1_true) ** 2 + (c2 - c2_true) ** 2) + auxdata.append([iter, k, c1, c2, k_dist, c_dist]) + + # Open the file in write mode + with open(filename, mode="w", newline="") as file: + # Create a csv.writer object + writer = csv.writer(file, delimiter="\t") + # Write data to the CSV file + writer.writerows(auxdata) + + +def main(): + parser = argparse.ArgumentParser( + prog="extract_aux.py", + description="Extract auxiliary parameters from pointsource_pde", + ) + + parser.add_argument("filename") + + args = parser.parse_args() + + extract_aux(args.filename) + + +if __name__ == "__main__": + main() diff -r c4f923017349 -r a0d61ac3719e include/dolfinx_access/minmax_p2.h --- a/include/dolfinx_access/minmax_p2.h Mon Jun 29 00:20:42 2026 +0300 +++ b/include/dolfinx_access/minmax_p2.h Sun Oct 04 13:02:26 2026 -0500 @@ -8,6 +8,9 @@ extern CoordValuePair max_Function_f64_p2(Function_f64 const* f); extern CoordValuePair minmax_Function_f64_p2(Function_f64 const* f, bool max); + void cell_dofs_Function_f64_p2(Function const* f, int32_t cell, + std::array& x, std::array& vals); + } // namespace dolfinx_access // extern "C" PyObject* PyInit_dolfinx_access(); diff -r c4f923017349 -r a0d61ac3719e misc/copy_results.sh --- a/misc/copy_results.sh Mon Jun 29 00:20:42 2026 +0300 +++ b/misc/copy_results.sh Sun Oct 04 13:02:26 2026 -0500 @@ -1,13 +1,15 @@ #!/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" +iter=10000 +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 +src1=results_alt14a/laser_and_mirrors_aux +src2=results_alt14a/laser_and_mirrors_aux2 ./plot.py ${src1}/radon_fb --simple --iter ${iter} --save-only +./extract_aux.py ${src1}/radon_fb ./plot.py ${src1}/radon_sliding_fb --simple --iter ${iter} --save-only +./extract_aux.py ${src1}/radon_sliding_fb 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 @@ -16,7 +18,9 @@ done ./plot.py ${src2}/radon_fb --simple --iter ${iter} --save-only +./extract_aux.py ${src2}/radon_sliding_fb ./plot.py ${src2}/radon_sliding_fb --simple --iter ${iter} --save-only +./extract_aux.py ${src2}/radon_fb 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 diff -r c4f923017349 -r a0d61ac3719e nanobind-sys/build.rs --- a/nanobind-sys/build.rs Mon Jun 29 00:20:42 2026 +0300 +++ b/nanobind-sys/build.rs Sun Oct 04 13:02:26 2026 -0500 @@ -28,7 +28,7 @@ }) .map_err(|e| anyhow::Error::from(std::io::Error::from(e)))?; - let prefix = prefix_override.ok_or_else(|| pyc.prefix_path()).unwrap(); + let prefix = prefix_override.unwrap_or_else(|| pyc.prefix_path().unwrap()); let version = pyc.semantic_version().unwrap(); let nanobind_root = ({ diff -r c4f923017349 -r a0d61ac3719e plot.py --- a/plot.py Mon Jun 29 00:20:42 2026 +0300 +++ b/plot.py Sun Oct 04 13:02:26 2026 -0500 @@ -14,6 +14,7 @@ import matplotlib.tri as tri import numpy as np from dolfinx import fem + from src.convection_diffusion import ConvectionDiffusion @@ -49,9 +50,10 @@ u0=lambda x: 0.0 * x[1], g=lambda x: 0.0 * x[1], # ((x[1] + 2) / 4) ** 2, # nonzero bcs # w0=lambda x: np.cos(5 * x[1]), - domain_size=[0, 0.5, 0, 0.5], - # nx=128, - # ny=128, + domain_dims=[[0, 0.5], [0, 0.5]], + nx=32, + ny=32, + num_steps=50, ) # Get coordinates @@ -135,9 +137,14 @@ true_μ_alpha = list(map(lambda x: x[2], true_μ)) if len(μ_alpha) == 0: - alpha_mi, alpha_ma, alpha_me = 0, 0, 0 + alpha_mi, alpha_ma, alpha_me, alpha_sum = 0, 0, 0, 0 else: - alpha_mi, alpha_ma, alpha_me = min(μ_alpha), max(μ_alpha), median(μ_alpha) + alpha_mi, alpha_ma, alpha_me, alpha_sum = ( + min(μ_alpha), + max(μ_alpha), + median(μ_alpha), + sum(μ_alpha), + ) true_alpha_mi, true_alpha_ma = min(true_μ_alpha), max(true_μ_alpha) alpha_base = min(true_alpha_mi, alpha_mi) @@ -314,13 +321,13 @@ ax.set_aspect("equal") plt.suptitle( - "Convection-Diffusion, iteration %d; len(μ) = %d; μ_min = %f, μ_max = %f, μ_median = %f;" + "Convection-Diffusion, iteration %d; len(μ) = %d; μ_min = %f, μ_max = %f, μ_sum = %f;" % ( iter, len(μ), alpha_mi, alpha_ma, - alpha_me, + alpha_sum, ), fontsize=14, ) @@ -441,17 +448,22 @@ # print("Saved: solution_evolution.png + solution_time.png") -parser = argparse.ArgumentParser( - prog="plot.py", - description="Plots results of pointsource_pde", -) +def main(): + 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") + 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() + args = parser.parse_args() + + plot(args.filename, args.simple, args.shaded, args.iter, args.save_only) -plot(args.filename, args.simple, args.shaded, args.iter, args.save_only) + +if __name__ == "__main__": + main() diff -r c4f923017349 -r a0d61ac3719e src/compose.py --- a/src/compose.py Mon Jun 29 00:20:42 2026 +0300 +++ b/src/compose.py Sun Oct 04 13:02:26 2026 -0500 @@ -43,6 +43,8 @@ def __init__(self, f, op): self.f = f self.op = op + self.override_lipschitz = None + self.override_lipschitz_pair = None def apply(self, *args): return self.f.apply(self.op.apply(*args)) @@ -93,7 +95,9 @@ 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"): + if self.override_lipschitz is not None: + return self.override_lipschitz + elif hasattr(self.op, "opnorm"): return self.f.diff_lipschitz_factor() * self.op.opnorm() ** 2 else: mdf = self.f.diff_bound() @@ -113,6 +117,9 @@ there is no special handling of linear opeartors. """ + if self.override_lipschitz_pair is not None: + return self.override_lipschitz_pair + mdf = self.f.diff_bound() lda1, lda2 = self.op.diff_chain_lipschitz_factor_pair() @@ -151,3 +158,31 @@ # raise Exception("Linear operators have unbounded range") # else: # return xbound + + +class ScaleAuxStepHack: + """ + This hack only scale the differential, to get different step length parameters wrt. k and c, + without rescaling the variables themselves. Essentially it amounts to computing the gradient wrt. + different inner products. + """ + + def __init__(self, scale): + self.scale = scale + + def apply(self, x): + # μ, (k, (c1, c2)) = x + # sck, scc = self.scale + # return μ, (k * sck, (c1 * scc, c2 * scc)) + return x + + def diff_adjdir(self, j, _x, apply_result=None): + dμ, (dk, (dc1, dc2)) = j + sck, scc = self.scale + return dμ, (dk * sck, (dc1 * scc, dc2 * scc)) + + def opnorm(self, *args): + return 1 + + def diff_bound(self): + return 1 diff -r c4f923017349 -r a0d61ac3719e src/convection_diffusion.py --- a/src/convection_diffusion.py Mon Jun 29 00:20:42 2026 +0300 +++ b/src/convection_diffusion.py Sun Oct 04 13:02:26 2026 -0500 @@ -64,7 +64,7 @@ u0, g, # w0, - domain_size=[-2, 2, -2, 2], + domain_dims=[[-2, 2], [-2, 2]], nx=32, ny=32, degree=2, @@ -89,11 +89,11 @@ self.C_emb = C_emb # embedding constant self.C_reg = C_reg if Delta_t is None: - self.Delta_t = T / num_steps # Matches your N=50 intent + self.Delta_t = T / num_steps else: self.Delta_t = Delta_t - self.domain_size = domain_size + self.domain_dims = domain_dims self.nx = nx self.ny = ny self.degree = degree @@ -112,19 +112,18 @@ domain = mesh.create_rectangle( MPI.COMM_SELF, [ - np.array([domain_size[0], domain_size[2]]), - np.array([domain_size[1], domain_size[3]]), + np.array([domain_dims[0][0], domain_dims[1][0]]), + np.array([domain_dims[0][1], domain_dims[1][1]]), ], [nx, ny], mesh.CellType.triangle, ) + # Compute the domain volume for the adjoint Lipschitz factor self.domain = domain - # Compute the domain volume for the adjoint Lipschitz factor - self.domain_size = domain_size - # Domain diameter from your domain_size + # Domain diameter from your domain self.diam = np.sqrt( - (domain_size[1] - domain_size[0]) ** 2 - + (domain_size[3] - domain_size[2]) ** 2 + (domain_dims[0][0] - domain_dims[0][1]) ** 2 + + (domain_dims[1][0] - domain_dims[1][1]) ** 2 ) V = fem.functionspace(domain, ("Lagrange", degree)) @@ -140,26 +139,36 @@ self.bcs = [] # Initialize BCS as LIST self.adjoint_bcs = [] # Initialize ADJOINT BCS as LIST (w = 0) - # Identify boundary facets for Dirichlet BC location (e.g. left/right boundaries at x=domain_size[0] or x=domain_size[1]) + # Identify boundary facets for Dirichlet BC location (e.g. left/right boundaries fdim = domain.topology.dim - 1 # Γ₁ facets: TOP/BOTTOM boundaries (y = ymin, ymax) - STATE Dirichlet u = g gamma1_facets = mesh.locate_entities_boundary( domain, fdim, lambda x: ( - np.isclose(x[1], domain_size[2]) # x[1] = y = ymin - | np.isclose(x[1], domain_size[3]) + np.isclose(x[1], domain_dims[1][0]) # x[1] = y = ymin + | np.isclose(x[1], domain_dims[1][1]) ), # x[1] = y = ymax ) - # Γ₂ facets: LEFT/RIGHT boundaries (x = xmin, xmax) - ADJOINT Dirichlet w = 0 + # Γ₂ facets: LEFT/RIGHT boundaries (x = xmin, xmax) gamma2_facets = mesh.locate_entities_boundary( domain, fdim, lambda x: ( - np.isclose(x[0], domain_size[0]) # x[0] = x = xmin - | np.isclose(x[0], domain_size[1]) + np.isclose(x[0], domain_dims[0][0]) # x[0] = x = xmin + | np.isclose(x[0], domain_dims[0][1]) ), # x[0] = x = xmax ) + # tag Γ2 facets (left/right) once in __init__ + facet_idx = np.concatenate([gamma1_facets, gamma2_facets]) + facet_vals = np.concatenate( + [ + np.full(len(gamma1_facets), 1, np.int32), + np.full(len(gamma2_facets), 2, np.int32), + ] + ) + order = np.argsort(facet_idx) + facet_tags = mesh.meshtags(domain, fdim, facet_idx[order], facet_vals[order]) # STATE BC: u = g on Γ₁ (top/bottom) dofs_gamma1 = fem.locate_dofs_topological(V, fdim, gamma1_facets) # Create BC and ADD TO LIST @@ -168,12 +177,12 @@ # ADJOINT BC: w = 0 on Γ₁ (same facets, homogenized) zero = fem.Constant(domain, PETSc.ScalarType(0.0)) - dofs_gamma2 = fem.locate_dofs_topological(V, fdim, gamma2_facets) - bc_adjoint = fem.dirichletbc(zero, dofs_gamma2, V) + # dofs_gamma2 = fem.locate_dofs_topological(V, fdim, gamma2_facets) + bc_adjoint = fem.dirichletbc(zero, dofs_gamma1, V) self.adjoint_bcs.append(bc_adjoint) self.fdim = fdim - self.boundary_facets = np.unique(np.concatenate((gamma1_facets, gamma2_facets))) + # self.boundary_facets = np.unique(np.concatenate((gamma1_facets, gamma2_facets))) ux = fem.Function(V) wpx = fem.Function(V) @@ -211,21 +220,21 @@ j_u = fem.Function(V) # New copy! j_u.name = "j_u" self.j_u = j_u + ds = ufl.Measure("ds", domain=domain, subdomain_data=facet_tags) + n = ufl.FacetNormal(domain) a2 = ( w * v * ufl.dx - + dt - * ( - ufl.dot(self.k * ufl.grad(w), ufl.grad(v)) - - (self.c0 * w.dx(0) * v + self.c1 * w.dx(1) * v) - ) - * ufl.dx + + dt * ufl.dot(self.k * ufl.grad(w), ufl.grad(v)) * ufl.dx + - dt * (self.c0 * w.dx(0) + self.c1 * w.dx(1)) * v * ufl.dx + + dt * (self.c0 * n[0] + self.c1 * n[1]) * w * v * ds(2) ) self.bilinear_form_w = fem.form(a2) # Step n=N-1: L2 uses w_n = w^N=0 → solves w^{N-1}. (∂t)_primal^* # Step n=N-2: L2 uses w_n = w^{N-1} → solves w^{N-2} # Replace your L2 with 2 simple forms (same notation) - L2 = w_n * v * ufl.dx + dt * j_u * v * ufl.dx # kown w_n (previous step) + # No dt scaling in j_u, as that's already included in it. + L2 = w_n * v * ufl.dx + j_u * v * ufl.dx # kown w_n (previous step) self.linear_form_w = fem.form(L2) # Solving forward PDE @@ -290,7 +299,7 @@ b_ps.setValuesLocal( cell_dofs, - alpha * phi, + alpha * phi * dt, addv=PETSc.InsertMode.ADD_VALUES, ) @@ -351,8 +360,6 @@ assert len(j) == num_steps, f"MISMATCH: j({len(j)}) != u_n_list({num_steps})" # Terminal condition: w(T) = 0 w_n = self.w_n - w_n.x.array[:] = 0.0 - w_n.x.scatter_forward() self.c0.value = c[0] self.c1.value = c[1] @@ -375,11 +382,12 @@ t = self.t0 + num_steps * dt # Start at time T # N+1 adjoint values to match N PDE steps - w_n_list = [None] * (num_steps + 1) + w_n_list = [None] * num_steps whT = fem.Function(V) whT.x.array[:] = 0.0 # w^N = 0 whT.x.scatter_forward() - w_n_list[num_steps] = whT # Store at FINAL index + w_n.x.array[:] = whT.x.array + w_n.x.scatter_forward() # backward steps**: compute w^{N-1}, ..., w^0 @@ -436,15 +444,7 @@ # 1. expr1: dual sensitivity w.r.t. μ = ∫₀^T w_p dt (always Function) wp_bar = fem.Function(V) - wp_bar.x.array[:] = ( - 0.5 - * dt - * ( - w_n_list[0].x.array - + w_n_list[-1].x.array - + 2 * np.sum([w.x.array for w in w_n_list[1:-1]], axis=0) - ) - ) + wp_bar.x.array[:] = dt * sum(w.x.array for w in w_n_list) ux = self.ux wpx = self.wpx @@ -462,8 +462,7 @@ wpx.x.array[:] = wp.x.array[:] wpx.x.scatter_forward() v2 = fem.assemble_scalar(form) - weight = 0.5 * dt if (n == 0 or n == len(u_n_list) - 1) else dt - expr2 -= weight * v2 + expr2 -= dt * v2 else: raise Exception("Unimplemented - out of date") @@ -475,8 +474,7 @@ for n, (u, wp) in enumerate(zip(u_n_list, w_n_list)): form = ufl.inner(ufl.grad(u), ufl.grad(wp)) * ufl.dx v2 = fem.assemble_scalar(fem.form(form)) - weight = 0.5 * dt if (n == 0 or n == len(u_n_list) - 1) else dt - expr2a -= weight * v2 + expr2a -= dt * v2 # 3. expr3: dual sensitivity w.r.t. c = (c1, c2) if is_scalar_c: @@ -492,8 +490,7 @@ wpx.x.array[:] = wp.x.array[:] wpx.x.scatter_forward() v3 = fem.assemble_scalar(forms[i]) - weight = 0.5 * dt if (n == 0 or n == len(u_n_list) - 1) else dt - expr3i -= weight * v3 + expr3i -= dt * v3 return expr3i expr3 = tuple(val(i) for i in range(dim)) @@ -509,8 +506,7 @@ for n, (u, wp) in enumerate(zip(u_n_list, w_n_list)): form = wp * ufl.grad(u)[i] * ufl.dx v3 = fem.assemble_scalar(fem.form(form)) - weight = 0.5 * dt if (n == 0 or n == len(u_n_list) - 1) else dt - expr3ia -= weight * v3 + expr3ia -= dt * v3 return expr3i expr3 = tuple(val(i) for i in range(dim)) @@ -781,8 +777,9 @@ def _own(self, x): return x if isinstance(x, float) else own(x) - def __init__(self, α, base=None): + def __init__(self, α, base=None, τscale=None): self.α = α + self.τscale = τscale if base is not None: (k, (c1, c2)) = base self.base = (self._own(k), (self._own(c1), self._own(c2))) @@ -846,6 +843,11 @@ else: αk, αc = self.α + if self.τscale is not None: + kscale, cscale = self.τscale + αk = αk * kscale + αc = αc * cscale + return ( self._prox1(τ * αk, k, kb), (self._prox1(τ * αc, c1, c1b), self._prox1(τ * αc, c2, c2b)), @@ -995,9 +997,9 @@ class BoxedQuadraticRegularisation: - def __init__(self, ω0, ω1, α, base=None): + def __init__(self, ω0, ω1, α, base=None, τscale=None): self.box = BoxConstraints(ω0, ω1) - self.quadratic = QuadraticRegularisation(α, base) + self.quadratic = QuadraticRegularisation(α, base, τscale) def apply(self, x): return self.box.apply(x) + self.quadratic.apply(x) diff -r c4f923017349 -r a0d61ac3719e src/dolfinx_access.rs --- a/src/dolfinx_access.rs Mon Jun 29 00:20:42 2026 +0300 +++ b/src/dolfinx_access.rs Sun Oct 04 13:02:26 2026 -0500 @@ -10,7 +10,6 @@ use anyhow::bail; use pyo3::ffi::PyImport_AppendInittab; use std::ffi::CString; -use std::mem::MaybeUninit; /// Import helper required for the linking to the sys crates to actually happen. #[allow(unused_imports)] @@ -98,6 +97,8 @@ max: bool, ) -> Result; + pub unsafe fn cell_dofs_Function_f64_p2(f: *const Function_f64, cell:i32, coords: &mut [f64; 18], vals: &mut [f64; 6]); + /// Checks if a Dolfinx `Function` is supported (P2 elements on triangular mesh) pub unsafe fn check_Function_f64(o: *const PyObject) -> bool; @@ -147,10 +148,9 @@ /// `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], - cell: i32, + unsafe fn minmax_dolfinx_p2( + simplex_coords: [f64; 18], + vals: [f64; 6], max: bool, ) -> CoordValuePair; } @@ -158,59 +158,25 @@ use ffi::CoordValuePair; -/// Convert an array of 2D oordinates obtained from Fenics, which hard-codes weird 3D coordinates. -#[inline] -fn coords_from_fenics(&[x1, x2, _, y1, y2, _, z1, z2, _]: &[f64; 9]) -> [Loc<2, f64>; 3] { - [Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])] -} -// fn coords_from_fenics(coords: &[[f64; 3]; N]) -> [Loc; N] { -// coords.map(|[x1, x2, _]| Loc([x1, x2])) -// } - -/// Convert an array of 2D oordinates to Fenics, which hard-codes weird 3D coordinates. -#[inline] -fn coords_to_fenics([Loc([x1, x2]), Loc([y1, y2]), Loc([z1, z2])]: [Loc<2, f64>; 3]) -> [f64; 9] { - [x1, x2, 0.0, y1, y2, 0.0, z1, z2, 0.0] -} -// fn coords_to_fenics(coords: [Loc; N]) -> [[f64; 3]; N] { -// coords.map(|Loc([x1, x2])| [x1, x2, 0.0]) -// } - /// Given a dolphinx `Function` `f`, as well as 6 coordinates on a simplex in its mesh /// (the corners and midpoints), construct a new internal P2 model on `cell`. /// If `max` is true, the model is negated compared to `f`. +#[inline] unsafe fn model_dolfinx_p2_cell( - f: *const ffi::Function_f64, - simplex_coords: &[f64; 9], - cell: i32, + [p0x, p0y, _, p1x, p1y, _, p2x, p2y, _, p12x, p12y, _, p20x, p20y, _, p01x, p01y, _]: [f64; 18], + [v0, v1, v2, v12, v20, v01]: [f64; 6], max: bool, ) -> (P2LocalModel, Simplex) { // Form simplex - let simplex = Simplex(coords_from_fenics(simplex_coords)); - // Form an array of all coordinates - let midpoints = simplex.midpoints(); - let mut all_coords_uninit: [MaybeUninit; 6 * 3] = [const { MaybeUninit::uninit() }; 6 * 3]; - all_coords_uninit - .iter_mut() - .zip( - simplex_coords - .iter() - .copied() - .chain(coords_to_fenics(midpoints).into_iter()), - ) - .for_each(|(a, b)| { - a.write(b); - }); - let all_coords = MaybeUninit::array_assume_init(all_coords_uninit); - // Evaluate the function, and invert the signs of maximising - let mut values = ffi::eval_Function_f64_cell_6(f, &all_coords, cell); - if max { - values.iter_mut().for_each(|v| *v = -*v); - } + let coords = [Loc([p0x, p0y]), Loc([p1x, p1y]), Loc([p2x, p2y])]; + let simplex = Simplex(coords); + let values = if max { + [-v0, -v1, -v2, -v01, -v12, -v20] + } else { + [v0, v1, v2, v01, v12, v20] + }; // Form local model - let [ref n0, ref n1, ref n2] = simplex.0; - let [ref n01, ref n12, ref n20] = midpoints; - let nodes = [*n0, *n1, *n2, *n01, *n12, *n20]; + let nodes = [Loc([p0x, p0y]), Loc([p1x, p1y]), Loc([p2x, p2y]), Loc([p01x, p01y]), Loc([p12x, p12y]), Loc([p20x, p20y])]; let model = P2LocalModel::::new(&nodes, &values); return (model, simplex); @@ -221,13 +187,12 @@ /// in `simplex_doords`, as well as edge midpoints). Then it minimises this /// model. The result should be accurate, given that `f` is assumed to have a /// P2 model on the simplex. -pub unsafe fn minmax_dolfinx_p2_cell( - f: *const ffi::Function_f64, - simplex_coords: &[f64; 9], - cell: i32, +pub unsafe fn minmax_dolfinx_p2( + simplex_coords: [f64; 18], + vals: [f64; 6], max: bool, ) -> ffi::CoordValuePair { - let (model, simplex) = model_dolfinx_p2_cell(f, simplex_coords, cell, max); + let (model, simplex) = model_dolfinx_p2_cell(simplex_coords, vals, max); let (x, v) = model.minimise(&simplex); CoordValuePair { x: x.into(), v } } diff -r c4f923017349 -r a0d61ac3719e src/dolfinx_access/function.rs --- a/src/dolfinx_access/function.rs Mon Jun 29 00:20:42 2026 +0300 +++ b/src/dolfinx_access/function.rs Sun Oct 04 13:02:26 2026 -0500 @@ -16,6 +16,7 @@ use pyo3::prelude::*; use pyo3::types::{PyDict, PyFunction}; use std::ffi::CStr; +use std::mem::MaybeUninit; /// A helper structure of dealing with dolfinx functions. /// `N` is the domain dimension, `O` the order, and `D` is the codomain dimension. @@ -274,12 +275,17 @@ let (model, _) = unsafe { // Find the coordinates of the corresponding triangle. We guarantee a // triangular mesh in the constructor of `DolfinxPyFunction_f64`. - let simplex_coords = ffi::cell_coords_Function_f64_triangle(f, cell); + //let simplex_coords = ffi::cell_coords_Function_f64_triangle(f, cell); // Construct a new P2 model on that cell, as this is easier than dealing // with the awful documentation of Dolphinx and C++, and probably having // to reimplemen our P2 stuff anyway, to do any optimisation with the basix // basis functions. - super::model_dolfinx_p2_cell(f, &simplex_coords, cell, false) + let coords_uninit: [MaybeUninit; 18] = [const { MaybeUninit::uninit() }; 18]; + let vals_uninit: [MaybeUninit; 6] = [const { MaybeUninit::uninit() }; 6]; + let mut coords = MaybeUninit::array_assume_init(coords_uninit); + let mut vals = MaybeUninit::array_assume_init(vals_uninit); + ffi::cell_dofs_Function_f64_p2(f, cell, &mut coords, &mut vals); + super::model_dolfinx_p2_cell(coords, vals, false) }; model.differential(&x_own) } diff -r c4f923017349 -r a0d61ac3719e src/dolfinx_access/minmax_p2.cc --- a/src/dolfinx_access/minmax_p2.cc Mon Jun 29 00:20:42 2026 +0300 +++ b/src/dolfinx_access/minmax_p2.cc Sun Oct 04 13:02:26 2026 -0500 @@ -83,53 +83,41 @@ // printf("%zd\n", d.size()); // } - // DOF coordinates (may be more than nodes, e.g., intermediate points) - auto gx = geom.x(); - // auto num_points = gx.size() / 3; - // Map from cells to DOF indices - auto geom_dofmap = geom.dofmap(); - assert(geom_dofmap.extent(0) == (size_t)n_cells); - assert(geom_dofmap.extent(1) == 3); + // DOF coordinates and dofmap for the function space itself (not the mesh geometry). + // For a scalar P2 space this gives us, per cell, the 6 nodal points (3 corners + 3 edge + // midpoints) directly in basix's local ordering, so we no longer need to reconstruct + // edge midpoints via f->eval. auto fp_dofmap = fp->dofmap(); - - // auto dof_coords = fp->tabulate_dof_coordinates(); + // transpose=false gives a flat vector laid out as (num_dofs, 3) row-major, i.e. + // dof_coords[3*d + component]. + auto dof_coords = fp->tabulate_dof_coordinates(false); + auto v_array = v->array(); // local coefficient vector CoordValuePair res{{{0.0, 0.0}}, INFINITY}; - // printf("%zd %zd\n", v->array().size(), gx.size()); + for (auto i = 0; i < n_cells; i++) { + auto dofs = fp_dofmap->cell_dofs(i); + assert(dofs.size() == 6); - for (auto i = 0; i < n_cells; i++) { - // Cell corners - auto j0 = geom_dofmap(i, 0); - auto j1 = geom_dofmap(i, 1); - auto j2 = geom_dofmap(i, 2); + // Coordinates of all 6 P2 nodes for this cell, in Fenics' hard-coded 3D layout. + std::array x{}; + // Nodal values of f at those same 6 points. + std::array vals{}; - // Construct list of coordinates. Due to f->eval, we construct this is a single list - // with all the Fenics weirdness of hard-coded 3D. - std::array x = {gx[3 * j0], gx[3 * j0 + 1], 0, - gx[3 * j1], gx[3 * j1 + 1], 0, - gx[3 * j2], gx[3 * j2 + 1], 0}; + for (std::size_t k = 0; k < 6; k++) { + auto d = dofs[k]; + vals[k] = v_array[d]; + x[3 * k] = dof_coords[3 * d]; + x[3 * k + 1] = dof_coords[3 * d + 1]; + x[3 * k + 2] = 0.0; + } // Find the solution on the simplex. - auto newres = minmax_dolfinx_p2_cell(f, x, i, max); + auto newres = minmax_dolfinx_p2(x, vals, max); - // printf("%f [%f, %f] [%f %f %f %f %f %f] <%f, %f> <%f, %f> <‸%f, %f>\n", newres.v, - // newres.x[0], newres.x[1], values[0], values[1], values[2], values[3], - // values[4], values[5], x[0], x[1], x[3], x[4], x[6], x[7]); if (newres.v < res.v) { res = newres; } - - // // Full DOF data - // auto d = fp_dofmap->cell_dofs(i); - // // TODO: array() gives just local part of vector - // assert(d.size() == 6); - // for (auto k = 0; k < d.size(); k++) { - // auto j = d[k]; - // // printf("%d %d\n", 3 * j + 1, nodes.size()); - // assert(2 * j + 1 <= nodes.size()); - // auto x = {nodes[3 * j], nodes[3 * j + 1]}; - // } } // Negate result if maximising. @@ -139,6 +127,25 @@ return res; } + void cell_dofs_Function_f64_p2(Function const* f, int32_t cell, + std::array& x, std::array& vals) { + auto fp = f->function_space(); + auto fp_dofmap = fp->dofmap(); + auto dofs = fp_dofmap->cell_dofs(cell); + assert(dofs.size() == 6); + auto dof_coords = fp->tabulate_dof_coordinates(false); + auto v = f->x(); + auto v_array = v->array(); // local coefficient vector + + for (std::size_t k = 0; k < 6; k++) { + auto d = dofs[k]; + vals[k] = v_array[d]; + x[3 * k] = dof_coords[3 * d]; + x[3 * k + 1] = dof_coords[3 * d + 1]; + x[3 * k + 2] = 0.0; + } + } + CoordValuePair min_Function_f64_p2(Function const* f) { return minmax_Function_f64_p2(f, false); } @@ -154,6 +161,7 @@ m.def("minmax_Function_f64_p2", &minmax_Function_f64_p2); m.def("eval_Function_f64", &eval_Function_f64); m.def("cell_Function_f64", &cell_Function_f64); + m.def("cell_dofs_Function_f64_p2", &cell_dofs_Function_f64_p2); // m.def("cell_Mesh_f64", &cell_Mesh_f64); m.def("cell_FunctionSpace_f64", &cell_FunctionSpace_f64); nanobind::class_(m, "CoordValuePair") diff -r c4f923017349 -r a0d61ac3719e src/experiments.rs --- 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>, + + /// Domain + #[pyo3(set, get)] + domain: Option>, + + /// Randon number generator seed + #[pyo3(set, get)] + seed: Option, } #[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, regterm: RegTermPy, @@ -277,8 +288,10 @@ auxinit: Option>, μinit: Option, f64>>, plot_factory: Option>, + domain: Option>, + seed: Option, ) -> 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>>( + &self, + cli: &CommandLineArgs, + cli_algorithm_overrides: &AlgorithmOverrides, + prefix: String, + problem: PyRef, + 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> = 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 for PythonExperiment { fn name(&self) -> &str { &self.name @@ -295,11 +482,8 @@ fn runall( &self, cli: &CommandLineArgs, - algs: Option>>>, + cli_algorithm_overrides: &AlgorithmOverrides, ) -> 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 = 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> = 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, _> = + 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 { - self.algorithm_overrides - .get(&alg) - .cloned() - .unwrap_or(Default::default()) + }) } } diff -r c4f923017349 -r a0d61ac3719e src/laser_sampling.py --- a/src/laser_sampling.py Mon Jun 29 00:20:42 2026 +0300 +++ b/src/laser_sampling.py Sun Oct 04 13:02:26 2026 -0500 @@ -1,12 +1,15 @@ +import basix import numpy as np -from dolfinx import fem, geometry -from dolfinx_access import cell_FunctionSpace_f64 +from dolfinx import fem from mpi4py import MPI from petsc4py import PETSc from pointsource_pde.dolfinx_extras import get_mass_matrix +from scipy.sparse import coo_matrix from scipy.sparse.linalg import svds from slepc4py import SLEPc +from dolfinx_access import cell_FunctionSpace_f64 + class LaserSampling: """ @@ -20,7 +23,6 @@ nx, ny, beams, - num_segments=100, domain_size=[-2, 2, -2, 2], ): self.V = V @@ -28,7 +30,10 @@ self.nx, self.ny = nx, ny self.M = len(beams) self.domain_size = domain_size - self.num_segments = num_segments + + # Basix element matching V, used to tabulate phi_h exactly at points + # self._element = self._create_basix_element() + self._element = self.V.element.basix_element # Build R matrix (single function does everything) self.R = self.matrix(beams) @@ -60,71 +65,133 @@ self._opnorm = s[0] # / sigma # print(f"Laser sampling opnorm = {self._opnorm:.6f}") - def matrix(self, beams): - """Combined: generate beams + segment + cell collision + dof weighting""" - dofmap = self.V.dofmap - tdim = self.domain.topology.dim + # ------------------------------------------------------------------ + # Exact per-cell line integration + # + # A beam is split into the ordered sub-segments where it actually + # crosses each element (via Cyrus-Beck clipping against the cell's + # polygon), and each sub-segment is integrated exactly using the + # 3-point Simpson rule (exact for phi_h of degree <= 2, since phi_h + # pulled back along a straight line stays degree <= 2 in the + # parameter t). This replaces the previous approach of sampling many + # midpoints and splitting weight equally across a cell's dofs. + # ------------------------------------------------------------------ - # Geometry trees - - bb_tree = geometry.bb_tree(self.domain, tdim) - - n_dofs = self.V.dofmap.index_map.size_global - R = np.zeros((self.M, n_dofs)) + def _cell_vertices_2d(self, cell_idx): + """Physical (x, y) vertices of a cell, ordered CCW around its centroid.""" + x_geom = self.domain.geometry.x + verts = self.domain.geometry.dofmap[cell_idx] + pts = x_geom[verts, :2] + centroid = pts.mean(axis=0) + angles = np.arctan2(pts[:, 1] - centroid[1], pts[:, 0] - centroid[0]) + return pts[np.argsort(angles)] - xmin, xmax, ymin, ymax = self.domain_size + def _clip_segment_to_cell(self, pts, x, y, eps=1e-12): + """ + Cyrus-Beck clip of the segment x->y against the convex polygon `cell_idx`. + Returns (t0, t1) in [0, 1] such that x + t*(y - x), t in [t0, t1], is the + part of the segment inside the cell -- or None if there's no overlap. + Works for any convex cell (triangle or quad); assumes straight-sided + (affine, P1) geometry. + """ + n_verts = len(pts) + d = y - x + t_lo, t_hi = 0.0, 1.0 + + for i in range(n_verts): + v0 = pts[i] + v1 = pts[(i + 1) % n_verts] + edge = v1 - v0 + normal = np.array([edge[1], -edge[0]]) # outward normal (CCW polygon) - for beam_idx, (start_pt, end_pt) in enumerate(beams): - # Segment into midpoints - vec = end_pt - start_pt - total_len = np.linalg.norm(vec) - seg_len = total_len / self.num_segments + denom = normal @ d + numer = normal @ (v0 - x) + + if abs(denom) < eps: + if numer < -eps: # segment parallel to edge, and entirely outside it + return None + continue # parallel and inside: this edge imposes no constraint - midpoints = [] - for i in range(self.num_segments): - t = (i + 0.5) / self.num_segments - midpoint = start_pt + t * vec - midpoints.append(midpoint) - # midpoints = np.array(midpoints, dtype=np.float64) # Force float64 + t_edge = numer / denom + if denom > 0: + t_hi = min(t_hi, t_edge) + else: + t_lo = max(t_lo, t_edge) + + if t_lo > t_hi + eps: + return None + + if t_lo > t_hi: + return None + return t_lo, t_hi - # Ensure shape (N, 3), C-contiguous, read-only for DOLFINx - # assert midpoints.shape[1] == 2, ( - # "Expected 2D points, got shape[1]=2 but need padding to 3D" - # ) - # midpoints = np.pad( - # midpoints, ((0, 0), (0, 1)), mode="constant" - # ) # (N, 3) with z=0 - # midpoints = np.ascontiguousarray(midpoints) # C-order - # midpoints.setflags(write=False) # Read-only requirement - - # Find colliding cells - # cell_candidates = geometry.compute_collisions_points(bb_tree, midpoints) - # colliding_cells = geometry.compute_colliding_cells( - # self.domain, cell_candidates, midpoints - # ) + def matrix(self, beams): + """Combined: exact beam/cell intersection + exact per-cell line integration.""" + n_dofs = self.V.dofmap.index_map.size_global + # R = np.zeros((self.M, n_dofs)) + rows, cols, data = [], [], [] - # Cell→DOFs→weight accumulation - # cell_lists = list(colliding_cells.array) - # for seg_i, cell_idx_raw in enumerate(cell_lists): - # cell_idx = int(cell_idx_raw) # Convert np.int32 → Python int - for point in midpoints: - pad_point = np.array([point[0], point[1], 0.0]) - cell_idx = cell_FunctionSpace_f64(self.V._cpp_object, pad_point) - if cell_idx >= 0: # Valid cell index - # Get global DOFs - cell_dofs_local = dofmap.cell_dofs(cell_idx) - # Convert local dofs to numpy int32 array (DOLFINx requirement) - # local_dofs_array = np.asarray(cell_dofs_local, dtype=np.int32) - # local_dofs_array.setflags(write=False) # Read-only requirement + for cell_idx in range(len(self.domain.geometry.dofmap)): + pts = self._cell_vertices_2d(cell_idx) + for beam_idx, (start_pt, end_pt) in enumerate(beams): + clip = self._clip_segment_to_cell(pts, start_pt, end_pt) + if clip is not None: + t0, t1 = clip + if t1 - t0 > 1e-10: + beam_v = end_pt - start_pt + x0 = start_pt + t0 * beam_v + x1 = start_pt + t1 * beam_v + midpt = start_pt + (0.5 * (t1 + t0)) * beam_v + """ + Exact int_{x0}^{x1} phi_h ds for phi_h in V restricted to `cell_idx`, as a + (global_dofs, row_values) pair such that row_values @ phi_h.x.array[global_dofs] + equals the integral. Uses the 3-point Simpson rule (endpoints + midpoint), + which is exact since phi_h has degree <= 2 along the line. + """ + x_geom = self.domain.geometry.x + verts = self.domain.geometry.dofmap[cell_idx] + v0, v1, v2 = ( + x_geom[verts[0], :2], + x_geom[verts[1], :2], + x_geom[verts[2], :2], + ) + J = np.column_stack([v1 - v0, v2 - v0]) + """Map a physical point p to reference coordinates on an affine triangle.""" + pts_ref = np.array( + [ + np.linalg.solve(J, x0 - v0), + np.linalg.solve(J, midpt - v0), + np.linalg.solve(J, x1 - v0), + ] + ) + tab = self._element.tabulate(0, pts_ref)[ + 0 + ] # (npoints, ndof_local, value_size) + tab = tab[ + ..., 0 + ] # scalar element -> drop trailing value_size=1 axis - # Batch convert all local indices to global in one call - cell_dofs_global = dofmap.index_map.local_to_global(cell_dofs_local) - cell_dofs = np.unique(cell_dofs_global) + length = (t1 - t0) * np.linalg.norm(beam_v) + weights = np.array([1.0, 4.0, 1.0]) * length / 6.0 + local_dofs = self.V.dofmap.cell_dofs(cell_idx) + # global_dofs = self.V.dofmap.index_map.local_to_global(local_dofs) + # R[beam_idx, local_dofs] += weights @ tab + rows.append(np.full(len(local_dofs), beam_idx, dtype=np.int64)) + cols.append(np.asarray(local_dofs, dtype=np.int64)) + data.append(weights @ tab) + + # return R - weight = seg_len / len(cell_dofs) - R[beam_idx, cell_dofs] += weight # ADD for overlaps! + if rows: + rows = np.concatenate(rows) + cols = np.concatenate(cols) + data = np.concatenate(data) + else: + rows = np.array([], dtype=np.int64) + cols = np.array([], dtype=np.int64) + data = np.array([], dtype=np.float64) - return R + return coo_matrix((data, (rows, cols)), shape=(self.M, n_dofs)).tocsr() def apply(self, x): """ @@ -143,15 +210,12 @@ # ⟨x.array,R₀^* v⟩_{ℝ^n} = ⟨Ax.array,[R^*v].array⟩_{ℝ^m} # Taking [R^*v].array=A^{-1}R₀^* v, the conjugate works correctly. tmp = fem.Function(self.V) - np.matmul(self.R.T, j, out=tmp.x.array) + res = fem.Function(self.V) + tmp.x.array[:] = self.R.T @ j tmp.x.scatter_forward() - return tmp - # TODO: is convection_diffusion actualy based on norms in ℝ^n? - # res = fem.Function(self.V) - # res.x.array[:] = 0.0 - # self.solver.solve(tmp.x.petsc_vec, res.x.petsc_vec) - # res.x.scatter_forward() - # return res + self.solver.solve(tmp.x.petsc_vec, res.x.petsc_vec) + res.x.scatter_forward() + return res def opnorm(self, *args): # raise NotImplementedError("Need mesh weighting?")