Mercurial > repos > pointsource_pde / file comparison
comparison: experiments/laser_and_mirrors_aux2.py
experiments/laser_and_mirrors_aux2.py
- changeset 12
- a0d61ac3719e
- parent 10
- 0c9e19020217
equal
deleted
inserted
replaced
| 11:c4f923017349 | 12:a0d61ac3719e |
|---|---|
| 2 # Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity. | 2 # Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity. |
| 3 # Alternative parametrisation | 3 # Alternative parametrisation |
| 4 # | 4 # |
| 5 import laser_and_mirrors_aux | 5 import laser_and_mirrors_aux |
| 6 import numpy as np | 6 import numpy as np |
| 7 from laser_and_mirrors_aux import generic_setup, relnoise | 7 from laser_and_mirrors_aux import generic_setup, reg_params |
| 8 from measures import DiscreteMeasure_2_f64 | 8 from measures import DiscreteMeasure_2_f64 |
| 9 from pointsource_pde import Problem, RegTerm | 9 from pointsource_pde import Problem, RegTerm |
| 10 from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation | 10 from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation |
| 11 | 11 |
| 12 # Give name to the problem | 12 # Give name to the problem |
| 15 # Override algorithm settings | 15 # Override algorithm settings |
| 16 algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides | 16 algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides |
| 17 | 17 |
| 18 | 18 |
| 19 # Setup the problem | 19 # Setup the problem |
| 20 def setup(prefix): | 20 def setup( |
| 21 rng = np.random.default_rng(seed=31337) | 21 prefix, |
| 22 alpha=reg_params["alpha"], | |
| 23 alpha_k=reg_params["alpha_k"], | |
| 24 alpha_c=reg_params["alpha_c"], | |
| 25 c_noise=reg_params["c_noise"], | |
| 26 k_noise=reg_params["k_noise"], | |
| 27 noise=reg_params["noise"], | |
| 28 seed=reg_params["seed"], | |
| 29 sim_oversample=reg_params["sim_oversample"], | |
| 30 sim_oversample_time=reg_params["sim_oversample_time"], | |
| 31 true_init=False, | |
| 32 true_target=False, | |
| 33 scale_cstep=reg_params["scale_cstep"], | |
| 34 ): | |
| 35 rng = np.random.default_rng(seed=seed) | |
| 22 | 36 |
| 23 dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup( | 37 dat, auxtrue, aux0, _μ_bound, μ_true, plot_factory, pde = generic_setup( |
| 24 prefix, | 38 prefix, |
| 25 rng=rng, | 39 rng=rng, |
| 26 μ_true=DiscreteMeasure_2_f64([([0.2, 0.3], 0.15), ([0.4, 0.1], 0.04)]), | 40 μ_true=DiscreteMeasure_2_f64([([0.2, 0.3], 0.15), ([0.4, 0.1], 0.04)]), |
| 27 k=0.02, | 41 k=0.02, |
| 28 r=0.1, | 42 r=0.1, |
| 29 θ=120 * np.pi / 180, | 43 θ=120 * np.pi / 180, |
| 44 c_noise=c_noise, | |
| 45 k_noise=k_noise, | |
| 46 noise=noise, | |
| 47 sim_oversample=sim_oversample, | |
| 48 sim_oversample_time=sim_oversample_time, | |
| 49 scale_cstep=scale_cstep, | |
| 30 ) | 50 ) |
| 31 | 51 |
| 32 # Override Lipschitz parameter | 52 # Override Lipschitz parameter |
| 33 pde.override_lipschitz = (4.0,) | 53 dat.override_lipschitz = 1.0 |
| 34 pde.override_lipschitz_pair = (4.0, 4.0) | 54 dat.override_lipschitz_pair = (0.0, 0.0) |
| 35 # print("diff_chain_lipschitz_factor (modified)", pde.diff_chain_lipschitz_factor()) | 55 # print("diff_chain_lipschitz_factor (modified)", pde.diff_chain_lipschitz_factor()) |
| 36 # l1, l2 = pde.diff_chain_lipschitz_factor_pair() | 56 # l1, l2 = pde.diff_chain_lipschitz_factor_pair() |
| 37 # print("diff_chain_lipschitz_factor_pair (modified) ", l1, ", ", l2) | 57 # print("diff_chain_lipschitz_factor_pair (modified) ", l1, ", ", l2) |
| 38 | 58 |
| 39 reg = RegTerm.NonnegRadon(1.5e-6) | 59 reg = RegTerm.NonnegRadon(alpha) |
| 40 | 60 |
| 41 μ0 = DiscreteMeasure_2_f64([]) | 61 μ0 = DiscreteMeasure_2_f64([]) |
| 42 (k, (c1, c2)) = auxtrue | 62 |
| 43 aux0 = ( | 63 print("aux init ", aux0) |
| 44 max(0.001, relnoise(k, 0.02, rng)), | 64 aux = BoxedQuadraticRegularisation( |
| 45 ( | 65 (0.001, -1.0), |
| 46 relnoise(c1, 0.2, rng), | 66 (1.0, 1.0), |
| 47 relnoise(c2, 0.2, rng), | 67 (alpha_k, alpha_c), |
| 48 ), | 68 auxtrue if true_target else aux0, |
| 69 (1, scale_cstep) if scale_cstep is not None else None, | |
| 49 ) | 70 ) |
| 50 print("aux init ", aux0) | 71 print("Initial value:", dat.apply((μ0, aux0)) + aux.apply(aux0)) |
| 51 aux = BoxedQuadraticRegularisation((0.001, -1.0), (1.0, 1.0), (3.0, 0.0005), aux0) | 72 print("Value at true μ:", dat.apply((μ_true, auxtrue)) + aux.apply(auxtrue)) |
| 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 | 73 |
| 57 return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory) | 74 return Problem( |
| 75 dat, | |
| 76 reg, | |
| 77 aux, | |
| 78 auxtrue if true_init else aux0, | |
| 79 μ0, | |
| 80 plot_factory=plot_factory, | |
| 81 domain=pde.domain_dims, | |
| 82 ) |
