src/bisection_tree/support.rs

/*!
Traits for representing the support of a [`Mapping`], and analysing the mapping on a [`Cube`].
*/
use super::aggregator::Bounds;
pub use crate::bounds::{GlobalAnalysis, LocalAnalysis};
use crate::loc::Loc;
use crate::mapping::{DifferentiableImpl, DifferentiableMapping, Instance, Mapping, Space};
use crate::maputil::map2;
use crate::norms::{Linfinity, Norm, L1, L2};
pub use crate::operator_arithmetic::{Constant, Weighted};
use crate::sets::Cube;
use crate::types::{Float, Num};
use serde::Serialize;
use std::ops::{DivAssign, MulAssign, Neg};

/// A trait for working with the supports of [`Mapping`]s.
///
/// `Mapping` is not a super-trait to allow more general use.
pub trait Support<const N: usize, F: Num>: Sized + Sync + Send + 'static {
    /// Return a cube containing the support of the function represented by `self`.
    ///
    /// The hint may be larger than the actual support, but must contain it.
    fn support_hint(&self) -> Cube<N, F>;

    /// Indicate whether `x` is in the support of the function represented by `self`.
    fn in_support(&self, x: &Loc<N, F>) -> bool;

    // Indicate whether `cube` is fully in the support of the function represented by `self`.
    //fn fully_in_support(&self, cube : &Cube<F,N>) -> bool;

    /// Return an optional hint for bisecting the support.
    ///
    /// The output along each axis a possible coordinate at which to bisect `cube`.
    ///
    /// This is useful for nonsmooth functions to make finite element models as used by
    /// [`BTFN`][super::btfn::BTFN] minimisation/maximisation compatible with points of
    /// non-differentiability.
    ///
    /// The default implementation returns `[None; N]`.
    #[inline]
    #[allow(unused_variables)]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        [None; N]
    }

    /// Translate `self` by `x`.
    #[inline]
    fn shift(self, x: Loc<N, F>) -> Shift<Self, N, F> {
        Shift {
            shift: x,
            base_fn: self,
        }
    }
}

/// Shift of [`Support`] and [`Mapping`]; output of [`Support::shift`].
#[derive(Copy, Clone, Debug, Serialize)] // Serialize! but not implemented by Loc.
pub struct Shift<T, const N: usize, F = f64> {
    shift: Loc<N, F>,
    base_fn: T,
}

impl<'a, T, V: Space, F: Float, const N: usize> Mapping<Loc<N, F>> for Shift<T, N, F>
where
    T: Mapping<Loc<N, F>, Codomain = V>,
{
    type Codomain = V;

    #[inline]
    fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
        self.base_fn.apply(x.own() - &self.shift)
    }
}

impl<'a, T, V: Space, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for Shift<T, N, F>
where
    T: DifferentiableMapping<Loc<N, F>, DerivativeDomain = V>,
{
    type Derivative = V;

    #[inline]
    fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative {
        self.base_fn.differential(x.own() - &self.shift)
    }
}

impl<'a, T, F: Float, const N: usize> Support<N, F> for Shift<T, N, F>
where
    T: Support<N, F>,
{
    #[inline]
    fn support_hint(&self) -> Cube<N, F> {
        self.base_fn.support_hint().shift(&self.shift)
    }

    #[inline]
    fn in_support(&self, x: &Loc<N, F>) -> bool {
        self.base_fn.in_support(&(x - &self.shift))
    }

    // fn fully_in_support(&self, _cube : &Cube<F,N>) -> bool {
    //     //self.base_fn.fully_in_support(cube.shift(&vectorneg(self.shift)))
    //     todo!("Not implemented, but not used at the moment")
    // }

    #[inline]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        let base_hint = self.base_fn.bisection_hint(cube);
        map2(base_hint, &self.shift, |h, s| h.map(|z| z + *s))
    }
}

impl<'a, T, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for Shift<T, N, F>
where
    T: LocalAnalysis<F, Bounds<F>, N>,
{
    #[inline]
    fn global_analysis(&self) -> Bounds<F> {
        self.base_fn.global_analysis()
    }
}

impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Shift<T, N, F>
where
    T: LocalAnalysis<F, Bounds<F>, N>,
{
    #[inline]
    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
        self.base_fn.local_analysis(&cube.shift(&(-self.shift)))
    }
}

macro_rules! impl_shift_norm {
    ($($norm:ident)*) => { $(
        impl<'a, T, F : Float, const N : usize> Norm<$norm, F> for Shift<T, N, F>
        where T : Norm<$norm, F> {
            #[inline]
            fn norm(&self, n : $norm) -> F {
                self.base_fn.norm(n)
            }
        }
    )* }
}

impl_shift_norm!(L1 L2 Linfinity);

impl<'a, T, F: Float, C, const N: usize> Support<N, F> for Weighted<T, C>
where
    T: Support<N, F>,
    C: Constant<Type = F>,
{
    #[inline]
    fn support_hint(&self) -> Cube<N, F> {
        self.base_fn.support_hint()
    }

    #[inline]
    fn in_support(&self, x: &Loc<N, F>) -> bool {
        self.base_fn.in_support(x)
    }

    // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool {
    //     self.base_fn.fully_in_support(cube)
    // }

    #[inline]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        self.base_fn.bisection_hint(cube)
    }
}

impl<'a, T, F: Float, C> GlobalAnalysis<F, Bounds<F>> for Weighted<T, C>
where
    T: GlobalAnalysis<F, Bounds<F>>,
    C: Constant<Type = F>,
{
    #[inline]
    fn global_analysis(&self) -> Bounds<F> {
        let Bounds(lower, upper) = self.base_fn.global_analysis();
        debug_assert!(lower <= upper);
        match self.weight.value() {
            w if w < F::ZERO => Bounds(w * upper, w * lower),
            w => Bounds(w * lower, w * upper),
        }
    }
}

impl<'a, T, F: Float, C, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Weighted<T, C>
where
    T: LocalAnalysis<F, Bounds<F>, N>,
    C: Constant<Type = F>,
{
    #[inline]
    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
        let Bounds(lower, upper) = self.base_fn.local_analysis(cube);
        debug_assert!(lower <= upper);
        match self.weight.value() {
            w if w < F::ZERO => Bounds(w * upper, w * lower),
            w => Bounds(w * lower, w * upper),
        }
    }
}

macro_rules! make_weighted_scalarop_rhs {
    ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
        impl<F: Float, T> std::ops::$trait_assign<F> for Weighted<T, F> {
            #[inline]
            fn $fn_assign(&mut self, t: F) {
                self.weight.$fn_assign(t);
            }
        }

        impl<'a, F: Float, T> std::ops::$trait<F> for Weighted<T, F> {
            type Output = Self;
            #[inline]
            fn $fn(mut self, t: F) -> Self {
                self.weight.$fn_assign(t);
                self
            }
        }

        impl<'a, F: Float, T> std::ops::$trait<F> for &'a Weighted<T, F>
        where
            T: Clone,
        {
            type Output = Weighted<T, F>;
            #[inline]
            fn $fn(self, t: F) -> Self::Output {
                Weighted {
                    weight: self.weight.$fn(t),
                    base_fn: self.base_fn.clone(),
                }
            }
        }
    };
}

make_weighted_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
make_weighted_scalarop_rhs!(Div, div, DivAssign, div_assign);

macro_rules! impl_weighted_norm {
    ($($norm:ident)*) => { $(
        impl<'a, T, F : Float> Norm<$norm, F> for Weighted<T,F>
        where T : Norm<$norm, F> {
            #[inline]
            fn norm(&self, n : $norm) -> F {
                self.base_fn.norm(n) * self.weight.abs()
            }
        }
    )* }
}

impl_weighted_norm!(L1 L2 Linfinity);

/// Normalisation of [`Support`] and [`Mapping`] to L¹ norm 1.
///
/// Currently only scalar-valued functions are supported.
#[derive(Copy, Clone, Debug, Serialize, PartialEq)]
pub struct Normalised<T>(
    /// The base [`Support`] or [`Mapping`].
    pub T,
);

impl<'a, T, F: Float, const N: usize> Mapping<Loc<N, F>> for Normalised<T>
where
    T: Norm<L1, F> + Mapping<Loc<N, F>, Codomain = F>,
{
    type Codomain = F;

    #[inline]
    fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
        let w = self.0.norm(L1);
        if w == F::ZERO {
            F::ZERO
        } else {
            self.0.apply(x) / w
        }
    }
}

impl<'a, T, F: Float, const N: usize> Support<N, F> for Normalised<T>
where
    T: Norm<L1, F> + Support<N, F>,
{
    #[inline]
    fn support_hint(&self) -> Cube<N, F> {
        self.0.support_hint()
    }

    #[inline]
    fn in_support(&self, x: &Loc<N, F>) -> bool {
        self.0.in_support(x)
    }

    // fn fully_in_support(&self, cube : &Cube<F,N>) -> bool {
    //     self.0.fully_in_support(cube)
    // }

    #[inline]
    fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
        self.0.bisection_hint(cube)
    }
}

impl<'a, T, F: Float> GlobalAnalysis<F, Bounds<F>> for Normalised<T>
where
    T: Norm<L1, F> + GlobalAnalysis<F, Bounds<F>>,
{
    #[inline]
    fn global_analysis(&self) -> Bounds<F> {
        let Bounds(lower, upper) = self.0.global_analysis();
        debug_assert!(lower <= upper);
        let w = self.0.norm(L1);
        debug_assert!(w >= F::ZERO);
        Bounds(w * lower, w * upper)
    }
}

impl<'a, T, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for Normalised<T>
where
    T: Norm<L1, F> + LocalAnalysis<F, Bounds<F>, N>,
{
    #[inline]
    fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
        let Bounds(lower, upper) = self.0.local_analysis(cube);
        debug_assert!(lower <= upper);
        let w = self.0.norm(L1);
        debug_assert!(w >= F::ZERO);
        Bounds(w * lower, w * upper)
    }
}

impl<'a, T, F: Float> Norm<L1, F> for Normalised<T>
where
    T: Norm<L1, F>,
{
    #[inline]
    fn norm(&self, _: L1) -> F {
        let w = self.0.norm(L1);
        if w == F::ZERO {
            F::ZERO
        } else {
            F::ONE
        }
    }
}

macro_rules! impl_normalised_norm {
    ($($norm:ident)*) => { $(
        impl<'a, T, F : Float> Norm<$norm, F> for Normalised<T>
        where T : Norm<$norm, F> + Norm<L1, F> {
            #[inline]
            fn norm(&self, n : $norm) -> F {
                let w = self.0.norm(L1);
                if w == F::ZERO { F::ZERO } else { self.0.norm(n) / w }
            }
        }
    )* }
}

impl_normalised_norm!(L2 Linfinity);

/*
impl<F : Num, S : Support< N, F>, const N : usize> LocalAnalysis<F, NullAggregator, N> for S {
    fn local_analysis(&self, _cube : &Cube<N, F>) -> NullAggregator { NullAggregator }
}

impl<F : Float, S : Bounded<F>, const N : usize> LocalAnalysis<F, Bounds<F>, N> for S {
    #[inline]
    fn local_analysis(&self, cube : &Cube<N, F>) -> Bounds<F> {
        self.bounds(cube)
    }
}*/

/// Generator of [`Support`]-implementing component functions based on low storage requirement
/// [ids][`Self::Id`].
pub trait SupportGenerator<const N: usize, F: Float = f64>:
    MulAssign<F> + DivAssign<F> + Neg<Output = Self> + Clone + Sync + Send + 'static
{
    /// The identification type
    type Id: 'static + Copy;
    /// The type of the [`Support`] (often also a [`Mapping`]).
    type SupportType: 'static + Support<N, F>;
    /// An iterator over all the [`Support`]s of the generator.
    type AllDataIter<'a>: Iterator<Item = (Self::Id, Self::SupportType)>
    where
        Self: 'a;

    /// Returns the component identified by `id`.
    ///
    /// Panics if `id` is an invalid identifier.
    fn support_for(&self, id: Self::Id) -> Self::SupportType;

    /// Returns the number of different components in this generator.
    fn support_count(&self) -> usize;

    /// Returns an iterator over all pairs of `(id, support)`.
    fn all_data(&self) -> Self::AllDataIter<'_>;
}