src/kernels/ball_indicator.rs

//! Implementation of the indicator function of a ball with respect to various norms.
use super::base::*;
use crate::types::*;
use alg_tools::bisection_tree::{Bounds, Constant, GlobalAnalysis, LocalAnalysis, Support};
use alg_tools::coefficients::factorial;
use alg_tools::euclidean::StaticEuclidean;
use alg_tools::instance::Instance;
use alg_tools::loc::Loc;
use alg_tools::mapping::{DifferentiableImpl, Differential, LipschitzDifferentiableImpl, Mapping};
use alg_tools::maputil::array_init;
use alg_tools::norms::*;
use alg_tools::sets::Cube;
use anyhow::anyhow;
use float_extras::f64::tgamma as gamma;
use numeric_literals::replace_float_literals;
use serde::Serialize;

/// Representation of the indicator of the ball $𝔹_q = \\{ x ∈ ℝ^N \mid \\|x\\|\_q ≤ r \\}$,
/// where $q$ is the `Exponent`, and $r$ is the radius [`Constant`] `C`.
#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
pub struct BallIndicator<C: Constant, Exponent: NormExponent, const N: usize> {
    /// The radius of the ball.
    pub r: C,
    /// The exponent $q$ of the norm creating the ball
    pub exponent: Exponent,
}

/// Alias for the representation of the indicator of the $∞$-norm-ball
/// $𝔹_∞ = \\{ x ∈ ℝ^N \mid \\|x\\|\_∞ ≤ c \\}$.
pub type CubeIndicator<C, const N: usize> = BallIndicator<C, Linfinity, N>;

#[replace_float_literals(C::Type::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
    Mapping<Loc<N, C::Type>> for BallIndicator<C, Exponent, N>
where
    Loc<N, F>: Norm<Exponent, F>,
{
    type Codomain = C::Type;

    #[inline]
    fn apply<I: Instance<Loc<N, C::Type>>>(&self, x: I) -> Self::Codomain {
        let r = self.r.value();
        let n = x.eval(|x| x.norm(self.exponent));
        if n <= r {
            1.0
        } else {
            0.0
        }
    }
}

impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
    DifferentiableImpl<Loc<N, C::Type>> for BallIndicator<C, Exponent, N>
where
    C: Constant,
    Loc<N, F>: Norm<Exponent, F>,
{
    type Derivative = Loc<N, C::Type>;

    #[inline]
    fn differential_impl<I: Instance<Loc<N, C::Type>>>(&self, _x: I) -> Self::Derivative {
        Self::Derivative::origin()
    }
}

impl<F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> Lipschitz<L2>
    for BallIndicator<C, Exponent, N>
where
    C: Constant,
    Loc<N, F>: Norm<Exponent, F>,
{
    type FloatType = C::Type;

    fn lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
        Err(anyhow!("Not a Lipschitz function"))
    }
}

impl<'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
    LipschitzDifferentiableImpl<Loc<N, F>, L2> for BallIndicator<C, Exponent, N>
where
    C: Constant,
    Loc<N, F>: Norm<Exponent, F>,
{
    type FloatType = C::Type;

    fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
        Err(anyhow!("Not a Lipschitz-differentiable function"))
    }
}

impl<'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> NormBounded<L2>
    for Differential<'b, Loc<N, F>, BallIndicator<C, Exponent, N>>
where
    C: Constant,
    Loc<N, F>: Norm<Exponent, F>,
{
    type FloatType = C::Type;

    fn norm_bound(&self, _l2: L2) -> C::Type {
        F::INFINITY
    }
}

impl<'a, 'b, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
    NormBounded<L2> for Differential<'b, Loc<N, F>, &'a BallIndicator<C, Exponent, N>>
where
    C: Constant,
    Loc<N, F>: Norm<Exponent, F>,
{
    type FloatType = C::Type;

    fn norm_bound(&self, _l2: L2) -> C::Type {
        F::INFINITY
    }
}

impl<'a, F: Float, C: Constant<Type = F>, Exponent, const N: usize> Support<N, C::Type>
    for BallIndicator<C, Exponent, N>
where
    Exponent: NormExponent + Sync + Send + 'static,
    Loc<N, F>: Norm<Exponent, F>,
    Linfinity: Dominated<F, Exponent, Loc<N, F>>,
{
    #[inline]
    fn support_hint(&self) -> Cube<N, F> {
        let r = Linfinity.from_norm(self.r.value(), self.exponent);
        array_init(|| [-r, r]).into()
    }

    #[inline]
    fn in_support(&self, x: &Loc<N, F>) -> bool {
        let r = Linfinity.from_norm(self.r.value(), self.exponent);
        x.norm(self.exponent) <= r
    }

    /// This can only really work in a reasonable fashion for N=1.
    #[inline]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        let r = Linfinity.from_norm(self.r.value(), self.exponent);
        cube.map(|a, b| symmetric_interval_hint(r, a, b))
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize>
    GlobalAnalysis<F, Bounds<F>> for BallIndicator<C, Exponent, N>
where
    Loc<N, F>: Norm<Exponent, F>,
{
    #[inline]
    fn global_analysis(&self) -> Bounds<F> {
        Bounds(0.0, 1.0)
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, Exponent: NormExponent, const N: usize> Norm<Linfinity, F>
    for BallIndicator<C, Exponent, N>
where
    Loc<N, F>: Norm<Exponent, F>,
{
    #[inline]
    fn norm(&self, _: Linfinity) -> F {
        1.0
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F> for BallIndicator<C, L1, N> {
    #[inline]
    fn norm(&self, _: L1) -> F {
        // Using https://en.wikipedia.org/wiki/Volume_of_an_n-ball#Balls_in_Lp_norms,
        // we have V_N^1(r) = (2r)^N / N!
        let r = self.r.value();
        if N == 1 {
            2.0 * r
        } else if N == 2 {
            r * r
        } else {
            (2.0 * r).powi(N as i32) * F::cast_from(factorial(N))
        }
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F> for BallIndicator<C, L2, N> {
    #[inline]
    fn norm(&self, _: L1) -> F {
        // See https://en.wikipedia.org/wiki/Volume_of_an_n-ball#The_volume.
        let r = self.r.value();
        let π = F::PI;
        if N == 1 {
            2.0 * r
        } else if N == 2 {
            π * (r * r)
        } else {
            let ndiv2 = F::cast_from(N) / 2.0;
            let γ = F::cast_from(gamma((ndiv2 + 1.0).as_()));
            π.powf(ndiv2) / γ * r.powi(N as i32)
        }
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, C: Constant<Type = F>, const N: usize> Norm<L1, F>
    for BallIndicator<C, Linfinity, N>
{
    #[inline]
    fn norm(&self, _: L1) -> F {
        let two_r = 2.0 * self.r.value();
        two_r.powi(N as i32)
    }
}

macro_rules! indicator_local_analysis {
    ($exponent:ident) => {
        impl<'a, F: Float, C: Constant<Type = F>, const N: usize> LocalAnalysis<F, Bounds<F>, N>
            for BallIndicator<C, $exponent, N>
        where
            Loc<N, F>: Norm<$exponent, F>,
            Linfinity: Dominated<F, $exponent, Loc<N, F>>,
        {
            #[inline]
            fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
                // 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)
            }
        }
    };
}

indicator_local_analysis!(L1);
indicator_local_analysis!(L2);
indicator_local_analysis!(Linfinity);

#[replace_float_literals(F::cast_from(literal))]
impl<'a, F: Float, R, const N: usize> Mapping<Loc<N, F>> for AutoConvolution<CubeIndicator<R, N>>
where
    R: Constant<Type = F>,
{
    type Codomain = F;

    #[inline]
    fn apply<I: Instance<Loc<N, F>>>(&self, y: I) -> F {
        let two_r = 2.0 * self.0.r.value();
        // This is just a product of one-dimensional versions
        y.decompose()
            .iter()
            .map(|&x| 0.0.max(two_r - x.abs()))
            .product()
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<F: Float, R, const N: usize> Support<N, F> for AutoConvolution<CubeIndicator<R, N>>
where
    R: Constant<Type = F>,
{
    #[inline]
    fn support_hint(&self) -> Cube<N, F> {
        let two_r = 2.0 * self.0.r.value();
        array_init(|| [-two_r, two_r]).into()
    }

    #[inline]
    fn in_support(&self, y: &Loc<N, F>) -> bool {
        let two_r = 2.0 * self.0.r.value();
        y.iter().all(|x| x.abs() <= two_r)
    }

    #[inline]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        let two_r = 2.0 * self.0.r.value();
        cube.map(|c, d| symmetric_interval_hint(two_r, c, d))
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<F: Float, R, const N: usize> GlobalAnalysis<F, Bounds<F>>
    for AutoConvolution<CubeIndicator<R, N>>
where
    R: Constant<Type = F>,
{
    #[inline]
    fn global_analysis(&self) -> Bounds<F> {
        Bounds(0.0, self.apply(Loc::ORIGIN))
    }
}

impl<F: Float, R, const N: usize> LocalAnalysis<F, Bounds<F>, N>
    for AutoConvolution<CubeIndicator<R, N>>
where
    R: Constant<Type = F>,
{
    #[inline]
    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
        // The function is maximised/minimised where the absolute value is minimised/maximised.
        let lower = self.apply(cube.maxnorm_point());
        let upper = self.apply(cube.minnorm_point());
        Bounds(lower, upper)
    }
}