pointsource_pde: comparison experiments/laser_and_mirrors

-:a4137aedcb3a
+:c3a4f4bb87f7
 from pointsource_pde.compose import (
 ComposeFnWithOperator,
 SumOfSeparableFunctions,
 )
 from pointsource_pde.convection_diffusion import (
+BoxedQuadraticRegularisation,
 ConvectionDiffusion,
 PlotFactory,
-QuadraticRegularisation,
 XBound,
 )
 from pointsource_pde.laser_sampling import LaserSampling
 from pointsource_pde.quadratic_dataterm import QuadraticDataTerm
 name = "laser_and_mirrors_aux"
 # 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 = {
 "RadonFB": alg_finetune,
 "RadonSlidingFB": alg_finetune,
 "RadonSlidingPDPS": alg_finetune,
 }
 # Setup shared between experiment with known/unknown conductivity and diffusivity
-def generic_setup(prefix, simple=True):
+def generic_setup(
+prefix,
-# Create true source
+simple=True,
-μ_true = DiscreteMeasure_2_f64(
+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(
 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),
 )
 # Set up lasers and mirrors
 lasers = [
 np.array([0.1, 0.1]),
 v0 = 1.0
 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):
 r = np.sqrt(x[0] ** 2 + x[1] ** 2 + 1e-10)
 return v0 * x[0] / r
 auxtrue = (k, (c1, c2))
 # 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(
 SumOfSeparableFunctions(
 [
 xbound=dat_bound,
 b_norm=norm(b),
 λ=pde.dt / len(beams),
 )
 for data in b
-]
+],
 ),
 pde,
-xbound=[xbound],
+)
-xbound_pair=[xbound],
-)
+# Save true and initial data
-# Plot true and initial data
 plot_factory = PlotFactory(pde)
 plotter = plot_factory.produce(prefix)
 ω, _ = dat.diff((DiscreteMeasure_2_f64([]), auxtrue))
 ω_true, _ = 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,
 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)
-reg = RegTerm.NonnegRadon(10e-7)
+dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup(prefix, rng=rng)
+# 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
+(k, (c1, c2)) = auxtrue
-aux = QuadraticRegularisation(0.001, aux0)
+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("Initial data term value:", inival + ax)
-print("Value at true μ:", dat.apply((μ_true, auxtrue)) + ax)
+print("Data term value at true μ:", dat.apply((μ_true, auxtrue)) + ax)
-# auxinit = np.c_[
+# No curvature bound given: θ is aboslute.
-#     np.zeros_like(k),
+dat.curvature_bound_components = lambda: (None, None)
-#     np.zeros_like(c1),
-#     np.zeros_like(c2),
-# ]
-# print(np.shape(auxinit))
-# 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
 return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)

comparison: experiments/laser_and_mirrors_aux.py

experiments/laser_and_mirrors_aux.py

Mercurial > repos > pointsource_pde / file comparison

comparison: experiments/laser_and_mirrors_aux.py

experiments/laser_and_mirrors_aux.py