src/kernels/mollifier.rs

//! Implementation of the standard mollifier

use alg_tools::bisection_tree::{Bounds, Constant, GlobalAnalysis, LocalAnalysis, Support};
use alg_tools::euclidean::Euclidean;
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 float_extras::f64::tgamma as gamma;
use numeric_literals::replace_float_literals;
use rgsl::hypergeometric::hyperg_U;
use serde::Serialize;

/// Reresentation of the (unnormalised) standard mollifier.
///
/// For the `width` parameter $ε>0$, this is
/// <div>$$
///     f(x)=\begin{cases}
///         e^{\frac{ε^2}{\|x\|_2^2-ε^2}}, & \|x\|_2 < ε, \\
///         0, & \text{otherwise}.
///     \end{cases}
/// $$</div>
#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
pub struct Mollifier<C: Constant, const N: usize> {
    /// The parameter $ε$ of the mollifier.
    pub width: C,
}

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

    #[inline]
    fn apply<I: Instance<Loc<N, C::Type>>>(&self, x: I) -> Self::Codomain {
        let ε = self.width.value();
        let ε2 = ε * ε;
        let n2 = x.eval(|x| x.norm2_squared());
        if n2 < ε2 {
            (n2 / (n2 - ε2)).exp()
        } else {
            0.0
        }
    }
}

impl<'a, C: Constant, const N: usize> Support<N, C::Type> for Mollifier<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.norm2() < self.width.value()
    }

    /*fn fully_in_support(&self, _cube : &Cube<C::Type,N>) -> bool {
        todo!("Not implemented, but not used at the moment")
    }*/
}

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, C: Constant, const N: usize> GlobalAnalysis<C::Type, Bounds<C::Type>> for Mollifier<C, N> {
    #[inline]
    fn global_analysis(&self) -> Bounds<C::Type> {
        // The function is maximised/minimised where the 2-norm is minimised/maximised.
        Bounds(0.0, 1.0)
    }
}

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

/// Calculate integral of the standard mollifier of width 1 in $ℝ^n$.
///
/// This is based on the formula from
/// [https://math.stackexchange.com/questions/4359683/integral-of-the-usual-mollifier-function-finding-its-necessary-constant]().
///
/// If `rescaled` is `true`, return the integral of the scaled mollifier that has value one at the
/// origin.
#[inline]
pub fn mollifier_norm1(n_: usize, rescaled: bool) -> f64 {
    assert!(n_ > 0);
    let n = n_ as f64;
    let q = 2.0;
    let p = 2.0;
    let base = (2.0*gamma(1.0 + 1.0/p)).powi(n_ as i32)
               /*/ gamma(1.0 + n / p)
               * gamma(1.0 + n / q)*/
               * hyperg_U(1.0 + n / q, 2.0, 1.0);
    if rescaled {
        base
    } else {
        base / f64::E
    }
}

impl<'a, C: Constant, const N: usize> Norm<L1, C::Type> for Mollifier<C, N> {
    #[inline]
    fn norm(&self, _: L1) -> C::Type {
        let ε = self.width.value();
        C::Type::cast_from(mollifier_norm1(N, true)) * ε.powi(N as i32)
    }
}

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