diff -r a4137aedcb3a -r c3a4f4bb87f7 experiments/laser_and_mirrors_aux.py
--- 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)