pointsource_algs: comparison src/kernels/base.rs

-:9738b51d90d7
+:4f468d35fa29
 //! Things for constructing new kernels from component kernels and traits for analysing them
+use numeric_literals::replace_float_literals;
 use serde::Serialize;
-use numeric_literals::replace_float_literals;
+use alg_tools::bisection_tree::Support;
+use alg_tools::bounds::{Bounded, Bounds, GlobalAnalysis, LocalAnalysis};
+use alg_tools::instance::{Instance, Space};
+use alg_tools::loc::Loc;
+use alg_tools::mapping::{
+DifferentiableImpl, DifferentiableMapping, LipschitzDifferentiableImpl, Mapping,
+};
+use alg_tools::maputil::{array_init, map1_indexed, map2};
+use alg_tools::norms::*;
+use alg_tools::sets::Cube;
+use alg_tools::sets::SetOrd;
 use alg_tools::types::*;
-use alg_tools::norms::*;
+use anyhow::anyhow;
-use alg_tools::loc::Loc;
-use alg_tools::sets::Cube;
-use alg_tools::bisection_tree::{
-Support,
-Bounds,
-LocalAnalysis,
-GlobalAnalysis,
-Bounded,
-};
-use alg_tools::mapping::{
-Mapping,
-DifferentiableImpl,
-DifferentiableMapping,
-Differential,
-};
-use alg_tools::instance::{Instance, Space};
-use alg_tools::maputil::{array_init, map2, map1_indexed};
-use alg_tools::sets::SetOrd;
 use crate::fourier::Fourier;
 use crate::types::*;
 /// Representation of the product of two kernels.
 ///
 /// The kernels typically implement [`Support`] and [`Mapping`].
 ///
 /// The implementation [`Support`] only uses the [`Support::support_hint`] of the first parameter!
-#[derive(Copy,Clone,Serialize,Debug)]
+#[derive(Copy, Clone, Serialize, Debug)]
 pub struct SupportProductFirst<A, B>(
 /// First kernel
 pub A,
 /// Second kernel
-pub B
+pub B,
 );
-impl<A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
+impl<A, B, F: Float, const N: usize> Mapping<Loc<N, F>> for SupportProductFirst<A, B>
-for SupportProductFirst<A, B>
+where
-where
+A: Mapping<Loc<N, F>, Codomain = F>,
-A : Mapping<Loc<F, N>, Codomain = F>,
+B: Mapping<Loc<N, F>, Codomain = F>,
-B : Mapping<Loc<F, N>, Codomain = F>,
 {
 type Codomain = F;
 #[inline]
-fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
-self.0.apply(x.ref_instance()) * self.1.apply(x)
+x.eval_ref(|r| self.0.apply(r)) * self.1.apply(x)
 }
 }
-impl<A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
+impl<A, B, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for SupportProductFirst<A, B>
-for SupportProductFirst<A, B>
+where
-where
+A: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>, Codomain = F>,
-A : DifferentiableMapping<
+B: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>, Codomain = F>,
-Loc<F, N>,
+{
-DerivativeDomain=Loc<F, N>,
+type Derivative = Loc<N, F>;
-Codomain = F
->,
+#[inline]
-B : DifferentiableMapping<
+fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative {
-Loc<F, N>,
+x.eval_ref(|xr| {
-DerivativeDomain=Loc<F, N>,
+self.0.differential(xr) * self.1.apply(xr) + self.1.differential(xr) * self.0.apply(xr)
-Codomain = F,
+})
->
+}
-{
+}
-type Derivative = Loc<F, N>;
+impl<A, B, M: Copy, F: Float> Lipschitz<M> for SupportProductFirst<A, B>
-#[inline]
+where
-fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+A: Lipschitz<M, FloatType = F> + Bounded<F>,
-let xr = x.ref_instance();
+B: Lipschitz<M, FloatType = F> + Bounded<F>,
-self.0.differential(xr) * self.1.apply(xr) + self.1.differential(xr) * self.0.apply(x)
+{
-}
+type FloatType = F;
-}
+#[inline]
+fn lipschitz_factor(&self, m: M) -> DynResult<F> {
-impl<A, B, M : Copy, F : Float> Lipschitz<M>
-for SupportProductFirst<A, B>
-where A : Lipschitz<M, FloatType = F> + Bounded<F>,
-B : Lipschitz<M, FloatType = F> + Bounded<F> {
-type FloatType = F;
-#[inline]
-fn lipschitz_factor(&self, m : M) -> Option<F> {
 // f(x)g(x) - f(y)g(y) = f(x)[g(x)-g(y)] - [f(y)-f(x)]g(y)
 let &SupportProductFirst(ref f, ref g) = self;
-f.lipschitz_factor(m).map(|l| l * g.bounds().uniform())
+Ok(f.lipschitz_factor(m)? * g.bounds().uniform()
-.zip(g.lipschitz_factor(m).map(|l| l * f.bounds().uniform()))
++ g.lipschitz_factor(m)? * f.bounds().uniform())
-.map(|(a, b)| a + b)
+}
 }
-}
+impl<'a, A, B, M: Copy, Domain, F: Float> LipschitzDifferentiableImpl<Domain, M>
-impl<'a, A, B, M : Copy, Domain, F : Float> Lipschitz<M>
+for SupportProductFirst<A, B>
-for Differential<'a, Domain, SupportProductFirst<A, B>>
+where
-where
+Domain: Space,
-Domain : Space,
+Self: DifferentiableImpl<Domain>,
-A : Clone + DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
+A: DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
-B : Clone + DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
+B: DifferentiableMapping<Domain> + Lipschitz<M, FloatType = F> + Bounded<F>,
-SupportProductFirst<A, B> :  DifferentiableMapping<Domain>,
+SupportProductFirst<A, B>: DifferentiableMapping<Domain>,
-for<'b> A::Differential<'b> : Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType=F>,
+for<'b> A::Differential<'b>: Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType = F>,
-for<'b> B::Differential<'b> : Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType=F>
+for<'b> B::Differential<'b>: Lipschitz<M, FloatType = F> + NormBounded<L2, FloatType = F>,
 {
 type FloatType = F;
 #[inline]
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
 // ∇[gf] = f∇g + g∇f
 // ⟹ ∇[gf](x) - ∇[gf](y) = f(x)∇g(x) + g(x)∇f(x) - f(y)∇g(y) + g(y)∇f(y)
 //                        = f(x)[∇g(x)-∇g(y)] + g(x)∇f(x) - [f(y)-f(x)]∇g(y) + g(y)∇f(y)
 //                        = f(x)[∇g(x)-∇g(y)] + g(x)[∇f(x)-∇f(y)]
 //                          - [f(y)-f(x)]∇g(y) + [g(y)-g(x)]∇f(y)
-let &SupportProductFirst(ref f, ref g) = self.base_fn();
+let &SupportProductFirst(ref f, ref g) = self;
 let (df, dg) = (f.diff_ref(), g.diff_ref());
-[
+Ok([
-df.lipschitz_factor(m).map(|l| l * g.bounds().uniform()),
+df.lipschitz_factor(m)? * g.bounds().uniform(),
-dg.lipschitz_factor(m).map(|l| l * f.bounds().uniform()),
+dg.lipschitz_factor(m)? * f.bounds().uniform(),
-f.lipschitz_factor(m).map(|l| l * dg.norm_bound(L2)),
+f.lipschitz_factor(m)? * dg.norm_bound(L2),
-g.lipschitz_factor(m).map(|l| l * df.norm_bound(L2))
+g.lipschitz_factor(m)? * df.norm_bound(L2),
-].into_iter().sum()
+]
-}
+.into_iter()
-}
+.sum())
+}
+}
-impl<'a, A, B, F : Float, const N : usize> Support<F, N>
-for SupportProductFirst<A, B>
+impl<'a, A, B, F: Float, const N: usize> Support<N, F> for SupportProductFirst<A, B>
 where
-A : Support<F, N>,
+A: Support<N, F>,
-B : Support<F, N>
+B: Support<N, F>,
 {
 #[inline]
-fn support_hint(&self) -> Cube<F, N> {
+fn support_hint(&self) -> Cube<N, F> {
 self.0.support_hint()
 }
 #[inline]
-fn in_support(&self, x : &Loc<F, N>) -> bool {
+fn in_support(&self, x: &Loc<N, F>) -> bool {
 self.0.in_support(x)
 }
 #[inline]
-fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
 self.0.bisection_hint(cube)
 }
 }
-impl<'a, A, B, F : Float> GlobalAnalysis<F, Bounds<F>>
+impl<'a, A, B, F: Float> GlobalAnalysis<F, Bounds<F>> for SupportProductFirst<A, B>
-for SupportProductFirst<A, B>
+where
-where A : GlobalAnalysis<F, Bounds<F>>,
+A: GlobalAnalysis<F, Bounds<F>>,
-B : GlobalAnalysis<F, Bounds<F>> {
+B: GlobalAnalysis<F, Bounds<F>>,
+{
 #[inline]
 fn global_analysis(&self) -> Bounds<F> {
 self.0.global_analysis() * self.1.global_analysis()
 }
 }
-impl<'a, A, B, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
+impl<'a, A, B, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N>
 for SupportProductFirst<A, B>
-where A : LocalAnalysis<F, Bounds<F>, N>,
+where
-B : LocalAnalysis<F, Bounds<F>, N> {
+A: LocalAnalysis<F, Bounds<F>, N>,
-#[inline]
+B: LocalAnalysis<F, Bounds<F>, N>,
-fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+{
+#[inline]
+fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
 self.0.local_analysis(cube) * self.1.local_analysis(cube)
 }
 }
 /// Representation of the sum of two kernels
 ///
 /// The kernels typically implement [`Support`] and [`Mapping`].
 ///
 /// The implementation [`Support`] only uses the [`Support::support_hint`] of the first parameter!
-#[derive(Copy,Clone,Serialize,Debug)]
+#[derive(Copy, Clone, Serialize, Debug)]
 pub struct SupportSum<A, B>(
 /// First kernel
 pub A,
 /// Second kernel
-pub B
+pub B,
 );
-impl<'a, A, B, F : Float, const N : usize> Mapping<Loc<F, N>>
+impl<'a, A, B, F: Float, const N: usize> Mapping<Loc<N, F>> for SupportSum<A, B>
-for SupportSum<A, B>
+where
-where
+A: Mapping<Loc<N, F>, Codomain = F>,
-A : Mapping<Loc<F, N>, Codomain = F>,
+B: Mapping<Loc<N, F>, Codomain = F>,
-B : Mapping<Loc<F, N>, Codomain = F>,
 {
 type Codomain = F;
 #[inline]
-fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Codomain {
+fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Codomain {
-self.0.apply(x.ref_instance()) + self.1.apply(x)
+x.eval_ref(|r| self.0.apply(r)) + self.1.apply(x)
 }
 }
-impl<'a, A, B, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
+impl<'a, A, B, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for SupportSum<A, B>
-for SupportSum<A, B>
+where
-where
+A: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>>,
-A : DifferentiableMapping<
+B: DifferentiableMapping<Loc<N, F>, DerivativeDomain = Loc<N, F>>,
-Loc<F, N>,
+{
-DerivativeDomain = Loc<F, N>
+type Derivative = Loc<N, F>;
->,
-B : DifferentiableMapping<
+#[inline]
-Loc<F, N>,
+fn differential_impl<I: Instance<Loc<N, F>>>(&self, x: I) -> Self::Derivative {
-DerivativeDomain = Loc<F, N>,
+x.eval_ref(|r| self.0.differential(r)) + self.1.differential(x)
->
+}
-{
+}
-type Derivative = Loc<F, N>;
+impl<'a, A, B, F: Float, const N: usize> Support<N, F> for SupportSum<A, B>
+where
-#[inline]
+A: Support<N, F>,
-fn differential_impl<I : Instance<Loc<F, N>>>(&self, x : I) -> Self::Derivative {
+B: Support<N, F>,
-self.0.differential(x.ref_instance()) + self.1.differential(x)
+Cube<N, F>: SetOrd,
-}
+{
-}
+#[inline]
+fn support_hint(&self) -> Cube<N, F> {
-impl<'a, A, B, F : Float, const N : usize> Support<F, N>
-for SupportSum<A, B>
-where A : Support<F, N>,
-B : Support<F, N>,
-Cube<F, N> : SetOrd {
-#[inline]
-fn support_hint(&self) -> Cube<F, N> {
 self.0.support_hint().common(&self.1.support_hint())
 }
 #[inline]
-fn in_support(&self, x : &Loc<F, N>) -> bool {
+fn in_support(&self, x: &Loc<N, F>) -> bool {
 self.0.in_support(x) || self.1.in_support(x)
 }
 #[inline]
-fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
-map2(self.0.bisection_hint(cube),
+map2(
-self.1.bisection_hint(cube),
+self.0.bisection_hint(cube),
-|a, b| a.or(b))
+self.1.bisection_hint(cube),
-}
+|a, b| a.or(b),
-}
+)
+}
-impl<'a, A, B, F : Float> GlobalAnalysis<F, Bounds<F>>
+}
-for SupportSum<A, B>
-where A : GlobalAnalysis<F, Bounds<F>>,
+impl<'a, A, B, F: Float> GlobalAnalysis<F, Bounds<F>> for SupportSum<A, B>
-B : GlobalAnalysis<F, Bounds<F>> {
+where
+A: GlobalAnalysis<F, Bounds<F>>,
+B: GlobalAnalysis<F, Bounds<F>>,
+{
 #[inline]
 fn global_analysis(&self) -> Bounds<F> {
 self.0.global_analysis() + self.1.global_analysis()
 }
 }
-impl<'a, A, B, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
+impl<'a, A, B, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for SupportSum<A, B>
-for SupportSum<A, B>
+where
-where A : LocalAnalysis<F, Bounds<F>, N>,
+A: LocalAnalysis<F, Bounds<F>, N>,
-B : LocalAnalysis<F, Bounds<F>, N>,
+B: LocalAnalysis<F, Bounds<F>, N>,
-Cube<F, N> : SetOrd {
+Cube<N, F>: SetOrd,
-#[inline]
+{
-fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+#[inline]
+fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
 self.0.local_analysis(cube) + self.1.local_analysis(cube)
 }
 }
-impl<F : Float, M : Copy, A, B> Lipschitz<M> for SupportSum<A, B>
+impl<F: Float, M: Copy, A, B> Lipschitz<M> for SupportSum<A, B>
-where A : Lipschitz<M, FloatType = F>,
+where
-B : Lipschitz<M, FloatType = F> {
+A: Lipschitz<M, FloatType = F>,
-type FloatType = F;
+B: Lipschitz<M, FloatType = F>,
+{
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+type FloatType = F;
-match (self.0.lipschitz_factor(m), self.1.lipschitz_factor(m)) {
-(Some(l0), Some(l1)) => Some(l0 + l1),
+fn lipschitz_factor(&self, m: M) -> DynResult<F> {
-_ => None
+Ok(self.0.lipschitz_factor(m)? * self.1.lipschitz_factor(m)?)
 }
 }
-}
+impl<'b, F: Float, M: Copy, A, B, Domain> LipschitzDifferentiableImpl<Domain, M>
-impl<'b, F : Float, M : Copy, A, B, Domain> Lipschitz<M>
+for SupportSum<A, B>
-for Differential<'b, Domain, SupportSum<A, B>>
+where
-where
+Domain: Space,
-Domain : Space,
+Self: DifferentiableImpl<Domain>,
-A : Clone + DifferentiableMapping<Domain, Codomain=F>,
+A: DifferentiableMapping<Domain, Codomain = F>,
-B : Clone + DifferentiableMapping<Domain, Codomain=F>,
+B: DifferentiableMapping<Domain, Codomain = F>,
-SupportSum<A, B> : DifferentiableMapping<Domain, Codomain=F>,
+SupportSum<A, B>: DifferentiableMapping<Domain, Codomain = F>,
-for<'a> A :: Differential<'a> : Lipschitz<M, FloatType = F>,
+for<'a> A::Differential<'a>: Lipschitz<M, FloatType = F>,
-for<'a> B :: Differential<'a> : Lipschitz<M, FloatType = F>
+for<'a> B::Differential<'a>: Lipschitz<M, FloatType = F>,
 {
 type FloatType = F;
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
-let base = self.base_fn();
+Ok(self.0.diff_ref().lipschitz_factor(m)? + self.1.diff_ref().lipschitz_factor(m)?)
-base.0.diff_ref().lipschitz_factor(m)
-.zip(base.1.diff_ref().lipschitz_factor(m))
-.map(|(a, b)| a + b)
 }
 }
 /// Representation of the convolution of two kernels.
 ///
 /// The kernels typically implement [`Support`]s and [`Mapping`].
 //
 /// Trait implementations have to be on a case-by-case basis.
-#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 pub struct Convolution<A, B>(
 /// First kernel
 pub A,
 /// Second kernel
-pub B
+pub B,
 );
-impl<F : Float, M, A, B> Lipschitz<M> for Convolution<A, B>
+impl<F: Float, M, A, B> Lipschitz<M> for Convolution<A, B>
-where A : Norm<F, L1> ,
+where
-B : Lipschitz<M, FloatType = F> {
+A: Norm<L1, F>,
-type FloatType = F;
+B: Lipschitz<M, FloatType = F>,
+{
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+type FloatType = F;
+fn lipschitz_factor(&self, m: M) -> DynResult<F> {
 // For [f * g](x) = ∫ f(x-y)g(y) dy we have
 // [f * g](x) - [f * g](z) = ∫ [f(x-y)-f(z-y)]g(y) dy.
 // Hence |[f * g](x) - [f * g](z)| ≤ ∫ |f(x-y)-f(z-y)|g(y)| dy.
 //                                 ≤ L|x-z| ∫ |g(y)| dy,
 // where L is the Lipschitz factor of f.
-self.1.lipschitz_factor(m).map(|l| l * self.0.norm(L1))
+Ok(self.1.lipschitz_factor(m)? * self.0.norm(L1))
 }
 }
-impl<'b, F : Float, M, A, B, Domain> Lipschitz<M>
+impl<'b, F: Float, M, A, B, Domain> LipschitzDifferentiableImpl<Domain, M> for Convolution<A, B>
-for Differential<'b, Domain, Convolution<A, B>>
+where
-where
+Domain: Space,
-Domain : Space,
+Self: DifferentiableImpl<Domain>,
-A : Clone + Norm<F, L1> ,
+A: Norm<L1, F>,
-Convolution<A, B> : DifferentiableMapping<Domain, Codomain=F>,
+Convolution<A, B>: DifferentiableMapping<Domain, Codomain = F>,
-B : Clone + DifferentiableMapping<Domain, Codomain=F>,
+B: DifferentiableMapping<Domain, Codomain = F>,
-for<'a> B :: Differential<'a> : Lipschitz<M, FloatType = F>
+for<'a> B::Differential<'a>: Lipschitz<M, FloatType = F>,
 {
 type FloatType = F;
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
 // For [f * g](x) = ∫ f(x-y)g(y) dy we have
 // ∇[f * g](x) - ∇[f * g](z) = ∫ [∇f(x-y)-∇f(z-y)]g(y) dy.
 // Hence |∇[f * g](x) - ∇[f * g](z)| ≤ ∫ |∇f(x-y)-∇f(z-y)|g(y)| dy.
 //                                 ≤ L|x-z| ∫ |g(y)| dy,
 // where L is the Lipschitz factor of ∇f.
-let base = self.base_fn();
+self.1
-base.1.diff_ref().lipschitz_factor(m).map(|l| l * base.0.norm(L1))
+.diff_ref()
+.lipschitz_factor(m)
+.map(|l| l * self.0.norm(L1))
 }
 }
 /// Representation of the autoconvolution of a kernel.
 ///
 /// The kernel typically implements [`Support`] and [`Mapping`].
 ///
 /// Trait implementations have to be on a case-by-case basis.
-#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 pub struct AutoConvolution<A>(
 /// The kernel to be autoconvolved
-pub A
+pub A,
 );
-impl<F : Float, M, C> Lipschitz<M> for AutoConvolution<C>
+impl<F: Float, M, C> Lipschitz<M> for AutoConvolution<C>
-where C : Lipschitz<M, FloatType = F> + Norm<F, L1> {
+where
-type FloatType = F;
+C: Lipschitz<M, FloatType = F> + Norm<L1, F>,
+{
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+type FloatType = F;
+fn lipschitz_factor(&self, m: M) -> DynResult<F> {
 self.0.lipschitz_factor(m).map(|l| l * self.0.norm(L1))
 }
 }
-impl<'b, F : Float, M, C, Domain> Lipschitz<M>
+impl<'b, F: Float, M, C, Domain> LipschitzDifferentiableImpl<Domain, M> for AutoConvolution<C>
-for Differential<'b, Domain, AutoConvolution<C>>
+where
-where
+Domain: Space,
-Domain : Space,
+Self: DifferentiableImpl<Domain>,
-C : Clone + Norm<F, L1> + DifferentiableMapping<Domain, Codomain=F>,
+C: Norm<L1, F> + DifferentiableMapping<Domain, Codomain = F>,
-AutoConvolution<C> : DifferentiableMapping<Domain, Codomain=F>,
+AutoConvolution<C>: DifferentiableMapping<Domain, Codomain = F>,
-for<'a> C :: Differential<'a> : Lipschitz<M, FloatType = F>
+for<'a> C::Differential<'a>: Lipschitz<M, FloatType = F>,
 {
 type FloatType = F;
-fn lipschitz_factor(&self, m : M) -> Option<F> {
+fn diff_lipschitz_factor(&self, m: M) -> DynResult<F> {
-let base = self.base_fn();
+self.0
-base.0.diff_ref().lipschitz_factor(m).map(|l| l * base.0.norm(L1))
+.diff_ref()
-}
+.lipschitz_factor(m)
-}
+.map(|l| l * self.0.norm(L1))
+}
+}
 /// Representation a multi-dimensional product of a one-dimensional kernel.
 ///
 /// For $G: ℝ → ℝ$, this is the function $F(x\_1, …, x\_n) := \prod_{i=1}^n G(x\_i)$.
 /// The kernel $G$ typically implements [`Support`] and [`Mapping`]
-/// on [`Loc<F, 1>`]. Then the product implements them on [`Loc<F, N>`].
+/// on [`Loc<1, F>`]. Then the product implements them on [`Loc<N, F>`].
-#[derive(Copy,Clone,Serialize,Debug,Eq,PartialEq)]
+#[derive(Copy, Clone, Serialize, Debug, Eq, PartialEq)]
 #[allow(dead_code)]
-struct UniformProduct<G, const N : usize>(
+struct UniformProduct<G, const N: usize>(
 /// The one-dimensional kernel
-G
+G,
 );
-impl<'a, G, F : Float, const N : usize> Mapping<Loc<F, N>>
+impl<'a, G, F: Float, const N: usize> Mapping<Loc<N, F>> for UniformProduct<G, N>
-for UniformProduct<G, N>
+where
-where
+G: Mapping<Loc<1, F>, Codomain = F>,
-G : Mapping<Loc<F, 1>, Codomain = F>
 {
 type Codomain = F;
 #[inline]
-fn apply<I : Instance<Loc<F, N>>>(&self, x : I) -> F {
+fn apply<I: Instance<Loc<N, F>>>(&self, x: I) -> F {
-x.cow().iter().map(|&y| self.0.apply(Loc([y]))).product()
+x.decompose()
-}
+.iter()
-}
+.map(|&y| self.0.apply(Loc([y])))
+.product()
+}
+}
-impl<'a, G, F : Float, const N : usize> DifferentiableImpl<Loc<F, N>>
-for UniformProduct<G, N>
+impl<'a, G, F: Float, const N: usize> DifferentiableImpl<Loc<N, F>> for UniformProduct<G, N>
 where
-G : DifferentiableMapping<
+G: DifferentiableMapping<Loc<1, F>, DerivativeDomain = F, Codomain = F>,
-Loc<F, 1>,
+{
-DerivativeDomain = F,
+type Derivative = Loc<N, F>;
-Codomain = F,
->
+#[inline]
-{
+fn differential_impl<I: Instance<Loc<N, F>>>(&self, x0: I) -> Loc<N, F> {
-type Derivative = Loc<F, N>;
-#[inline]
-fn differential_impl<I : Instance<Loc<F, N>>>(&self, x0 : I) -> Loc<F, N> {
 x0.eval(|x| {
 let vs = x.map(|y| self.0.apply(Loc([y])));
 product_differential(x, &vs, |y| self.0.differential(Loc([y])))
 })
 }
 /// Helper function to calulate the differential of $f(x)=∏_{i=1}^N g(x_i)$.
 ///
 /// The vector `x` is the location, `vs` consists of the values `g(x_i)`, and
 /// `gd` calculates the derivative `g'`.
 #[inline]
-pub(crate) fn product_differential<F : Float, G : Fn(F) -> F, const N : usize>(
+pub(crate) fn product_differential<F: Float, G: Fn(F) -> F, const N: usize>(
-x : &Loc<F, N>,
+x: &Loc<N, F>,
-vs : &Loc<F, N>,
+vs: &Loc<N, F>,
-gd : G
+gd: G,
-) -> Loc<F, N> {
+) -> Loc<N, F> {
 map1_indexed(x, |i, &y| {
-gd(y) * vs.iter()
+gd(y)
-.zip(0..)
+* vs.iter()
-.filter_map(|(v, j)| (j != i).then_some(*v))
+.zip(0..)
-.product()
+.filter_map(|(v, j)| (j != i).then_some(*v))
-}).into()
+.product()
+})
+.into()
 }
 /// Helper function to calulate the Lipschitz factor of $∇f$ for $f(x)=∏_{i=1}^N g(x_i)$.
 ///
 /// The parameter `bound` is a bound on $|g|_∞$, `lip` is a Lipschitz factor for $g$,
 /// `dbound` is a bound on $|∇g|_∞$, and `dlip` a Lipschitz factor for $∇g$.
 #[inline]
-pub(crate) fn product_differential_lipschitz_factor<F : Float, const N : usize>(
+pub(crate) fn product_differential_lipschitz_factor<F: Float, const N: usize>(
-bound : F,
+bound: F,
-lip : F,
+lip: F,
-dbound : F,
+dbound: F,
-dlip : F
+dlip: F,
-) -> F {
+) -> DynResult<F> {
 // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
 // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{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 M_g^{n-1}|x-y|, where L_g is the Lipschitz factor of g, and
 // M_g a bound on it.
 // With $ψ_i=g for all i, j, it follows that
 // |∇ψ(x) - ∇ψ(y)| ≤ ∑_i L_{∇g} M_g^{n-1} + ∑_{k ≠ i} L_g M_g^{n-2} M_{∇g}
 //                 = n [L_{∇g} M_g^{n-1} + (n-1) L_g M_g^{n-2} M_{∇g}].
 //                 = n M_g^{n-2}[L_{∇g} M_g + (n-1) L_g M_{∇g}].
 if N >= 2 {
-F::cast_from(N) * bound.powi((N-2) as i32)
+Ok(F::cast_from(N)
-* (dlip * bound  + F::cast_from(N-1) * lip * dbound)
+* bound.powi((N - 2) as i32)
-} else if N==1 {
+* (dlip * bound + F::cast_from(N - 1) * lip * dbound))
-dlip
+} else if N == 1 {
+Ok(dlip)
 } else {
-panic!("Invalid dimension")
+Err(anyhow!("Invalid dimension"))
 }
 }
-impl<G, F : Float, const N : usize> Support<F, N>
+impl<G, F: Float, const N: usize> Support<N, F> for UniformProduct<G, N>
-for UniformProduct<G, N>
+where
-where G : Support<F, 1> {
+G: Support<1, F>,
-#[inline]
+{
-fn support_hint(&self) -> Cube<F, N> {
+#[inline]
-let [a] : [[F; 2]; 1] = self.0.support_hint().into();
+fn support_hint(&self) -> Cube<N, F> {
+let [a]: [[F; 2]; 1] = self.0.support_hint().into();
 array_init(|| a.clone()).into()
 }
 #[inline]
-fn in_support(&self, x : &Loc<F, N>) -> bool {
+fn in_support(&self, x: &Loc<N, F>) -> bool {
 x.iter().all(|&y| self.0.in_support(&Loc([y])))
 }
 #[inline]
-fn bisection_hint(&self, cube : &Cube<F, N>) -> [Option<F>; N] {
+fn bisection_hint(&self, cube: &Cube<N, F>) -> [Option<F>; N] {
 cube.map(|a, b| {
 let [h] = self.0.bisection_hint(&[[a, b]].into());
 h
 })
 }
 }
-impl<G, F : Float, const N : usize> GlobalAnalysis<F, Bounds<F>>
+impl<G, F: Float, const N: usize> GlobalAnalysis<F, Bounds<F>> for UniformProduct<G, N>
-for UniformProduct<G, N>
+where
-where G : GlobalAnalysis<F, Bounds<F>> {
+G: GlobalAnalysis<F, Bounds<F>>,
+{
 #[inline]
 fn global_analysis(&self) -> Bounds<F> {
 let g = self.0.global_analysis();
 (0..N).map(|_| g).product()
 }
 }
-impl<G, F : Float, const N : usize> LocalAnalysis<F, Bounds<F>, N>
+impl<G, F: Float, const N: usize> LocalAnalysis<F, Bounds<F>, N> for UniformProduct<G, N>
-for UniformProduct<G, N>
+where
-where G : LocalAnalysis<F, Bounds<F>, 1> {
+G: LocalAnalysis<F, Bounds<F>, 1>,
-#[inline]
+{
-fn local_analysis(&self, cube : &Cube<F, N>) -> Bounds<F> {
+#[inline]
-cube.iter_coords().map(
+fn local_analysis(&self, cube: &Cube<N, F>) -> Bounds<F> {
-|&[a, b]| self.0.local_analysis(&([[a, b]].into()))
+cube.iter_coords()
-).product()
+.map(|&[a, b]| self.0.local_analysis(&([[a, b]].into())))
+.product()
 }
 }
 macro_rules! product_lpnorm {
 ($lp:ident) => {
-impl<G, F : Float, const N : usize> Norm<F, $lp>
+impl<G, F: Float, const N: usize> Norm<$lp, F> for UniformProduct<G, N>
-for UniformProduct<G, N>
+where
-where G : Norm<F, $lp> {
+G: Norm<$lp, F>,
+{
 #[inline]
-fn norm(&self, lp : $lp) -> F {
+fn norm(&self, lp: $lp) -> F {
 self.0.norm(lp).powi(N as i32)
 }
 }
-}
+};
 }
 product_lpnorm!(L1);
 product_lpnorm!(L2);
 product_lpnorm!(Linfinity);
 /// Trait for bounding one kernel with respect to another.
 ///
 /// The type `F` is the scalar field, and `T` another kernel to which `Self` is compared.
-pub trait BoundedBy<F : Num, T> {
+pub trait BoundedBy<F: Num, T> {
 /// Calclate a bounding factor $c$ such that the Fourier transforms $ℱ\[v\] ≤ c ℱ\[u\]$ for
 /// $v$ `self` and $u$ `other`.
 ///
 /// If no such factors exits, `None` is returned.
-fn bounding_factor(&self, other : &T) -> Option<F>;
+fn bounding_factor(&self, other: &T) -> DynResult<F>;
 }
 /// This [`BoundedBy`] implementation bounds $(uv) * (uv)$ by $(ψ * ψ) u$.
 #[replace_float_literals(F::cast_from(literal))]
-impl<F, C, BaseP>
+impl<F, C, BaseP> BoundedBy<F, SupportProductFirst<AutoConvolution<C>, BaseP>>
-BoundedBy<F, SupportProductFirst<AutoConvolution<C>, BaseP>>
+for AutoConvolution<SupportProductFirst<C, BaseP>>
-for AutoConvolution<SupportProductFirst<C, BaseP>>
+where
-where F : Float,
+F: Float,
-C : Clone + PartialEq,
+C: PartialEq,
-BaseP :  Fourier<F> + PartialOrd, // TODO: replace by BoundedBy,
+BaseP: Fourier<F> + PartialOrd, // TODO: replace by BoundedBy,
-<BaseP as Fourier<F>>::Transformed : Bounded<F> + Norm<F, L1> {
+<BaseP as Fourier<F>>::Transformed: Bounded<F> + Norm<L1, F>,
+{
-fn bounding_factor(&self, kernel : &SupportProductFirst<AutoConvolution<C>, BaseP>) -> Option<F> {
+fn bounding_factor(
+&self,
+kernel: &SupportProductFirst<AutoConvolution<C>, BaseP>,
+) -> DynResult<F> {
 let SupportProductFirst(AutoConvolution(ref cutoff2), base_spread2) = kernel;
 let AutoConvolution(SupportProductFirst(ref cutoff, ref base_spread)) = self;
 let v̂ = base_spread.fourier();
 // Verify that the cut-off and ideal physical model (base spread) are the same.
-if cutoff == cutoff2
+if cutoff == cutoff2 && base_spread <= base_spread2 && v̂.bounds().lower() >= 0.0 {
-&& base_spread <= base_spread2
-&& v̂.bounds().lower() >= 0.0 {
 // Calculate the factor between the convolution approximation
 // `AutoConvolution<SupportProductFirst<C, BaseP>>` of $A_*A$ and the
 // kernel of the seminorm. This depends on the physical model P being
 // `SupportProductFirst<C, BaseP>` with the kernel `K` being
 // a `SupportSum` involving `SupportProductFirst<AutoConvolution<C>, BaseP>`.
-Some(v̂.norm(L1))
+Ok(v̂.norm(L1))
 } else {
 // We cannot compare
-None
+Err(anyhow!("Incomprable kernels"))
 }
 }
 }
-impl<F : Float, A, B, C> BoundedBy<F, SupportSum<B, C>> for A
+impl<F: Float, A, B, C> BoundedBy<F, SupportSum<B, C>> for A
-where A : BoundedBy<F, B>,
+where
-C : Bounded<F> {
+A: BoundedBy<F, B>,
+C: Bounded<F>,
+{
 #[replace_float_literals(F::cast_from(literal))]
-fn bounding_factor(&self, SupportSum(ref kernel1, kernel2) : &SupportSum<B, C>) -> Option<F> {
+fn bounding_factor(&self, SupportSum(ref kernel1, kernel2): &SupportSum<B, C>) -> DynResult<F> {
 if kernel2.bounds().lower() >= 0.0 {
 self.bounding_factor(kernel1)
 } else {
-None
+Err(anyhow!("Component kernel not lower-bounded by zero"))
 }
 }
 }
 /// Generates on $[a, b]$ [`Support::support_hint`] for a symmetric interval $[-r, r]$.
 ///
 /// It will attempt to place the subdivision point at $-r$ or $r$.
 /// If neither of these points lies within $[a, b]$, `None` is returned.
 #[inline]
-pub(super) fn symmetric_interval_hint<F : Float>(r : F, a : F, b : F) -> Option<F> {
+pub(super) fn symmetric_interval_hint<F: Float>(r: F, a: F, b: F) -> Option<F> {
 if a < -r && -r < b {
 Some(-r)
 } else if a < r && r < b {
 Some(r)
 } else {
 /// gives the longer length for the shorter of the two subintervals. If none of these points
 /// lies within $[a, b]$, or the resulting interval would be shorter than $0.3r$, `None` is
 /// returned.
 #[replace_float_literals(F::cast_from(literal))]
 #[inline]
-pub(super) fn symmetric_peak_hint<F : Float>(r : F, a : F, b : F) -> Option<F> {
+pub(super) fn symmetric_peak_hint<F: Float>(r: F, a: F, b: F) -> Option<F> {
 let stage1 = if a < -r {
 if b <= -r {
 None
 } else if a + r < -b {
 Some(-r)
 };
 // Ignore stage1 hint if either side of subdivision would be just a small fraction of the
 // interval
 match stage1 {
-Some(h) if (h - a).min(b-h) >= 0.3 * r => Some(h),
+Some(h) if (h - a).min(b - h) >= 0.3 * r => Some(h),
-_ => None
+_ => None,
 }
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/base.rs

src/kernels/base.rs