Mercurial > repos > pointsource_pde / file comparison
comparison: experiments/laser_and_mirrors_aux2.py
experiments/laser_and_mirrors_aux2.py
- changeset 3
- c3a4f4bb87f7
equal
deleted
inserted
replaced
| 1:a4137aedcb3a | 3:c3a4f4bb87f7 |
|---|---|
| 1 # | |
| 2 # Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity. | |
| 3 # Alternative parametrisation | |
| 4 # | |
| 5 import laser_and_mirrors_aux | |
| 6 import numpy as np | |
| 7 from laser_and_mirrors_aux import generic_setup, relnoise | |
| 8 from measures import DiscreteMeasure_2_f64 | |
| 9 from pointsource_pde import Problem, RegTerm | |
| 10 from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation | |
| 11 | |
| 12 # Give name to the problem | |
| 13 name = "laser_and_mirrors_aux2" | |
| 14 | |
| 15 # Override algorithm settings | |
| 16 algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides | |
| 17 | |
| 18 | |
| 19 # Setup the problem | |
| 20 def setup(prefix): | |
| 21 rng = np.random.default_rng(seed=31337) | |
| 22 | |
| 23 dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup( | |
| 24 prefix, | |
| 25 rng=rng, | |
| 26 μ_true=DiscreteMeasure_2_f64([([0.2, 0.3], 0.15), ([0.4, 0.1], 0.04)]), | |
| 27 k=0.02, | |
| 28 r=0.1, | |
| 29 θ=120 * np.pi / 180, | |
| 30 ) | |
| 31 | |
| 32 # Override Lipschitz parameter | |
| 33 pde.override_lipschitz = (4.0,) | |
| 34 pde.override_lipschitz_pair = (4.0, 4.0) | |
| 35 # print("diff_chain_lipschitz_factor (modified)", pde.diff_chain_lipschitz_factor()) | |
| 36 # l1, l2 = pde.diff_chain_lipschitz_factor_pair() | |
| 37 # print("diff_chain_lipschitz_factor_pair (modified) ", l1, ", ", l2) | |
| 38 | |
| 39 reg = RegTerm.NonnegRadon(1.5e-6) | |
| 40 | |
| 41 μ0 = DiscreteMeasure_2_f64([]) | |
| 42 (k, (c1, c2)) = auxtrue | |
| 43 aux0 = ( | |
| 44 max(0.001, relnoise(k, 0.02, rng)), | |
| 45 ( | |
| 46 relnoise(c1, 0.2, rng), | |
| 47 relnoise(c2, 0.2, rng), | |
| 48 ), | |
| 49 ) | |
| 50 print("aux init ", aux0) | |
| 51 aux = BoxedQuadraticRegularisation((0.001, -1.0), (1.0, 1.0), (3.0, 0.0005), aux0) | |
| 52 ax = aux.apply(aux0) | |
| 53 inival = dat.apply((μ0, aux0)) | |
| 54 print("Initial data term value:", inival + ax) | |
| 55 print("Data term value at true μ:", dat.apply((μ_true, auxtrue)) + ax) | |
| 56 | |
| 57 # No curvature bound given: θ is aboslute. | |
| 58 dat.curvature_bound_components = lambda: (None, None) | |
| 59 | |
| 60 return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory) |
