unthreaded/src/norms.rs

/*!
Norms, projections, etc.
*/

use serde::Serialize;
use crate::types::*;
use crate::euclidean::*;

//
// Abstract norms
//

/// An exponent for norms.
///
// Just a collection of desirabl attributes for a marker type
pub trait NormExponent : Copy + Send + Sync + 'static {}


/// Exponent type for the 1-[`Norm`].
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct L1;
impl NormExponent for L1 {}

/// Exponent type for the 2-[`Norm`].
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct L2;
impl NormExponent for L2 {}

/// Exponent type for the ∞-[`Norm`].
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct Linfinity;
impl NormExponent for Linfinity {}

/// Exponent type for 2,1-[`Norm`].
/// (1-norm over a domain Ω, 2-norm of a vector at each point of the domain.)
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct L21;
impl NormExponent for L21 {}

/// A Huber/Moreau–Yosida smoothed [`L1`] norm. (Not a norm itself.)
///
/// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher
/// values more smoothing. Behaviour with γ < 0 is undefined.
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct HuberL1<F : Float>(pub F);
impl<F : Float> NormExponent for HuberL1<F> {}

/// A Huber/Moreau–Yosida smoothed [`L21`] norm. (Not a norm itself.)
///
/// The parameter γ of this type is the smoothing factor. Zero means no smoothing, and higher
/// values more smoothing. Behaviour with γ < 0 is undefined.
#[derive(Copy,Debug,Clone,Serialize,Eq,PartialEq)]
pub struct HuberL21<F : Float>(pub F);
impl<F : Float> NormExponent for HuberL21<F> {}


/// A normed space (type) with exponent or other type `Exponent` for the norm.
///
/// Use as
/// ```
/// # use alg_tools::norms::{Norm, L1, L2, Linfinity};
/// # use alg_tools::loc::Loc;
/// let x = Loc([1.0, 2.0, 3.0]);
///
/// println!("{}, {} {}", x.norm(L1), x.norm(L2), x.norm(Linfinity))
/// ```
pub trait Norm<F : Num, Exponent : NormExponent> {
    /// Calculate the norm.
    fn norm(&self, _p : Exponent) -> F;
}

/// Indicates that the `Self`-[`Norm`] is dominated by the `Exponent`-`Norm` on the space
/// `Elem` with the corresponding field `F`.
pub trait Dominated<F : Num, Exponent : NormExponent, Elem> {
    /// Indicates the factor $c$ for the inequality $‖x‖ ≤ C ‖x‖_p$.
    fn norm_factor(&self, p : Exponent) -> F;
    /// Given a norm-value $‖x‖_p$, calculates $C‖x‖_p$ such that $‖x‖ ≤ C‖x‖_p$
    #[inline]
    fn from_norm(&self, p_norm : F, p : Exponent) -> F {
        p_norm * self.norm_factor(p)
    }
}

/// Trait for distances with respect to a norm.
pub trait Dist<F : Num, Exponent : NormExponent> : Norm<F, Exponent> {
    /// Calculate the distance
    fn dist(&self, other : &Self, _p : Exponent) -> F;
}

/// Trait for Euclidean projections to the `Exponent`-[`Norm`]-ball.
///
/// Use as
/// ```
/// # use alg_tools::norms::{Projection, L2, Linfinity};
/// # use alg_tools::loc::Loc;
/// let x = Loc([1.0, 2.0, 3.0]);
///
/// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity));
/// ```
pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Euclidean<F>
where F : Float {
    /// Projection of `self` to the `q`-norm-ball of radius ρ.
    fn proj_ball(mut self, ρ : F, q : Exponent) -> Self {
        self.proj_ball_mut(ρ, q);
        self
    }

    /// In-place projection of `self` to the `q`-norm-ball of radius ρ.
    fn proj_ball_mut(&mut self, ρ : F, _q : Exponent);
}

/*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E {
    #[inline]
    fn norm(&self, _p : L2) -> F { self.norm2() }

    fn dist(&self, other : &Self, _p : L2) -> F { self.dist2(other) }
}*/

impl<F : Float, E : Euclidean<F> + Norm<F, L2>> Projection<F, L2> for E {
    #[inline]
    fn proj_ball(self, ρ : F, _p : L2) -> Self { self.proj_ball2(ρ) }

    #[inline]
    fn proj_ball_mut(&mut self, ρ : F, _p : L2) { self.proj_ball2_mut(ρ) }
}

impl<F : Float> HuberL1<F> {
    fn apply(self, xnsq : F) -> F {
        let HuberL1(γ) = self;
        let xn = xnsq.sqrt();
        if γ == F::ZERO {
            xn
        } else {
            if xn > γ {
                xn-γ / F::TWO
            } else if xn<(-γ) {
                -xn-γ / F::TWO
            } else {
                xnsq / (F::TWO * γ)
            }
        }
    }
}

impl<F : Float, E : Euclidean<F>> Norm<F, HuberL1<F>> for E {
    fn norm(&self, huber : HuberL1<F>) -> F {
        huber.apply(self.norm2_squared())
    }
}

impl<F : Float, E : Euclidean<F>> Dist<F, HuberL1<F>> for E {
    fn dist(&self, other : &Self, huber : HuberL1<F>) -> F {
        huber.apply(self.dist2_squared(other))
    }
}