Mercurial > repos > pointsource_algs / file diff

--- a/src/forward_model.rs	Tue Aug 01 10:32:12 2023 +0300
+++ b/src/forward_model.rs	Thu Aug 29 00:00:00 2024 -0500
@@ -23,6 +23,7 @@
 use alg_tools::nalgebra_support::ToNalgebraRealField;
 use alg_tools::tabledump::write_csv;
 use alg_tools::error::DynError;
+use alg_tools::maputil::map2;
 
 use crate::types::*;
 use crate::measures::*;
@@ -594,15 +595,25 @@
 #[replace_float_literals(F::cast_from(literal))]
 impl<F, BT, S, P, const N : usize> TransportLipschitz<L2Squared>
 for SensorGrid<F, S, P, BT, N>
-where F : Float + nalgebra::RealField + ToNalgebraRealField,
+where F : Float + ToNalgebraRealField,
       BT : SensorGridBT<F, S, P, N>,
       S : Sensor<F, N>,
       P : Spread<F, N>,
-      Convolution<S, P> : Spread<F, N> + Lipschitz<L2> {
+      Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F> {
     type FloatType = F;
 
     fn transport_lipschitz_factor(&self, L2Squared : L2Squared) -> Self::FloatType {
-        todo!("Unimplemented")
+        // We estimate the factor by N_ψL^2, where L is the 2-norm Lipschitz factor of
+        // the base sensor (sensor * base_spread), and N_ψ the maximum overlap.
+        let l = self.base_sensor.lipschitz_factor(L2).unwrap();
+        let w = self.base_sensor.support_hint().width();
+        let d = map2(self.domain.width(), &self.sensor_count, |wi, &i| wi/F::cast_from(i));
+        let n = w.iter()
+                 .zip(d.iter())
+                 .map(|(&wi, &di)| (wi/di).ceil())
+                 .reduce(F::mul)
+                 .unwrap();
+        2.0 * n * l.powi(2)
     }
 }