--- a/src/kernels/ball_indicator.rs Sun Apr 27 15:03:51 2025 -0500
+++ b/src/kernels/ball_indicator.rs Thu Feb 26 11:38:43 2026 -0500
@@ -1,56 +1,44 @@
-
//! 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;
-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,
-};
-use alg_tools::mapping::{
- Mapping,
- Differential,
- DifferentiableImpl,
-};
-use alg_tools::instance::Instance;
-use alg_tools::euclidean::StaticEuclidean;
-use alg_tools::maputil::array_init;
-use alg_tools::coefficients::factorial;
-use crate::types::*;
-use super::base::*;
/// 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> {
+#[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,
+ pub r: C,
/// The exponent $q$ of the norm creating the ball
- pub exponent : Exponent,
+ 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>;
+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<C::Type, N>>
-for BallIndicator<C, Exponent, N>
+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<F, N> : Norm<F, Exponent>
+ Loc<N, F>: Norm<Exponent, F>,
{
type Codomain = C::Type;
#[inline]
- fn apply<I : Instance<Loc<C::Type, N>>>(&self, x : I) -> Self::Codomain {
+ 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 {
@@ -61,114 +49,105 @@
}
}
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-DifferentiableImpl<Loc<C::Type, N>>
-for BallIndicator<C, Exponent, N>
+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<F, N> : Norm<F, Exponent>
+ C: Constant,
+ Loc<N, F>: Norm<Exponent, F>,
{
- type Derivative = Loc<C::Type, N>;
+ type Derivative = Loc<N, C::Type>;
#[inline]
- fn differential_impl<I : Instance<Loc<C::Type, N>>>(&self, _x : I) -> Self::Derivative {
+ 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<F, N> : Norm<F, Exponent> {
+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) -> Option<C::Type> {
- None
+ 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>
-Lipschitz<L2>
-for Differential<'b, Loc<F, N>, BallIndicator<C, Exponent, N>>
-where C : Constant,
- Loc<F, N> : Norm<F, Exponent> {
+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 lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
- None
+ fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<C::Type> {
+ Err(anyhow!("Not a Lipschitz-differentiable function"))
}
}
-impl<'a, 'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Lipschitz<L2>
-for Differential<'b, Loc<F, N>, &'a BallIndicator<C, Exponent, N>>
-where C : Constant,
- Loc<F, N> : Norm<F, Exponent> {
+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 lipschitz_factor(&self, _l2 : L2) -> Option<C::Type> {
- None
- }
-}
-
-
-impl<'b, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-NormBounded<L2>
-for Differential<'b, Loc<F, N>, BallIndicator<C, Exponent, N>>
-where C : Constant,
- Loc<F, N> : Norm<F, Exponent> {
- type FloatType = C::Type;
-
- fn norm_bound(&self, _l2 : L2) -> 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<F, N>, &'a BallIndicator<C, Exponent, N>>
-where C : Constant,
- Loc<F, N> : Norm<F, Exponent> {
+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 {
+ fn norm_bound(&self, _l2: L2) -> C::Type {
F::INFINITY
}
}
-
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Support<C::Type, N>
-for BallIndicator<C, Exponent, N>
-where Loc<F, N> : Norm<F, Exponent>,
- Linfinity : Dominated<F, Exponent, Loc<F, N>> {
-
+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<F,N> {
+ 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<F,N>) -> bool {
+ 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<F, N>) -> [Option<F>; N] {
+ 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<F, N> : Norm<F, Exponent> {
+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)
@@ -176,29 +155,28 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, Exponent : NormExponent, const N : usize>
-Norm<F, Linfinity>
-for BallIndicator<C, Exponent, N>
-where Loc<F, N> : Norm<F, Exponent> {
+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 {
+ 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<F, L1>
-for BallIndicator<C, L1, N> {
+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 {
+ 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 {
+ if N == 1 {
2.0 * r
- } else if N==2 {
- r*r
+ } else if N == 2 {
+ r * r
} else {
(2.0 * r).powi(N as i32) * F::cast_from(factorial(N))
}
@@ -206,17 +184,15 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-Norm<F, L1>
-for BallIndicator<C, L2, N> {
+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 {
+ 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 {
+ if N == 1 {
2.0 * r
- } else if N==2 {
+ } else if N == 2 {
π * (r * r)
} else {
let ndiv2 = F::cast_from(N) / 2.0;
@@ -227,94 +203,100 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, C : Constant<Type=F>, const N : usize>
-Norm<F, L1>
-for BallIndicator<C, Linfinity, N> {
+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 {
+ 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<F, N> : Norm<F, $exponent>,
- Linfinity : Dominated<F, $exponent, Loc<F, N>> {
+ 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<F, N>) -> Bounds<F> {
+ 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<F, N>>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+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<F, N>>>(&self, y : I) -> F {
+ 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.cow().iter().map(|&x| {
- 0.0.max(two_r - x.abs())
- }).product()
+ 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<F, N>
-for AutoConvolution<CubeIndicator<R, N>>
-where R : Constant<Type=F> {
+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<F, N> {
+ 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<F, N>) -> bool {
+ 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<F, N>) -> [Option<F>; N] {
+ 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> {
+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> {
+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<F, N>) -> Bounds<F> {
+ 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());