src/kernels/hat.rs

//! Implementation of the hat function

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

/// 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<N, C::Type>> for Hat<C, N> {
    type Codomain = C::Type;

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

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

    #[inline]
    fn in_support(&self, x: &Loc<N, C::Type>) -> 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<N, C::Type>) -> [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) -> DynResult<C::Type> {
        Ok(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) -> DynResult<C::Type> {
        self.lipschitz_factor(L1)
            .map(|l1| <L2 as Dominated<C::Type, L1, Loc<N, C::Type>>>::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<N, C::Type>) -> 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<Linfinity, C::Type> for Hat<C, N> {
    #[inline]
    fn norm(&self, _: Linfinity) -> C::Type {
        self.bounds().upper()
    }
}