pointsource_algs: comparison src/kernels/hat

-:6105b5cd8d89
+:f0e8704d3f0e
 Bounds,
 LocalAnalysis,
 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;
 /// Hat convolution kernel.
 ///
 ///         -2 y^3-2 y^2+\frac{1}{3} & -\frac{1}{2}<y\leq 0, \\\\
 ///         2 y^3-2 y^2+\frac{1}{3} & 0<y<\frac{1}{2}, \\\\
 ///         -\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.
 pub radius : S,
 }
 pub fn radius(&self) -> S::Type {
 self.radius.value()
 }
 }
-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 {
+#[inline]
+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 Output = S::Type;
+type FloatType = S::Type;
 #[inline]
-fn apply(&self, y : Loc<S::Type, N>) -> Self::Output {
+fn lipschitz_factor(&self, L1 : L1) -> Option<Self::FloatType> {
-self.apply(&y)
+// 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 Derivative = Loc<S::Type, N>;
+#[inline]
+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() / (σ*σ),
+))
 }
 }
 #[replace_float_literals(S::Type::cast_from(literal))]
 - (8.0/3.0) * (y - 1.0).powi(3)
 } else /* 0 ≤ y ≤ 0.5 */ {
 (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>
 where S : Constant {
 #[inline]
 self.bounds().upper()
 }
 }
 #[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
 // $$
 //      [χ_{-β,β} * u_σ](x)
 //      = ∫_{x-β}^{x+β} u_σ(z) d z
 self.value_1d_σ1(x / σ, β / σ)
 })
 }
 }
-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 {
+fn differential_impl<I : Instance<Loc<F, N>>>(&self, y0 : I) -> Loc<F, N> {
-self.apply(&y)
+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
 let a = x - β;
 let b = x + β;
 fn pow4<F : Float>(x : F) -> F {
 let y = x * x;
 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,
 // -2 y^3 - 2 y^2 + 1/3  on  -1/2 < y ≤ 0
 // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y
 )
 )
 )
 )
 }
-}
+/// 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>,
 C : Constant<Type=F> {
 }
 #[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};
 use super::HatConv;

comparison: src/kernels/hat_convolution.rs

src/kernels/hat_convolution.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/hat_convolution.rs

src/kernels/hat_convolution.rs