pointsource_algs: comparison src/kernels/gaussian.rs

-:bd13c2ae3450
+:56c8adc32b09
 LocalAnalysis,
 GlobalAnalysis,
 Weighted,
 Bounded,
 };
-use alg_tools::mapping::Apply;
+use alg_tools::mapping::{Apply, Differentiable};
 use alg_tools::maputil::array_init;
+use crate::types::Lipschitz;
 use crate::fourier::Fourier;
 use super::base::*;
 use super::ball_indicator::CubeIndicator;
 /// Storage presentation of the the anisotropic gaussian kernel of `variance` $σ^2$.
 }
 impl<S, const N : usize> Apply<Loc<S::Type, N>> for Gaussian<S, N>
 where S : Constant {
 type Output = S::Type;
+#[inline]
+fn apply(&self, x : Loc<S::Type, N>) -> Self::Output {
+self.apply(&x)
+}
+}
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<'a, S, const N : usize> Differentiable<&'a Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+type Output = Loc<S::Type, N>;
+#[inline]
+fn differential(&self, x : &'a Loc<S::Type, N>) -> Self::Output {
+x * (self.apply(x) / self.variance.value())
+}
+}
+impl<S, const N : usize> Differentiable<Loc<S::Type, N>> for Gaussian<S, N>
+where S : Constant {
+type Output = Loc<S::Type, N>;
 // This is not normalised to neither to have value 1 at zero or integral 1
 // (unless the cut-off ε=0).
 #[inline]
-fn apply(&self, x : Loc<S::Type, N>) -> Self::Output {
+fn differential(&self, x : Loc<S::Type, N>) -> Self::Output {
-self.apply(&x)
+x * (self.apply(&x) / self.variance.value())
 }
 }
+#[replace_float_literals(S::Type::cast_from(literal))]
+impl<S, const N : usize> Lipschitz<L2> for Gaussian<S, N>
+where S : Constant {
+type FloatType = S::Type;
+fn lipschitz_factor(&self, L2 : L2) -> Option<Self::FloatType> {
+// f(x)=f_1(‖x‖_2/σ) * √(2π) / √(2πσ)^N, where f_1 is one-dimensional Gaussian with
+// variance 1. The Lipschitz factor of f_1 is e^{-1/2}/√(2π), see, e.g.,
+// https://math.stackexchange.com/questions/3630967/is-the-gaussian-density-lipschitz-continuous
+// Thus the Lipschitz factor we want is e^{-1/2} / (√(2πσ)^N * σ).
+Some((-0.5).exp() / (self.scale() * self.variance.value().sqrt()))
+}
+}
 #[replace_float_literals(S::Type::cast_from(literal))]
 impl<'a, S, const N : usize> Gaussian<S, N>
 where S : Constant {
 /// where $a>0$ and $f$ is a gaussian kernel on $ℝ^n$.
 pub type BasicCutGaussian<C, S, const N : usize> = SupportProductFirst<CubeIndicator<C, N>,
 Gaussian<S, N>>;
-/// This implements $χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$
+/// This implements $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$ and $f$ is
-/// where $a,b>0$ and $f$ is a gaussian kernel on $ℝ^n$.
+/// a gaussian kernel on $ℝ^n$. For an expression for $g$, see Lemma 3.9 in the manuscript.
 #[replace_float_literals(F::cast_from(literal))]
 impl<'a, F : Float, R, C, S, const N : usize> Apply<&'a Loc<F, N>>
 for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,
 C : Constant<Type=F>,
 fn apply(&self, y : Loc<F, N>) -> F {
 self.apply(&y)
 }
 }
+/// This implements the differential of $g := χ\_{[-b, b]^n} \* (f χ\_{[-a, a]^n})$ where $a,b>0$
+/// and $f$ is a gaussian kernel on $ℝ^n$. For an expression for the value of $g$, from which the
+/// derivative readily arises (at points of differentiability), see Lemma 3.9 in the manuscript.
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Differentiable<&'a Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+C : Constant<Type=F>,
+S : Constant<Type=F> {
+type Output = Loc<F, N>;
+#[inline]
+fn differential(&self, y : &'a Loc<F, N>) -> Loc<F, N> {
+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 c = σ * (8.0/π).sqrt();
+let cd = (8.0).sqrt(); // σ * (8.0/π).sqrt() / t * (√2/π)
+// Calculate the values for all component functions of the
+// product. This is just the loop from apply above.
+let unscaled_vs = y.map(|x| {
+let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+let c2 = a.min(b - x);
+if c1 >= c2 {
+0.0
+} else {
+let e1 = F::cast_from(erf((c1 / t).as_()));
+let e2 = F::cast_from(erf((c2 / t).as_()));
+debug_assert!(e2 >= e1);
+c * (e2 - e1)
+}
+});
+// This computes the gradient for each coordinate
+product_differential(y, &unscaled_vs, |x| {
+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
+// from the chain rule
+let de1 = (-(c1/t).powi(2)).exp();
+let de2 = (-(c2/t).powi(2)).exp();
+cd * (de1 - de2)
+}
+}) / gaussian.scale()
+}
+}
+impl<F : Float, R, C, S, const N : usize> Differentiable<Loc<F, N>>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+C : Constant<Type=F>,
+S : Constant<Type=F> {
+type Output = Loc<F, N>;
+#[inline]
+fn differential(&self, y : Loc<F, N>) -> Loc<F, N> {
+self.differential(&y)
+}
+}
+#[replace_float_literals(F::cast_from(literal))]
+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;
+fn lipschitz_factor(&self, L2 : L2) -> Option<F> {
+todo!("This requirement some error function work.")
+}
+}
 impl<F : Float, R, C, S, const N : usize>
 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,
 C : Constant<Type=F>,
 S : Constant<Type=F> {

Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/gaussian.rs

src/kernels/gaussian.rs