pointsource_algs: comparison src/kernels/gaussian.rs

-:aec67cdd6b14
+:efa60bc4f743
 let c1 = -(a.min(b + x)); //(-a).max(-x-b);
 let c2 = a.min(b - x);
 if c1 >= c2 {
 0.0
 } else {
-// erf'(z) = (2/√π)*exp(-z^2), and we get extra factor -1/√(2*σ) = -1/t
+// erf'(z) = (2/√π)*exp(-z^2), and we get extra factor -1/(√2*σ) = -1/t
-// from the chain rule
+// from the chain rule (the minus comes from inside c_1 or c_2).
 let de1 = (-(c1/t).powi(2)).exp();
 let de2 = (-(c2/t).powi(2)).exp();
 c_div_t * (de1 - de2)
 }
 })
 self.differential(&y)
 }
 }
 #[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Lipschitz<L1>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+C : Constant<Type=F>,
+S : Constant<Type=F> {
+type FloatType = F;
+fn lipschitz_factor(&self, L1 : L1) -> Option<F> {
+// To get the product Lipschitz factor, we note that 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_i |φ_i|
+//
+// Thus we need 1D Lipschitz factors, and the maximum for φ = θ * ψ.
+//
+// We have
+// θ * ψ(x) = 0 if c_1(x) ≥ c_2(x)
+//          = (1/2)[erf(c_2(x)/(√2σ)) - erf(c_1(x)/(√2σ))] if c_1(x) < c_2(x),
+// where c_1(x) = max{-x-b,-a} = -min{b+x,a} and c_2(x)=min{b-x,a}, C is the Gaussian
+// normalisation factor, and erf(s) = (2/√π) ∫_0^s e^{-t^2} dt.
+// Thus, if c_1(x) < c_2(x) and c_1(y) < c_2(y), we have
+// θ * ψ(x) - θ * ψ(y) = (1/√π)[∫_{c_1(x)/(√2σ)}^{c_1(y)/(√2σ) e^{-t^2} dt
+//                       - ∫_{c_2(x)/(√2σ)}^{c_2(y)/(√2σ)] e^{-t^2} dt]
+// Thus
+// |θ * ψ(x) - θ * ψ(y)| ≤ (1/√π)/(√2σ)(|c_1(x)-c_1(y)|+|c_2(x)-c_2(y)|)
+//                       ≤ 2(1/√π)/(√2σ)|x-y|
+//                       ≤ √2/(√πσ)|x-y|.
+//
+// For the product we also need the value θ * ψ(0), which is
+// (1/2)[erf(min{a,b}/(√2σ))-erf(max{-b,-a}/(√2σ)]
+//  = (1/2)[erf(min{a,b}/(√2σ))-erf(-min{a,b}/(√2σ))]
+//  = erf(min{a,b}/(√2σ))
+//
+// If c_1(x) ≥ c_2(x), then x ∉ [-(a+b), a+b]. If also y is outside that range,
+// θ * ψ(x) = θ * ψ(y). If only y is in the range [-(a+b), a+b], we can replace
+// x by -(a+b) or (a+b), either of which is closer to y and still θ * ψ(x)=0.
+// Thus same calculations as above work for the Lipschitz factor.
+let Convolution(ref ind,
+SupportProductFirst(ref cut,
+ref gaussian)) = self;
+let a = cut.r.value();
+let b = ind.r.value();
+let σ = gaussian.variance.value().sqrt();
+let π = F::PI;
+let t = F::SQRT_2 * σ;
+let l1d = F::SQRT_2 / (π.sqrt() * σ);
+let e0 = F::cast_from(erf((a.min(b) / t).as_()));
+Some(l1d * e0.powi(N as i32-1))
+}
+}
 impl<'a, F : Float, R, C, S, const N : usize> Lipschitz<L2>
 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,
 C : Constant<Type=F>,
 S : Constant<Type=F> {
 type FloatType = F;
+#[inline]
-fn lipschitz_factor(&self, L2 : L2) -> Option<F> {
+fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
-todo!("This requirement some error function work.")
+self.lipschitz_factor(L1).map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
 }
 }
 impl<F : Float, R, C, S, const N : usize>
 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>

Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/gaussian.rs

src/kernels/gaussian.rs