src/subproblem/l1squared_nonneg.rs

/*!
Iterative algorithms for solving the finite-dimensional subproblem with constraint.
*/

use itertools::izip;
use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1};
use numeric_literals::replace_float_literals;

use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState};
use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField};
use alg_tools::norms::{Dist, L1};

use super::l1squared_unconstrained::l1squared_prox;
use super::nonneg::nonneg_soft_thresholding;
use super::{InnerMethod, InnerSettings};
use crate::types::*;

/// Return maximum of `dist` and distnce of inteval `[lb, ub]` to zero.
#[replace_float_literals(F::cast_from(literal))]
#[inline]
pub(super) fn max_interval_dist_to_zero<F: Float>(dist: F, lb: F, ub: F) -> F {
    if lb < 0.0 {
        if ub > 0.0 {
            dist
        } else {
            dist.max(-ub)
        }
    } else
    /* lb ≥ 0.0*/
    {
        dist.max(lb)
    }
}

/// Returns the ∞-norm minimal subdifferential of $x ↦ (β/2)|x-y|_1^2 - g^⊤ x + λ\|x\|₁ +δ_{≥0}(x)$ at $x$.
///
/// `v` will be modified and cannot be trusted to contain useful values afterwards.
#[replace_float_literals(F::cast_from(literal))]
fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>(
    y: &Vector<F, Dyn, S1>,
    x: &Vector<F, Dyn, S2>,
    g: &Vector<F, Dyn, S3>,
    λ: F,
) -> F
where
    S1: Storage<F, Dyn>,
    S2: Storage<F, Dyn>,
    S3: Storage<F, Dyn>,
    ShapeConstraint: StridesOk<F, Dyn, U1, S2>,
{
    let mut val = 0.0;
    let tmp = y.dist(x, L1);
    for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) {
        let (mut lb, mut ub) = (-g_i + λ, -g_i + λ);
        if num_traits::abs(x_i - *y_i) < F::EPSILON {
            lb -= tmp;
            ub += tmp
        } else if x_i > *y_i {
            lb += tmp;
            ub += tmp
        } else {
            lb -= tmp;
            ub -= tmp
        }
        if x_i < F::EPSILON {
            lb = F::NEG_INFINITY;
        }
        val = max_interval_dist_to_zero(val, lb, ub);
    }
    val
}

// #[replace_float_literals(F::cast_from(literal))]
// fn lbd_soft_thresholding<F: Float>(v: F, λ: F, b: F) -> F {
//     match (b >= 0.0, v >= b) {
//         (true, false) => b,
//         (true, true) => b.max(v - λ), // soft-to-b from above
//         (false, true) => super::unconstrained::soft_thresholding(v, λ),
//         (false, false) => 0.0.min(b.max(v + λ)), // soft-to-0 with lower bound
//     }
// }

/// Calculates $\prox\_f(x)$ for $f(x)=λ\abs{x-y} + δ\_{≥0}(x)$.
/// This is the first output. The second output is the first output $-y$, i..e, the prox
/// of $f\_0$ for the shift function $f\_0(z) = f(z+y) = λ\abs{z} + δ\_{≥-y}(z)$
/// satisfying
/// $$
///     \prox_f(x) = \prox\_{f\_0}(x - y) + y,
/// $$
/// which is also used internally.
/// The third output indicates whether the output is “locked” to either $0$ (resp. $-y$ in
/// the reformulation), or $y$ ($0$ in the reformulation).
/// The fourth output indicates the next highest $λ$ where such locking would happen.
#[replace_float_literals(F::cast_from(literal))]
#[inline]
fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, Option<F>) {
    let z = x - y; // The shift to f_0
    if -y >= 0.0 {
        if x > λ {
            (x - λ, z - λ, false, Some(x))
        } else {
            (0.0, -y, true, None)
        }
    } else if z < 0.0 {
        if λ < -x {
            (0.0, -y, true, Some(-x))
        } else if z + λ < 0.0 {
            (x + λ, z + λ, false, Some(-z))
        } else {
            (y, 0.0, true, None)
        }
    } else if z > λ {
        (x - λ, z - λ, false, Some(z))
    } else {
        (y, 0.0, true, None)
    }
}

/// Calculate $\prox\_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}\_1\^2 + δ\_{≥0}(x)$.
///
/// To derive an algorithm for this, we can use
/// $$
///     \prox_f(x) = \prox\_{f\_0}(x - y) + y
///     \quad\text{for}\quad
///     f\_0(z)=\frac{β}{2}\norm{z}\_1\^2 + δ\_{≥-y}(z).
/// $$
/// Now, the optimality conditions for $w = \prox\_{f\_0}(x)$ are
/// $$\tag{*}
///     x ∈ w + β\norm{w}\_1\sign w + N\_{≥ -y}(w).
/// $$
/// If we know $\norm{w}\_1$, then this is easily solved by lower-bounded soft-thresholding.
/// We find this by sorting the elements by the distance to the 'locked' lower-bounded
/// soft-thresholding target ($0$ or $-y_i$).
/// Then we loop over this sorted vector, increasing our estimate of $\norm{w}\_1$ as we decide
/// that the soft-thresholding parameter $β\norm{w}\_1$ has to be such that the passed elements
/// will reach their locked value (after which they cannot change anymore, for a larger
/// soft-thresholding parameter. This has to be slightly more fine-grained for account
/// for the case that $-y\_i<0$ and $x\_i < -y\_i$.
///
/// Indeed, denoting by $x'$ and $w'$ the subset of elements such that $w\_i ≠ 0$ and
/// $w\_i > -y\_i$, we can calculate by applying $⟨\cdot, \sign w'⟩$ to the corresponding
/// lines of (*) that
/// $$
///     \norm{x'} = \norm{w'} + β \norm{w}\_1 m,
/// $$
/// where $m$ is the number of unlocked components.
/// We can now calculate the mass $t = \norm{w}-\norm{w'}$ of locked components, and so obtain
/// $$
///     \norm{x'} + t = (1+β m)\norm{w}\_1,
/// $$
/// from where we can calculate the soft-thresholding parameter $λ=β\norm{w}\_1$.
/// Since we do not actually know the unlocked elements, but just loop over all the possibilities
/// for them, we have to check that $λ$ is above the current lower bound for this parameter
/// (`shift` in the code), and below a value that would cause changes in the locked set
/// (`max_shift` in the code).
#[replace_float_literals(F::cast_from(literal))]
pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField, S1, S2>(
    x: &mut Vector<F, Dyn, S1>,
    y: &Vector<F, Dyn, S2>,
    β: F,
) where
    S2: Storage<F, Dyn>,
    S1: StorageMut<F, Dyn>,
{
    // nalgebra double-definition bullshit workaround
    let abs = alg_tools::NumTraitsFloat::abs;
    let min = alg_tools::NumTraitsFloat::min;

    let mut λ = 0.0;
    loop {
        let mut w_locked = 0.0;
        let mut n_unlocked = 0; // m
        let mut x_prime = 0.0;
        let mut max_shift = F::INFINITY;
        for (&x_i, &y_i) in izip!(x.iter(), y.iter()) {
            let (_, a_shifted, locked, next_lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ);

            if let Some(t) = next_lock {
                max_shift = min(max_shift, t);
                assert!(max_shift > λ);
            }
            if locked {
                w_locked += abs(a_shifted);
            } else {
                n_unlocked += 1;
                x_prime += abs(x_i - y_i);
            }
        }

        // We need ‖x'‖ = ‖w'‖ + β m ‖w‖, i.e. ‖x'‖ + (‖w‖-‖w'‖)= (1 + β m)‖w‖.
        let λ_new = (x_prime + w_locked) / (1.0 / β + F::cast_from(n_unlocked));

        if λ_new > max_shift {
            λ = max_shift;
        } else {
            assert!(λ_new >= λ);
            // success
            x.zip_apply(y, |x_i, y_i| {
                let (a, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new);
                //*x_i = y_i + lbd_soft_thresholding(*x_i, λ_new, -y_i)
                *x_i = a;
            });
            return;
        }
    }
}

/// Proximal point method implementation of [`l1squared_nonneg`].
/// For detailed documentation of the inputs and outputs, refer to there.
///
/// The `λ` component of the model is handled in the proximal step instead of the gradient step
/// for potential performance improvements.
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
pub fn l1squared_nonneg_pp<F, I, S1, S2, S3>(
    y: &Vector<F::MixedType, Dyn, S1>,
    g: &Vector<F::MixedType, Dyn, S2>,
    λ_: F,
    x: &mut Vector<F::MixedType, Dyn, S3>,
    τ_: F,
    θ_: F,
    iterator: I,
) -> usize
where
    F: Float + ToNalgebraRealField,
    I: AlgIteratorFactory<F>,
    S1: Storage<F::MixedType, Dyn>,
    S2: Storage<F::MixedType, Dyn>,
    S3: StorageMut<F::MixedType, Dyn>,
    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
    let λ = λ_.to_nalgebra_mixed();
    let mut τ = τ_.to_nalgebra_mixed();
    let θ = θ_.to_nalgebra_mixed();
    let mut iters = 0;

    iterator.iterate(|state| {
        // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2+δ_{≥0}+}(x^k - τ(λ𝟙^⊤-g))
        x.apply(|x_i| *x_i -= τ * λ);
        x.axpy(τ, g, 1.0);
        l1squared_nonneg_prox(x, y, τ);

        iters += 1;
        // This gives O(1/N^2) rates due to monotonicity of function values.
        // Higher acceleration does not seem to be numerically stable.
        τ += θ;

        // This gives O(1/N^3) rates due to monotonicity of function values.
        // Higher acceleration does not seem to be numerically stable.
        //τ + = F::cast_from(iters).to_nalgebra_mixed()*θ;

        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
    });

    iters
}

/// PDPS implementation of [`l1squared_nonneg`].
/// For detailed documentation of the inputs and outputs, refer to there.
///
/// The `λ` component of the model is handled in the proximal step instead of the gradient step
/// for potential performance improvements.
/// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model.
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
pub fn l1squared_nonneg_pdps<F, I, S1, S2, S3>(
    y: &Vector<F::MixedType, Dyn, S1>,
    g: &Vector<F::MixedType, Dyn, S2>,
    λ_: F,
    x: &mut Vector<F::MixedType, Dyn, S3>,
    τ_: F,
    σ_: F,
    θ_: F,
    iterator: I,
) -> usize
where
    F: Float + ToNalgebraRealField,
    I: AlgIteratorFactory<F>,
    S1: Storage<F::MixedType, Dyn>,
    S2: Storage<F::MixedType, Dyn>,
    S3: StorageMut<F::MixedType, Dyn>,
    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
    let λ = λ_.to_nalgebra_mixed();
    let τ = τ_.to_nalgebra_mixed();
    let σ = σ_.to_nalgebra_mixed();
    let θ = θ_.to_nalgebra_mixed();
    let mut w = DVector::zeros(x.len());
    let mut tmp = DVector::zeros(x.len());
    let mut xprev = x.clone_owned();
    let mut iters = 0;

    iterator.iterate(|state| {
        // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2}(x^k - τ (w^k - g))
        x.axpy(-τ * θ, &w, 1.0);
        x.axpy(τ, g, 1.0);
        l1squared_prox(&mut tmp, x, y, τ);

        // Dual step: w^{k+1} = proj_{[-∞,λ]}(w^k + σ(2x^{k+1}-x^k))
        w.axpy(2.0 * σ * θ, x, 1.0);
        w.axpy(-σ * θ, &xprev, 1.0);
        w.apply(|w_i| *w_i = w_i.min(λ));
        xprev.copy_from(x);

        iters += 1;

        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
    });

    iters
}

/// Alternative PDPS implementation of [`l1squared_nonneg`].
/// For detailed documentation of the inputs and outputs, refer to there.
///
/// By not dualising the 1-norm, this should produce more sparse solutions than
/// [`l1squared_nonneg_pdps`].
///
/// The `λ` component of the model is handled in the proximal step instead of the gradient step
/// for potential performance improvements.
/// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model.
/// We rewrite
/// <div>$$
///     \begin{split}
///     & \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) \\
///     & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x)
///      - \left(x ↦ \frac{1}{2θ} |x-y|_1^2 \right)^*(w).
///     \end{split}
/// $$</div>
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
pub fn l1squared_nonneg_pdps_alt<F, I, S1, S2, S3>(
    y: &Vector<F::MixedType, Dyn, S1>,
    g: &Vector<F::MixedType, Dyn, S2>,
    λ_: F,
    x: &mut Vector<F::MixedType, Dyn, S3>,
    τ_: F,
    σ_: F,
    θ_: F,
    iterator: I,
) -> usize
where
    F: Float + ToNalgebraRealField,
    I: AlgIteratorFactory<F>,
    S1: Storage<F::MixedType, Dyn>,
    S2: Storage<F::MixedType, Dyn>,
    S3: StorageMut<F::MixedType, Dyn>,
    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
    let λ = λ_.to_nalgebra_mixed();
    let τ = τ_.to_nalgebra_mixed();
    let σ = σ_.to_nalgebra_mixed();
    let θ = θ_.to_nalgebra_mixed();
    let σθ = σ * θ;
    let τλ = τ * λ;
    let one_div_σθ = 1.0 / σθ;
    let mut w = DVector::zeros(x.len());
    let mut tmp = DVector::zeros(x.len());
    let mut xprev = x.clone_owned();
    let mut iters = 0;

    let mut y_scale = y.clone_owned();
    y_scale *= σ;

    iterator.iterate(|state| {
        // Primal step: x^{k+1} = nonnegsoft_τλ(x^k - τ(θ w^k -g))
        if θ == 1.0 {
            for (x_i, xprev_i, w_i, &g_i) in
                izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter())
            {
                *x_i = nonneg_soft_thresholding(*x_i - τ * (*w_i - g_i), τλ);
                // Fused dual part from below
                *w_i += σ * (2.0 * *x_i - *xprev_i);
                *xprev_i = *w_i;
            }
        } else {
            for (x_i, xprev_i, w_i, &g_i) in
                izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter())
            {
                *x_i = nonneg_soft_thresholding(*x_i - τ * (θ * *w_i - g_i), τλ);
                // Fused dual part from below
                *w_i += σ * (2.0 * *x_i - *xprev_i);
                *xprev_i = *w_i;
            }
        }
        // This is numerically unstable:
        // x.axpy(-τθ, &w, 1.0);
        // x.axpy(τ, g, 1.0);
        // x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ * λ));

        // Dual step: with g(x) = (β/(2θ))‖x-y‖₁² and q = w^k + σ(2x^{k+1}-x^k),
        // we compute w^{k+1} = prox_{σg^*}(q) for
        //                    = q - σ prox_{g/σ}(q/σ)
        //                    = q - σ prox_{(β/(2θσ))‖.-y‖₁²}(q/σ)
        //                    = σ(q/σ - prox_{(β/(2θσ))‖.-y‖₁²}(q/σ))
        // ALT
        //                    = q - prox_{(β/(2θσ))‖.-σy‖₁²}(q)
        // where q/σ = w^k/σ + (2x^{k+1}-x^k),

        // This has been fused into the loop below
        // for (xprev_i, w_i, &x_i) in izip!(xprev.iter_mut(), w.iter_mut(), x.iter()) {
        //     *w_i += σ * (2.0 * x_i - *xprev_i);
        //     *xprev_i = *w_i;
        // }
        // xprev.axpy(2.0, x, -1.0);
        // w.axpy(σ, &xprev, 1.0);
        // xprev.copy_from(&w); // use xprev as temporary variable
        //l1squared_prox(&mut tmp, &mut xprev, &y_scale, β_div_σθ);
        l1squared_prox(&mut tmp, &mut xprev, &y_scale, one_div_σθ);
        w -= &xprev;
        xprev.copy_from(x);

        iters += 1;

        state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
    });

    iters
}

/// This is an exact solver for
/// <div>$$
///     \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x).
/// $$</div>
/// i.e.,
/// <div>$$
///     \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 + (λ - g)^⊤ x + δ_{≥ 0}(x),
/// $$</div>
/// which has the optimality conditions
/// <div>$$
///     0 ∈ |x-y|_1\sign(x-y)_i + λ - g_i + δ_{≥ 0}(x_i)
///     \quad\text{for all}\quad i.
/// $$</div>
///
/// If $x_i > 0$ and $x_i ≠ y_i$, then this forces $|g_i -λ| = |x-y|_1$.
/// Let $i^*$ be such an index. Then if $i ≠ i^*$ does not have the same value of $|g_i-λ|$,
/// we *must* have $x_i=y_i$ or $x_i=0$. In fact, even $x_i=y_i$ is impossible for $y_i>0$ if
/// $|g_i-λ| > |g_{i^*}-λ|$, whereas, otherwise $x_i=y_i$ is possible if $y_i>=0$.
/// If $|g_i-λ| > |g_{i^*}-λ|$, we must, therefore, either have $x_i=0$ or be able to take $i^*=i$.
///
/// If $|g_i-λ| ≤ |g_{i^*}-λ|$, we can have either $x_i=0$ or $x_i=y_i$. If $y_i<=0$, then,
/// of course $x_i=0$. If $y_i>0$, we can take $x_i=y_i$.
///
/// We, therefore, sort the |g_i-λ|, and look for an index $i^*$ that gives a non-contradictory
/// $β|x-y|_1=|g_i -λ|$, noting that contributions to β|x-y|_1= come from, besides $i^*$, from
/// indices $i$ such that $x_i= 0 ≠ y_i$.
///
/// Finally, if no index $i$ reaches $|g_i -λ| = |x-y|_1$, we have to $|x-y|_1$ being in an
/// intermediate interval. This can similarly done using sorting.
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
pub(super) fn l1squared_nonneg_solve_exact<F, S1, S2, S3>(
    y: &Vector<F::MixedType, Dyn, S1>,
    g: &Vector<F::MixedType, Dyn, S2>,
    λ_: F,
    x: &mut Vector<F::MixedType, Dyn, S3>,
) -> bool
where
    F: Float + ToNalgebraRealField,
    S1: Storage<F::MixedType, Dyn>,
    S2: Storage<F::MixedType, Dyn>,
    S3: StorageMut<F::MixedType, Dyn>,
    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
    let λ = λ_.to_nalgebra_mixed();

    assert_eq!(y.len(), g.len());
    assert_eq!(x.len(), g.len());
    assert!(x.len() <= u32::MAX as usize);

    #[derive(Debug)]
    struct Tmp<F> {
        y: F,
        d: F,
        contrib: F,
        i: u32, // Shrink size for sort
        λ_le_g: bool,
    }
    let mut sorted = Vec::from_iter(izip!(g.iter(), y.iter(), 0..).map(|(&g_i, &y_i, i)| {
        // We already precompute the comparison λ <= g_i here, since we need it for the abs,
        // and would need to store g_i for that in any case.
        let (d, λ_le_g) = if λ <= g_i {
            (g_i - λ, true)
        } else {
            (λ - g_i, false)
        };
        Tmp { d, λ_le_g, y: y_i, i, contrib: 0.0 }
    }));
    sorted
        .as_mut_slice()
        .sort_unstable_by(|a, b| b.d.total_cmp(&a.d));

    // Reverse-compute contribs
    sorted.iter_mut().rev().fold(0.0, |contrib, a| {
        a.contrib = contrib;
        if a.y < 0.0 {
            contrib - a.y
        } else {
            contrib
        }
    });

    let mut contrib0 = 0.0;

    x.fill(0.0);

    let mut it = sorted.iter();
    let mut found = false;

    // We first try to find an index m such that y_m ≠ x_m > 0.0 that satisfies ‖x-y‖₁ = |λ - g_m|.
    'search: while let Some(m) = it.next() {
        let d_m = m.d;
        let y_m = m.y;
        let contrib = m.contrib + contrib0;
        let δ = d_m - contrib;
        if δ > 0.0 {
            // If λ < g[m], we must have x[m]≥y[m] so (contrib + x[m]-y[m])=g[m]-λ=d
            // If λ = g_m, then also d_m=0, so contrib=0, and this gtives *x_m=y_m.
            // If λ > g[m], We must have x[m]≤y[m] so -(contrib + y[m]-x[m])=g[m]-λ=-d.
            let x_m_prime = if m.λ_le_g { y_m + δ } else { y_m - δ };
            if x_m_prime > 0.0 {
                let x_m = unsafe { x.get_unchecked_mut(m.i as usize) };
                *x_m = x_m_prime;
                found = true;
                break 'search;
            }
        }
        contrib0 += num_traits::abs(y_m);
    }

    if !found {
        // Every x_i is either zero or equal to y_i.
        // We scan intervals (bound,prev_bound) between the d-values of each component,
        // for the value of ‖x-y‖₁. Then it can be decided whether x_i should be one or equal to
        // y_i. Due to the sorting and pre-calculation of posterior contributions, this simplifies
        // into the following:
        let mut prev_bound = F::MixedType::INFINITY;
        let mut contrib0 = 0.0;
        it = sorted.iter();
        'search_degenerate: while let Some(m) = it.next() {
            let bound = m.d;
            let contrib = m.contrib + contrib0 - m.y.min(0.0);
            if prev_bound >= contrib && contrib >= bound {
                if m.y > 0.0 {
                    let x_m = unsafe { x.get_unchecked_mut(m.i as usize) };
                    *x_m = m.y;
                }
                found = true;
                break 'search_degenerate;
            }
            prev_bound = bound;
            if m.y > 0.0 {
                contrib0 += m.y;
            }
        }
        // Check final interval from last element to -∞, if nothing found.
        if !found && !(prev_bound >= contrib0) {
            panic!("l1squared_nonneg_solve_exact failure")
        }
    }

    // Set remaining components to their now known values
    while let Some(a) = it.next() {
        // Safety: size checked above.
        let y_i = a.y;
        if y_i >= 0.0 {
            // We can leave unchanged, as ‖x-y‖₁ is guaranteed large enough if
            // the current maximising index attempt succeeds.
            let x_i = unsafe { x.get_unchecked_mut(a.i as usize) };
            *x_i = y_i
        }
        // Zero fill in `else` case done above in init already
    }

    debug_assert!(min_subdifferential(y, x, g, λ) <= F::EPSILON.to_nalgebra_mixed() * 10.0);

    return true;
}

/// This function applies an iterative method for the solution of the problem
/// <div>$$
///     \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x).
/// $$</div>
///
/// This function returns the number of iterations taken.
#[replace_float_literals(F::cast_from(literal))]
pub fn l1squared_nonneg<F, I, S1, S2, S3>(
    y: &Vector<F::MixedType, Dyn, S1>,
    g: &Vector<F::MixedType, Dyn, S2>,
    λ: F,
    x: &mut Vector<F::MixedType, Dyn, S3>,
    inner: &InnerSettings<F>,
    iterator: I,
) -> usize
where
    F: Float + ToNalgebraRealField,
    I: AlgIteratorFactory<F>,
    S1: Storage<F::MixedType, Dyn>,
    S2: Storage<F::MixedType, Dyn>,
    S3: StorageMut<F::MixedType, Dyn>,
    ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
    if let InnerMethod::Exact = inner.method {
        // Try exact solution, fall back to PDPS if it does not work.
        if l1squared_nonneg_solve_exact(y, g, λ, x) {
            return 1;
        }
    }

    match inner.method {
        InnerMethod::PDPS | InnerMethod::Exact => {
            let inner_θ = 1.0;
            //Estimate of ‖K‖ for K=θ\Id.
            let normest = inner_θ;
            let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest);
            l1squared_nonneg_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator)
        }
        InnerMethod::PP | InnerMethod::FB => {
            let inner_τ = inner.pp_τ.0;
            let inner_θ = inner.pp_τ.1;
            l1squared_nonneg_pp(y, g, λ, x, inner_τ, inner_θ, iterator)
        }
        other => unimplemented!("${other:?} is unimplemented"),
    }
}