Mercurial > repos > pointsource_algs / file diff

--- 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)