Mercurial > repos > pointsource_algs / file diff

--- a/src/forward_model/sensor_grid.rs	Thu Jan 23 23:39:40 2025 +0100
+++ b/src/forward_model/sensor_grid.rs	Fri Jan 31 18:00:17 2025 -0500
@@ -38,12 +38,10 @@
     AutoConvolution,
     BoundedBy,
 };
-use crate::types::L2Squared;
-use crate::transport::TransportLipschitz;
 use crate::preadjoint_helper::PreadjointHelper;
 use super::{
     ForwardModel,
-    LipschitzValues,
+    BoundedCurvature,
     AdjointProductBoundedBy
 };
 use crate::frank_wolfe::FindimQuadraticModel;
@@ -309,6 +307,7 @@
     }
 }
 
+/*
 #[replace_float_literals(F::cast_from(literal))]
 impl<'a, F, S, P, BT, const N : usize> LipschitzValues
 for SensorGridPreadjoint<'a, SensorGrid<F, S, P, BT, N>, F, N>
@@ -340,6 +339,50 @@
         fw.base_sensor.diff_ref().lipschitz_factor(L2).map(|l| (2.0 * n).sqrt() * l)
     }
 }
+*/
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<'a, F, S, P, BT, const N : usize> BoundedCurvature
+for SensorGrid<F, S, P, BT, N>
+where F : Float,
+      BT : SensorGridBT<F, S, P, N>,
+      S : Sensor<F, N>,
+      P : Spread<F, N>,
+      Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType=F> + DifferentiableMapping<Loc<F,N>> + LocalAnalysis<F, BT::Agg, N>,
+      for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<F,N>>>::Differential<'b> : Lipschitz<L2, FloatType=F>,
+{
+
+    type FloatType = F;
+
+    /// Returns a bound $ℓ_F$ on the curvature
+    /// $$
+    ///     𝒦_F(μ, γ) = ∫ B_{F'(μ)} dγ +  B_F(μ, μ+Δ).
+    /// $$
+    /// such that $𝒦_F(μ, γ)  ≤ ℓ_F ∫  c_2 d|γ|$.
+    ///
+    /// For $F(μ)=(1/2)‖Aμ-b‖^2$, we have $B_F(μ, μ+Δ)=(1/2)‖AΔ‖^2$, where $Δ = (π_♯^1-π_♯^0)γ$.
+    /// So we use Lemma 3.8 for that, bounding
+    /// $(1/2)‖AΔ‖^2 ≤ Θ ∫ c_2 dγ$ for $Θ=2N_ψML_ψ^2$, where
+    ///   * $L_ψ$ is the 2-norm Lipschitz factor of $ψ$ (sensor * base_spread), and
+    ///   * $N_ψ$ the maximum overlap,
+    ///   * M is a bound on $|γ|(Ω^2)$.
+    ///
+    /// We also have $B_{F'(μ)}(x, y) = v(y) - v(x) ⟨∇v(x), x-y⟩$ for $v(x)=A^*(Aμ-b)$.
+    /// This we want the Lipschitz factor of $∇v$.
+    /// By Example 4.15, it makes sense to estimate this by $√(2N_ψ)L_{∇ψ}‖b‖$, where
+    /// $L_{∇ψ}$ is the Lipshitz factor of $∇ψ$.
+    fn curvature_bound_components(&self) -> (Option<Self::FloatType>, Option<Self::FloatType>) {
+        let n_ψ = self.max_overlapping();
+        let ψ_diff_lip = self.base_sensor.diff_ref().lipschitz_factor(L2);
+        let ψ_lip = self.base_sensor.lipschitz_factor(L2);
+        let a = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l);
+        let b = ψ_lip.map(|l| 2.0 * n_ψ * l.powi(2));
+
+        (a, b)
+    }
+}
+
+
 
 #[derive(Clone,Debug)]
 pub struct SensorGridSupportGenerator<F, S, P, const N : usize>
@@ -491,28 +534,6 @@
     }
 }
 
-#[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 + ToNalgebraRealField,
-      BT : SensorGridBT<F, S, P, N>,
-      S : Sensor<F, N>,
-      P : Spread<F, N>,
-      Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F>
-{
-    type FloatType = F;
-
-    fn transport_lipschitz_factor(&self, L2Squared : L2Squared) -> Self::FloatType {
-        // 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.
-        // The factors two comes from Lipschitz estimates having two possible
-        // points of overlap.
-        let l = self.base_sensor.lipschitz_factor(L2).unwrap();
-        2.0 * self.max_overlapping() * l.powi(2)
-    }
-}
-
-
 macro_rules! make_sensorgridsupportgenerator_scalarop_rhs {
     ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
         impl<F, S, P, const N : usize>
@@ -630,3 +651,4 @@
       P : Spread<F, N>,
       Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>,
 { }
+