src/weight_optim.rs

//! Weight optimisation routines

use numeric_literals::replace_float_literals;
use alg_tools::nalgebra_support::ToNalgebraRealField;
use alg_tools::loc::Loc;
use alg_tools::iterate::AlgIteratorFactory;
use crate::types::*;
use crate::subproblem::*;
use crate::measures::*;
use crate::forward_model::ForwardModel;

/// Helper struct for pre-initialising the finite-dimensional subproblems solver
/// [`prepare_optimise_weights_l2`].
///
/// The pre-initialisation is done by [`prepare_optimise_weights_l2`].
pub struct FindimData<F : Float> {
    opAnorm_squared : F
}

/// Return a pre-initialisation struct for [`prepare_optimise_weights_l2`].
///
/// The parameter `opA` is the forward operator $A$.
pub fn prepare_optimise_weights<F, A, const N : usize>(opA : &A) -> FindimData<F>
where F : Float + ToNalgebraRealField,
      A : ForwardModel<Loc<F, N>, F> {
    FindimData{
        opAnorm_squared : opA.opnorm_bound().powi(2)
    }
}

/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
/// point source localisation problem.
///
/// That is, we minimise
/// <div>$$
///     μ ↦ \frac{1}{2}\|Aμ-b\|_2^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
/// $$</div>
/// only with respect to the weights of $μ$.
///
/// The parameter `μ` is the discrete measure whose weights are to be optimised.
/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
/// save intermediate iteration states as images. The parameter `findim_data` should be
/// prepared using [`prepare_optimise_weights`]:
///
/// Returns the number of iterations taken by the method configured in `inner`.
pub fn optimise_weights_l2<'a, F, A, I, const N : usize>(
    μ : &mut DiscreteMeasure<Loc<F, N>, F>,
    opA : &'a A,
    b : &A::Observable,
    α : F,
    findim_data : &FindimData<F>,
    inner : &InnerSettings<F>,
    iterator : I
) -> usize
where F : Float + ToNalgebraRealField,
      I : AlgIteratorFactory<F>,
      A : ForwardModel<Loc<F, N>, F>
{
    // Form and solve finite-dimensional subproblem.
    let (Ã, g̃) = opA.findim_quadratic_model(&μ, b);
    let mut x = μ.masses_dvector();

    // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
    // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
    // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
    // 2-norm, we estimate
    //      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
    //                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
    // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
    // square root is needed when we scale:
    let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len()));
    let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator);
    // Update masses of μ based on solution of finite-dimensional subproblem.
    μ.set_masses_dvector(&x);

    iters
}

/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity
/// point source localisation problem.
///
/// That is, we minimise
/// <div>$$
///     μ ↦ \|Aμ-b\|_1 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
/// $$</div>
/// only with respect to the weights of $μ$.
///
/// The parameter `μ` is the discrete measure whose weights are to be optimised.
/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the
/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while
/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to
/// save intermediate iteration states as images.
///
/// Returns the number of iterations taken by the method configured in `inner`.
#[replace_float_literals(F::cast_from(literal))]
pub fn optimise_weights_l1<'a, F, A, I, const N : usize>(
    μ : &mut DiscreteMeasure<Loc<F, N>, F>,
    opA : &'a A,
    b : &A::Observable,
    α : F,
    findim_data : &FindimData<F>,
    inner : &InnerSettings<F>,
    iterator : I
) -> usize
where F : Float + ToNalgebraRealField,
      I : AlgIteratorFactory<F>,
      A : ForwardModel<Loc<F, N>, F>
{
    // Form and solve finite-dimensional subproblem.
    let Ã = opA.findim_model(&μ);
    let mut x = μ.masses_dvector();

    // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to
    // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the
    // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to
    // 2-norm, we estimate
    //      ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2
    //                 = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2},
    // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no
    // square root is needed when we scale:
    let l = (findim_data.opAnorm_squared * F::cast_from(μ.len())).sqrt();
    let inner_σ = (0.99 / inner.τ0) / l;
    let inner_τ = inner.τ0 / l;
    let iters = l1_nonneg(L1InnerMethod::PDPS, &Ã, b, α, &mut x, inner_τ, inner_σ, iterator);
    // Update masses of μ based on solution of finite-dimensional subproblem.
    μ.set_masses_dvector(&x);

    iters
}