--- 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)