experiments/laser_and_mirrors_aux2.py

Fri, 08 May 2026 17:16:34 -0500

author
Tuomo Valkonen <tuomov@iki.fi>
date
Fri, 08 May 2026 17:16:34 -0500
changeset 4
49b062acace9
parent 3
c3a4f4bb87f7
permissions
-rw-r--r--

Do not directly depend on ndarray, but through numpy

#
# Laser and mirrors convection-diffusion experiment with known diffusivity and convectivity.
# Alternative parametrisation
#
import laser_and_mirrors_aux
import numpy as np
from laser_and_mirrors_aux import generic_setup, relnoise
from measures import DiscreteMeasure_2_f64
from pointsource_pde import Problem, RegTerm
from pointsource_pde.convection_diffusion import BoxedQuadraticRegularisation

# Give name to the problem
name = "laser_and_mirrors_aux2"

# Override algorithm settings
algorithm_overrides = laser_and_mirrors_aux.algorithm_overrides


# Setup the problem
def setup(prefix):
    rng = np.random.default_rng(seed=31337)

    dat, auxtrue, _μ_bound, μ_true, plot_factory, pde = generic_setup(
        prefix,
        rng=rng,
        μ_true=DiscreteMeasure_2_f64([([0.2, 0.3], 0.15), ([0.4, 0.1], 0.04)]),
        k=0.02,
        r=0.1,
        θ=120 * np.pi / 180,
    )

    # 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([])
    (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 data term value:", inival + ax)
    print("Data term value at true μ:", dat.apply((μ_true, auxtrue)) + ax)

    # No curvature bound given: θ is aboslute.
    dat.curvature_bound_components = lambda: (None, None)

    return Problem(dat, reg, aux, aux0, μ0, plot_factory=plot_factory)

mercurial