--- a/src/kernels/hat_convolution.rs Tue Aug 01 10:25:09 2023 +0300
+++ b/src/kernels/hat_convolution.rs Mon Feb 17 13:54:53 2025 -0500
@@ -14,9 +14,15 @@
GlobalAnalysis,
Bounded,
};
-use alg_tools::mapping::Apply;
+use alg_tools::mapping::{
+ Mapping,
+ Instance,
+ DifferentiableImpl,
+ Differential,
+};
use alg_tools::maputil::array_init;
+use crate::types::Lipschitz;
use super::base::*;
use super::ball_indicator::CubeIndicator;
@@ -38,6 +44,31 @@
/// -\frac{2}{3} (y-1)^3 & \frac{1}{2}\leq y<1. \\\\
/// \end{cases}
/// $$
+// Hence
+// $$
+// (h\*h)'(y) =
+// \begin{cases}
+// 2 (y+1)^2 & -1<y\leq -\frac{1}{2}, \\\\
+// -6 y^2-4 y & -\frac{1}{2}<y\leq 0, \\\\
+// 6 y^2-4 y & 0<y<\frac{1}{2}, \\\\
+// -2 (y-1)^2 & \frac{1}{2}\leq y<1. \\\\
+// \end{cases}
+// $$
+// as well as
+// $$
+// (h\*h)''(y) =
+// \begin{cases}
+// 4 (y+1) & -1<y\leq -\frac{1}{2}, \\\\
+// -12 y-4 & -\frac{1}{2}<y\leq 0, \\\\
+// 12 y-4 & 0<y<\frac{1}{2}, \\\\
+// -4 (y-1) & \frac{1}{2}\leq y<1. \\\\
+// \end{cases}
+// $$
+// This is maximised at y=±1/2 with value 2, and minimised at y=0 with value -4.
+// Now observe that
+// $$
+// [∇f(x\_1, …, x\_n)]_j = \frac{4}{σ} (h\*h)'(x\_j/σ) \prod\_{j ≠ i} \frac{4}{σ} (h\*h)(x\_i/σ)
+// $$
#[derive(Copy,Clone,Debug,Serialize,Eq)]
pub struct HatConv<S : Constant, const N : usize> {
/// The parameter $σ$ of the kernel.
@@ -60,24 +91,85 @@
}
}
-impl<'a, S, const N : usize> Apply<&'a Loc<S::Type, N>> for HatConv<S, N>
+impl<'a, S, const N : usize> Mapping<Loc<S::Type, N>> for HatConv<S, N>
where S : Constant {
- type Output = S::Type;
+ type Codomain = S::Type;
+
#[inline]
- fn apply(&self, y : &'a Loc<S::Type, N>) -> Self::Output {
+ fn apply<I : Instance<Loc<S::Type, N>>>(&self, y : I) -> Self::Codomain {
let σ = self.radius();
- y.product_map(|x| {
+ y.cow().product_map(|x| {
self.value_1d_σ1(x / σ) / σ
})
}
}
-impl<'a, S, const N : usize> Apply<Loc<S::Type, N>> for HatConv<S, N>
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<S, const N : usize> Lipschitz<L1> for HatConv<S, N>
+where S : Constant {
+ type FloatType = S::Type;
+ #[inline]
+ fn lipschitz_factor(&self, L1 : L1) -> Option<Self::FloatType> {
+ // For any ψ_i, we have
+ // ∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)
+ // = [ψ_1(x_1)-ψ_1(y_1)] ∏_{i=2}^N ψ_i(x_i)
+ // + ψ_1(y_1)[ ∏_{i=2}^N ψ_i(x_i) - ∏_{i=2}^N ψ_i(y_i)]
+ // = ∑_{j=1}^N [ψ_j(x_j)-ψ_j(y_j)]∏_{i > j} ψ_i(x_i) ∏_{i < j} ψ_i(y_i)
+ // Thus
+ // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)|
+ // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_j |ψ_j|
+ let σ = self.radius();
+ let l1d = self.lipschitz_1d_σ1() / (σ*σ);
+ let m1d = self.value_1d_σ1(0.0) / σ;
+ Some(l1d * m1d.powi(N as i32 - 1))
+ }
+}
+
+impl<S, const N : usize> Lipschitz<L2> for HatConv<S, N>
+where S : Constant {
+ type FloatType = S::Type;
+ #[inline]
+ fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
+ self.lipschitz_factor(L1).map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
+ }
+}
+
+
+impl<'a, S, const N : usize> DifferentiableImpl<Loc<S::Type, N>> for HatConv<S, N>
where S : Constant {
- type Output = S::Type;
+ type Derivative = Loc<S::Type, N>;
+
#[inline]
- fn apply(&self, y : Loc<S::Type, N>) -> Self::Output {
- self.apply(&y)
+ fn differential_impl<I : Instance<Loc<S::Type, N>>>(&self, y0 : I) -> Self::Derivative {
+ let y = y0.cow();
+ let σ = self.radius();
+ let σ2 = σ * σ;
+ let vs = y.map(|x| {
+ self.value_1d_σ1(x / σ) / σ
+ });
+ product_differential(&*y, &vs, |x| {
+ self.diff_1d_σ1(x / σ) / σ2
+ })
+ }
+}
+
+
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<'a, F : Float, S, const N : usize> Lipschitz<L2>
+for Differential<'a, Loc<F, N>, HatConv<S, N>>
+where S : Constant<Type=F> {
+ type FloatType = F;
+
+ #[inline]
+ fn lipschitz_factor(&self, _l2 : L2) -> Option<F> {
+ let h = self.base_fn();
+ let σ = h.radius();
+ Some(product_differential_lipschitz_factor::<F, N>(
+ h.value_1d_σ1(0.0) / σ,
+ h.lipschitz_1d_σ1() / (σ*σ),
+ h.maxabsdiff_1d_σ1() / (σ*σ),
+ h.lipschitz_diff_1d_σ1() / (σ*σ),
+ ))
}
}
@@ -97,6 +189,54 @@
(4.0/3.0) + 8.0 * y * y * (y - 1.0)
}
}
+
+ /// Computes the differential of the kernel for $n=1$ with $σ=1$.
+ #[inline]
+ fn diff_1d_σ1(&self, x : F) -> F {
+ let y = x.abs();
+ if y >= 1.0 {
+ 0.0
+ } else if y > 0.5 {
+ - 8.0 * (y - 1.0).powi(2)
+ } else /* 0 ≤ y ≤ 0.5 */ {
+ (24.0 * y - 16.0) * y
+ }
+ }
+
+ /// Computes the Lipschitz factor of the kernel for $n=1$ with $σ=1$.
+ #[inline]
+ fn lipschitz_1d_σ1(&self) -> F {
+ // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis
+ 2.0
+ }
+
+ /// Computes the maximum absolute differential of the kernel for $n=1$ with $σ=1$.
+ #[inline]
+ fn maxabsdiff_1d_σ1(&self) -> F {
+ // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis
+ 2.0
+ }
+
+ /// Computes the second differential of the kernel for $n=1$ with $σ=1$.
+ #[inline]
+ #[allow(dead_code)]
+ fn diff2_1d_σ1(&self, x : F) -> F {
+ let y = x.abs();
+ if y >= 1.0 {
+ 0.0
+ } else if y > 0.5 {
+ - 16.0 * (y - 1.0)
+ } else /* 0 ≤ y ≤ 0.5 */ {
+ 48.0 * y - 16.0
+ }
+ }
+
+ /// Computes the differential of the kernel for $n=1$ with $σ=1$.
+ #[inline]
+ fn lipschitz_diff_1d_σ1(&self) -> F {
+ // Maximal absolute second differential achieved at 0 by diff2_1d_σ1 analysis
+ 16.0
+ }
}
impl<'a, S, const N : usize> Support<S::Type, N> for HatConv<S, N>
@@ -159,21 +299,21 @@
}
#[replace_float_literals(F::cast_from(literal))]
-impl<'a, F : Float, R, C, const N : usize> Apply<&'a Loc<F, N>>
+impl<'a, F : Float, R, C, const N : usize> Mapping<Loc<F, N>>
for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
where R : Constant<Type=F>,
C : Constant<Type=F> {
- type Output = F;
+ type Codomain = F;
#[inline]
- fn apply(&self, y : &'a Loc<F, N>) -> F {
+ fn apply<I : Instance<Loc<F, N>>>(&self, y : I) -> F {
let Convolution(ref ind, ref hatconv) = self;
let β = ind.r.value();
let σ = hatconv.radius();
// This is just a product of one-dimensional versions
- y.product_map(|x| {
+ y.cow().product_map(|x| {
// With $u_σ(x) = u_1(x/σ)/σ$ the normalised hat convolution
// we have
// $$
@@ -188,24 +328,66 @@
}
}
-impl<'a, F : Float, R, C, const N : usize> Apply<Loc<F, N>>
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, const N : usize> DifferentiableImpl<Loc<F, N>>
for Convolution<CubeIndicator<R, N>, HatConv<C, N>>
where R : Constant<Type=F>,
C : Constant<Type=F> {
- type Output = F;
+ type Derivative = Loc<F, N>;
#[inline]
- fn apply(&self, y : Loc<F, N>) -> F {
- self.apply(&y)
+ fn differential_impl<I : Instance<Loc<F, N>>>(&self, y0 : I) -> Loc<F, N> {
+ let y = y0.cow();
+ let Convolution(ref ind, ref hatconv) = self;
+ let β = ind.r.value();
+ let σ = hatconv.radius();
+ let σ2 = σ * σ;
+
+ let vs = y.map(|x| {
+ self.value_1d_σ1(x / σ, β / σ)
+ });
+ product_differential(&*y, &vs, |x| {
+ self.diff_1d_σ1(x / σ, β / σ) / σ2
+ })
}
}
+/// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
+/// If $b > d$, add $g()$ to the result.
+#[inline]
+#[replace_float_literals(F::cast_from(literal))]
+fn i<F: Float>(a : F, b : F, c : F, d : F, f : impl Fn(F) -> F,
+ g : impl Fn() -> F) -> F {
+ if b < c {
+ 0.0
+ } else if b <= d {
+ if a <= c {
+ f(b) - f(c)
+ } else {
+ f(b) - f(a)
+ }
+ } else /* b > d */ {
+ g() + if a <= c {
+ f(d) - f(c)
+ } else if a < d {
+ f(d) - f(a)
+ } else {
+ 0.0
+ }
+ }
+}
+
#[replace_float_literals(F::cast_from(literal))]
impl<F : Float, C, R, const N : usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
where R : Constant<Type=F>,
C : Constant<Type=F> {
+
+ /// Calculates the value of the 1D hat convolution further convolved by a interval indicator.
+ /// As both functions are piecewise polynomials, this is implemented by explicit integral over
+ /// all subintervals of polynomiality of the cube indicator, using easily formed
+ /// antiderivatives.
#[inline]
pub fn value_1d_σ1(&self, x : F, β : F) -> F {
// The integration interval
@@ -218,34 +400,10 @@
y * y
}
- /// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
- /// If $b > d$, add $g()$ to the result.
- #[inline]
- fn i<F: Float>(a : F, b : F, c : F, d : F, f : impl Fn(F) -> F,
- g : impl Fn() -> F) -> F {
- if b < c {
- 0.0
- } else if b <= d {
- if a <= c {
- f(b) - f(c)
- } else {
- f(b) - f(a)
- }
- } else /* b > d */ {
- g() + if a <= c {
- f(d) - f(c)
- } else if a < d {
- f(d) - f(a)
- } else {
- 0.0
- }
- }
- }
-
// Observe the factor 1/6 at the front from the antiderivatives below.
// The factor 4 is from normalisation of the original function.
(4.0/6.0) * i(a, b, -1.0, -0.5,
- // (2/3) (y+1)^3 on -1 < y ≤ - 1/2
+ // (2/3) (y+1)^3 on -1 < y ≤ -1/2
// The antiderivative is (2/12)(y+1)^4 = (1/6)(y+1)^4
|y| pow4(y+1.0),
|| i(a, b, -0.5, 0.0,
@@ -266,8 +424,53 @@
)
)
}
+
+ /// Calculates the derivative of the 1D hat convolution further convolved by a interval
+ /// indicator. The implementation is similar to [`Self::value_1d_σ1`], using the fact that
+ /// $(θ * ψ)' = θ * ψ'$.
+ #[inline]
+ pub fn diff_1d_σ1(&self, x : F, β : F) -> F {
+ // The integration interval
+ let a = x - β;
+ let b = x + β;
+
+ // The factor 4 is from normalisation of the original function.
+ 4.0 * i(a, b, -1.0, -0.5,
+ // (2/3) (y+1)^3 on -1 < y ≤ -1/2
+ |y| (2.0/3.0) * (y + 1.0).powi(3),
+ || i(a, b, -0.5, 0.0,
+ // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0
+ |y| -2.0*(y + 1.0) * y * y + (1.0/3.0),
+ || i(a, b, 0.0, 0.5,
+ // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+ |y| 2.0*(y - 1.0) * y * y + (1.0/3.0),
+ || i(a, b, 0.5, 1.0,
+ // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
+ |y| -(2.0/3.0) * (y - 1.0).powi(3),
+ || 0.0
+ )
+ )
+ )
+ )
+ }
}
+/*
+impl<'a, F : Float, R, C, const N : usize> Lipschitz<L2>
+for Differential<Loc<F, N>, Convolution<CubeIndicator<R, N>, HatConv<C, N>>>
+where R : Constant<Type=F>,
+ C : Constant<Type=F> {
+
+ type FloatType = F;
+
+ #[inline]
+ fn lipschitz_factor(&self, _l2 : L2) -> Option<F> {
+ dbg!("unimplemented");
+ None
+ }
+}
+*/
+
impl<F : Float, R, C, const N : usize>
Convolution<CubeIndicator<R, N>, HatConv<C, N>>
where R : Constant<Type=F>,
@@ -409,7 +612,7 @@
#[cfg(test)]
mod tests {
use alg_tools::lingrid::linspace;
- use alg_tools::mapping::Apply;
+ use alg_tools::mapping::Mapping;
use alg_tools::norms::Linfinity;
use alg_tools::loc::Loc;
use crate::kernels::{BallIndicator, CubeIndicator, Convolution};