Mercurial > repos > pointsource_algs / file revision

//! Implementation of the hat function

use numeric_literals::replace_float_literals;
use serde::Serialize;
use alg_tools::types::*;
use alg_tools::norms::*;
use alg_tools::loc::Loc;
use alg_tools::sets::Cube;
use alg_tools::bisection_tree::{
    Support,
    Constant,
    Bounds,
    LocalAnalysis,
    GlobalAnalysis,
    Bounded,
};
use alg_tools::mapping::{Mapping, Instance};
use alg_tools::maputil::array_init;
use crate::types::Lipschitz;

/// Representation of the hat function $f(x)=1-\\|x\\|\_1/ε$ of `width` $ε$ on $ℝ^N$.
#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
pub struct Hat<C : Constant, const N : usize> {
    /// The parameter $ε>0$.
    pub width : C,
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize> Mapping<Loc<C::Type, N>> for Hat<C, N> {
    type Codomain = C::Type;

    #[inline]
    fn apply<I : Instance<Loc<C::Type, N>>>(&self, x : I) -> Self::Codomain {
        let ε = self.width.value();
        0.0.max(1.0-x.cow().norm(L1)/ε)
    }
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize> Support<C::Type, N> for Hat<C, N> {
    #[inline]
    fn support_hint(&self) -> Cube<C::Type,N> {
        let ε = self.width.value();
        array_init(|| [-ε, ε]).into()
    }

    #[inline]
    fn in_support(&self, x : &Loc<C::Type,N>) -> bool {
        x.norm(L1) < self.width.value()
    }
    
    /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool {
        todo!("Not implemented, but not used at the moment")
    }*/

    #[inline]
    fn bisection_hint(&self, cube : &Cube<C::Type,N>) -> [Option<C::Type>; N] {
        let ε = self.width.value();
        cube.map(|a, b| {
            if a < 1.0 {
                if 1.0 < b {
                    Some(1.0)
                } else {
                    if a < -ε {
                        if b > -ε { Some(-ε) } else { None }
                    } else {
                        None
                    }
                }
            } else {
                if b > ε { Some(ε) } else { None }
            }
        });
        todo!("also diagonals")
    }
}


#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize>
GlobalAnalysis<C::Type, Bounds<C::Type>>
for Hat<C, N> {
    #[inline]
    fn global_analysis(&self) -> Bounds<C::Type> {
        Bounds(0.0, 1.0)
    }
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize> Lipschitz<L1> for Hat<C, N> {
    type FloatType = C::Type;

    fn lipschitz_factor(&self, _l1 : L1) -> Option<C::Type> {
        Some(1.0/self.width.value())
    }
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize> Lipschitz<L2> for Hat<C, N> {
    type FloatType = C::Type;

    fn lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
        self.lipschitz_factor(L1).map(|l1|
            <L2 as Dominated<C::Type, L1, Loc<C::Type,N>>>::from_norm(&L2, l1, L1)
        )
    }
}

impl<'a, C : Constant, const N : usize>
LocalAnalysis<C::Type, Bounds<C::Type>, N>
for Hat<C, N> {
    #[inline]
    fn local_analysis(&self, cube : &Cube<C::Type, N>) -> Bounds<C::Type> {
        // The function is maximised/minimised where the 1-norm is minimised/maximised.
        let lower = self.apply(cube.maxnorm_point());
        let upper = self.apply(cube.minnorm_point());
        Bounds(lower, upper)
    }
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C : Constant, const N : usize>
Norm<C::Type, Linfinity>
for Hat<C, N> {
    #[inline]
    fn norm(&self, _ : Linfinity) -> C::Type {
        self.bounds().upper()
    }
}