--- a/src/kernels/hat_convolution.rs Fri May 15 14:40:02 2026 -0500
+++ b/src/kernels/hat_convolution.rs Sun Jul 19 07:34:39 2026 +0200
@@ -97,42 +97,30 @@
#[inline]
fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain {
let σ = self.radius();
- y.decompose().product_map(|x| self.value_1d_σ1(x / σ) / σ)
+ y.decompose().product_map(|x| self.value_1d_σ1(x / σ)) / σ.powi(N as i32)
}
}
#[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) -> DynResult<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) / σ;
- Ok(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) -> DynResult<Self::FloatType> {
- self.lipschitz_factor(L1)
- .map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
+ fn lipschitz_factor(&self, _L2: L2) -> DynResult<Self::FloatType> {
+ // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
+ // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j < i} ψ_j(y_j) ∏_{j > i} ψ_j(x_j)
+ // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have
+ // |ψ(x) - ψ(y)| ≤ ∑_i L_g|x_i-y_i| M_g^{n-1} ≤ √n L_g M_g^{n-1} |x-y|₂
+ // where L_g is the Lipschitz factor of g, and M_g a bound on it.
+ let σ = self.radius();
+ let l1d = self.lipschitz_1d_σ1();
+ let m1d = self.value_1d_σ1(0.0);
+ Ok(
+ S::Type::cast_from(N).sqrt() * l1d * m1d.powi((N - 1) as i32)
+ / σ.powi((N - 1 + 2) as i32),
+ )
}
}
@@ -146,9 +134,9 @@
fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative {
let y = y0.decompose();
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)
+ let tmp = y.map(|x| x / σ);
+ let vs = tmp.map(|x| self.value_1d_σ1(x));
+ product_differential(&tmp, &vs, |x| self.diff_1d_σ1(x)) / σ.powi((N + 1) as i32)
}
}
@@ -164,11 +152,14 @@
fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> {
let σ = self.radius();
product_differential_lipschitz_factor::<F, N>(
- self.value_1d_σ1(0.0) / σ,
- self.lipschitz_1d_σ1() / (σ * σ),
- self.maxabsdiff_1d_σ1() / (σ * σ),
- self.lipschitz_diff_1d_σ1() / (σ * σ),
+ self.value_1d_σ1(0.0),
+ self.lipschitz_1d_σ1(),
+ self.maxabsdiff_1d_σ1(),
+ self.lipschitz_diff_1d_σ1(),
)
+ .map(|l| l / σ.powi((N + 2) as i32))
+ // We have f(x)=f₀(x/σ)/σ^N, where f₀(z)=∏_{i=1}^N f₁(z_i).
+ // Thus |f'(x)-f'(y)| = |f₀'(x/σ)-f₀'(y/σ)|/σ^{N+1} ≤ L_{f₀'}/σ^{N+2}.
}
}
@@ -195,15 +186,15 @@
/// 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 {
+ if x >= 1.0 || x <= -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
+ } else if x > 0.5 {
+ -8.0 * (x - 1.0).powi(2)
+ } else if x < -0.5 {
+ 8.0 * (x + 1.0).powi(2)
+ } else {
+ /* 0 ≤ y ≤ 0.5 */
+ (24.0 * x.abs() - 16.0) * x
}
}
@@ -221,22 +212,6 @@
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 {
@@ -332,7 +307,7 @@
// = ∫_{x-β}^{x+β} u_σ(z) d z
// = (1/σ)∫_{x-β}^{x+β} u_1(z/σ) d z
// = ∫_{(x-β)/σ}^{(x+β)/σ} u_1(z) d z
- // = [χ_{-β/σ, β/σ} * u_1](x/σ)
+ // = (χ_{[-β/σ, β/σ]} * u_1)(x/σ)
// $$
self.value_1d_σ1(x / σ, β / σ)
})
@@ -354,14 +329,14 @@
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)
+ product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ)) / σ
}
}
-/// Integrate $f$, whose support is $[c, d]$, on $[a, b]$.
+/// Integrate $f'$, whose support is $[c, d]$, on $[a, b]$.
+/// The value $f$ is given, not the derivative.
/// If $b > d$, add $g()$ to the result.
#[inline]
#[replace_float_literals(F::cast_from(literal))]
@@ -403,12 +378,6 @@
let a = x - β;
let b = x + β;
- #[inline]
- fn pow4<F: Float>(x: F) -> F {
- let y = x * x;
- y * y
- }
-
// 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)
@@ -419,37 +388,28 @@
-0.5,
// (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),
+ |y| (y + 1.0).powi(4),
|| {
i(
a,
b,
-0.5,
- 0.0,
+ 0.5,
// -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
- |y| y * (-y * y * (y * 3.0 + 4.0) + 2.0),
+ // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+ // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
+ |y| y * (y * y * (y.abs() * 3.0 - 4.0) + 2.0),
|| {
i(
a,
b,
- 0.0,
0.5,
- // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
- // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
- |y| y * (y * y * (y * 3.0 - 4.0) + 2.0),
- || {
- i(
- a,
- b,
- 0.5,
- 1.0,
- // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
- // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4
- |y| -pow4(y - 1.0),
- || 0.0,
- )
- },
+ 1.0,
+ // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
+ // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4
+ |y| -(y - 1.0).powi(4),
+ || 0.0,
)
},
)
@@ -466,41 +426,30 @@
let a = x - β;
let b = x + β;
- // The factor 4 is from normalisation of the original function.
- 4.0 * i(
+ 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),
+ |y| (8.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),
+ 0.5,
+ // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
+ |y| 8.0 * (y.abs() - 1.0) * y * y + (4.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,
- )
- },
+ 1.0,
+ // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
+ |y| -(8.0 / 3.0) * (y - 1.0).powi(3),
+ || 0.0,
)
},
)
@@ -525,6 +474,7 @@
}
*/
+#[replace_float_literals(F::cast_from(literal))]
impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
where
R: Constant<Type = F>,