Mercurial > repos > pointsource_algs / changeset

--- a/src/kernels/base.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/kernels/base.rs	Fri Oct 07 13:11:08 2022 +0300
@@ -13,6 +13,7 @@
     LocalAnalysis,
     GlobalAnalysis,
     Bounded,
+    Taylor2ModelParams
 };
 use alg_tools::mapping::Apply;
 use alg_tools::maputil::{array_init, map2};
@@ -402,3 +403,62 @@
         _ => None
     }
 }
+
+pub trait Product2Taylor<F : Float, const N : usize> : Sized {
+    fn product2taylor(self) -> Taylor2ModelParams<F, N>;
+
+    fn product2taylor_scaled(self, sc : F) -> Taylor2ModelParams<F, N> {
+        let (a, b, c, d) = self.product2taylor();
+        (a / sc, b / sc, c.map(|x| x / sc), d / sc)
+    }
+}
+
+impl<F: Float> Product2Taylor<F, 1> for Loc<Loc<F, 3>, 1> {
+    fn product2taylor(self) -> Taylor2ModelParams<F, 1> {
+        let Loc([Loc([a, b, c])]) = self.into();
+        (a, Loc([b]), Loc([Loc([c])]), F::ZERO)
+    }
+}
+
+impl<F: Float> Product2Taylor<F, 2> for Loc<Loc<F, 3>, 2> {
+    fn product2taylor(self) -> Taylor2ModelParams<F, 2> {
+        let Loc([Loc([a1, b1, c1]), Loc([a2, b2, c2])]) = self;
+        (a1 * a2,
+         Loc([b1 * a2, a1 * b2]),
+         Loc([Loc([c1 * a2, b1 * b2]), Loc([b1 * b2, a1 * c2])]),
+         F::ZERO)
+    }
+}
+
+        
+        // // Zeroeth order term
+        // let a = factors.iter().map(|&(p0, _p1, _p2)| p0).product() / sc;
+        // // First order term
+        // let b = factors.iter().zip(0..).map(|(&(_p0, p1, _p2), i)| {
+        //     (p1 / sc) * factors.iter()
+        //                        .zip(0..)
+        //                        .filter(|&(_, j)| j != i)
+        //                        .map(|(&(q0, _q1, _q2), _)| q0)
+        //                        .product()
+        // }).into();
+
+        // // Second order term
+        // let c = factors.iter().zip(0..).map(|(&(_p0, p1, p2), i)| {
+        //     factors.iter().zip(0..).map(|(&(_q0, q1, _q2), j)| {
+        //         if i == j {
+        //             (p2 /sc) * factors.iter()
+        //                               .zip(0..)
+        //                               .filter(|&(_, k)| k != i)
+        //                               .map(|(&(w0, _w1, _w2), _)| w0)
+        //                               .product()
+        //         } else {
+        //             (p1 * q1 / sc) * factors.iter()
+        //                                     .zip(0..)
+        //                                     .filter(|&(_, k)| k != i && k != j)
+        //                                     .map(|(&(w0, _w1, _w2), _)| w0)
+        //                                     .product()
+        //         }
+        //     }).into()
+        // }).into();
+        //
+        // (a, b, c, 0.0)

--- a/src/kernels/gaussian.rs	Thu Dec 01 23:07:35 2022 +0200
+++ b/src/kernels/gaussian.rs	Fri Oct 07 13:11:08 2022 +0300
@@ -16,6 +16,8 @@
     GlobalAnalysis,
     Weighted,
     Bounded,
+    Taylor2Model,
+    Taylor2ModelParams,
 };
 use alg_tools::mapping::Apply;
 use alg_tools::maputil::array_init;
@@ -224,6 +226,49 @@
     }
 }
 
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F : Float, R, C, S, const N : usize> Taylor2Model<F, N>
+for Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
+where R : Constant<Type=F>,
+      C : Constant<Type=F>,
+      S : Constant<Type=F>,
+      Loc<Loc<F, 3>, N> : Product2Taylor<F, N> {
+
+    fn taylor2model(&self, y : Loc<F, N>) -> Taylor2ModelParams<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 sc = gaussian.scale();
+        
+        // This is just a product of one-dimensional versions
+        y.map(|x| {
+            let c1 = -(a.min(b + x)); //(-a).max(-x-b);
+            let c2 = a.min(b - x);
+            if c1 >= c2 {
+                Loc([0.0, 0.0, 0.0])
+            } else {
+                let s1 = c1 / t;
+                let s2 = c2 / t;
+                let e1 = F::cast_from(erf(s1.as_()));
+                let e2 = F::cast_from(erf(s2.as_()));
+                let d1 = F::FRAC_2_SQRT_PI * (-s1*s1).exp();
+                let d2 = F::FRAC_2_SQRT_PI * (-s2*s2).exp();
+                let dd1 = -2.0 * s1 * d1;
+                let dd2 = -2.0 * s2 * s2;
+                debug_assert!(e2 >= e1);
+                Loc([c * (e2 - e1), c * (d2 - d1), c * (dd2 - dd1)])
+            }
+        }).product2taylor_scaled(sc)
+    }
+
+}
+
 impl<F : Float, R, C, S, const N : usize>
 Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>>
 where R : Constant<Type=F>,