Mercurial > repos > pointsource_pde / file revision

#
# 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, reg_params
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,
    alpha=reg_params["alpha"],
    alpha_k=reg_params["alpha_k"],
    alpha_c=reg_params["alpha_c"],
    c_noise=reg_params["c_noise"],
    k_noise=reg_params["k_noise"],
    noise=reg_params["noise"],
    seed=reg_params["seed"],
    sim_oversample=reg_params["sim_oversample"],
    sim_oversample_time=reg_params["sim_oversample_time"],
    true_init=False,
    true_target=False,
    scale_cstep=reg_params["scale_cstep"],
):
    rng = np.random.default_rng(seed=seed)

    dat, auxtrue, aux0, _μ_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,
        c_noise=c_noise,
        k_noise=k_noise,
        noise=noise,
        sim_oversample=sim_oversample,
        sim_oversample_time=sim_oversample_time,
        scale_cstep=scale_cstep,
    )

    # Override Lipschitz parameter
    dat.override_lipschitz = 1.0
    dat.override_lipschitz_pair = (0.0, 0.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(alpha)

    μ0 = DiscreteMeasure_2_f64([])

    print("aux init ", aux0)
    aux = BoxedQuadraticRegularisation(
        (0.001, -1.0),
        (1.0, 1.0),
        (alpha_k, alpha_c),
        auxtrue if true_target else aux0,
        (1, scale_cstep) if scale_cstep is not None else None,
    )
    print("Initial value:", dat.apply((μ0, aux0)) + aux.apply(aux0))
    print("Value at true μ:", dat.apply((μ_true, auxtrue)) + aux.apply(auxtrue))

    return Problem(
        dat,
        reg,
        aux,
        auxtrue if true_init else aux0,
        μ0,
        plot_factory=plot_factory,
        domain=pde.domain_dims,
    )