diff -r 9738b51d90d7 -r 4f468d35fa29 src/forward_model/sensor_grid.rs --- a/src/forward_model/sensor_grid.rs Sun Apr 27 15:03:51 2025 -0500 +++ b/src/forward_model/sensor_grid.rs Thu Feb 26 11:38:43 2026 -0500 @@ -2,13 +2,17 @@ Sensor grid forward model */ -use nalgebra::base::{DMatrix, DVector}; -use numeric_literals::replace_float_literals; -use std::iter::Zip; -use std::ops::RangeFrom; - +use super::{BasicCurvatureBoundEstimates, ForwardModel}; +use crate::frank_wolfe::FindimQuadraticModel; +use crate::kernels::{AutoConvolution, BoundedBy, Convolution}; +use crate::measures::{DiscreteMeasure, Radon, RNDM}; +use crate::preadjoint_helper::PreadjointHelper; +use crate::seminorms::{ConvolutionOp, SimpleConvolutionKernel}; +use crate::types::*; use alg_tools::bisection_tree::*; -use alg_tools::error::DynError; +use alg_tools::bounds::Bounded; +use alg_tools::error::{DynError, DynResult}; +use alg_tools::euclidean::Euclidean; use alg_tools::instance::Instance; use alg_tools::iter::{MapX, Mappable}; use alg_tools::lingrid::*; @@ -18,79 +22,74 @@ use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::{Linfinity, Norm, L1, L2}; use alg_tools::tabledump::write_csv; +use anyhow::anyhow; +use nalgebra::base::{DMatrix, DVector}; +use numeric_literals::replace_float_literals; +use std::iter::Zip; +use std::ops::RangeFrom; -use super::{AdjointProductBoundedBy, BoundedCurvature, ForwardModel}; -use crate::frank_wolfe::FindimQuadraticModel; -use crate::kernels::{AutoConvolution, BoundedBy, Convolution}; -use crate::measures::{DiscreteMeasure, Radon}; -use crate::preadjoint_helper::PreadjointHelper; -use crate::seminorms::{ConvolutionOp, SimpleConvolutionKernel}; -use crate::types::*; - -type RNDM = DiscreteMeasure, F>; - -pub type ShiftedSensor = Shift, F, N>; +pub type ShiftedSensor = Shift, N, F>; /// Trait for physical convolution models. Has blanket implementation for all cases. -pub trait Spread: - 'static + Clone + Support + RealMapping + Bounded +pub trait Spread: + 'static + Clone + Support + RealMapping + Bounded { } -impl Spread for T +impl Spread for T where F: Float, - T: 'static + Clone + Support + Bounded + RealMapping, + T: 'static + Clone + Support + Bounded + RealMapping, { } /// Trait for compactly supported sensors. Has blanket implementation for all cases. -pub trait Sensor: - Spread + Norm + Norm +pub trait Sensor: + Spread + Norm + Norm { } -impl Sensor for T +impl Sensor for T where F: Float, - T: Spread + Norm + Norm, + T: Spread + Norm + Norm, { } pub trait SensorGridBT: - Clone + BTImpl> + Clone + BTImpl> where F: Float, - S: Sensor, - P: Spread, + S: Sensor, + P: Spread, { } impl SensorGridBT for T where - T: Clone + BTImpl>, + T: Clone + BTImpl>, F: Float, - S: Sensor, - P: Spread, + S: Sensor, + P: Spread, { } // We need type alias bounds to access associated types #[allow(type_alias_bounds)] pub type SensorGridBTFN, const N: usize> = - BTFN, BT, N>; + BTFN, BT, N>; /// Sensor grid forward model #[derive(Clone)] pub struct SensorGrid where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, BT: SensorGridBT, { - domain: Cube, + domain: Cube, sensor_count: [usize; N], sensor: S, spread: P, @@ -102,16 +101,16 @@ where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, { /// Create a new sensor grid. /// /// The parameter `depth` indicates the search depth of the created [`BT`]s /// for the adjoint values. pub fn new( - domain: Cube, + domain: Cube, sensor_count: [usize; N], sensor: S, spread: P, @@ -141,12 +140,12 @@ where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { /// Return the grid of sensor locations. - pub fn grid(&self) -> LinGrid { + pub fn grid(&self) -> LinGrid { lingrid_centered(&self.domain, &self.sensor_count) } @@ -157,7 +156,7 @@ /// Constructs a sensor shifted by `x`. #[inline] - fn shifted_sensor(&self, x: Loc) -> ShiftedSensor { + fn shifted_sensor(&self, x: Loc) -> ShiftedSensor { self.base_sensor.clone().shift(x) } @@ -180,44 +179,44 @@ } } -impl Mapping> for SensorGrid +impl Mapping> for SensorGrid where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { type Codomain = DVector; #[inline] - fn apply>>(&self, μ: I) -> DVector { + fn apply>>(&self, μ: I) -> DVector { let mut y = self._zero_observable(); self.apply_add(&mut y, μ); y } } -impl Linear> for SensorGrid +impl Linear> for SensorGrid where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { } #[replace_float_literals(F::cast_from(literal))] -impl GEMV, DVector> for SensorGrid +impl GEMV, DVector> for SensorGrid where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { - fn gemv>>(&self, y: &mut DVector, α: F, μ: I, β: F) { + fn gemv>>(&self, y: &mut DVector, α: F, μ: I, β: F) { let grid = self.grid(); if β == 0.0 { y.fill(0.0) @@ -227,38 +226,42 @@ if α == 1.0 { self.apply_add(y, μ) } else { - for δ in μ.ref_instance() { - for &d in self.bt.iter_at(&δ.x) { - let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); - y[d] += sensor.apply(&δ.x) * (α * δ.α); + μ.eval_ref(|μr| { + for δ in μr { + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + y[d] += sensor.apply(&δ.x) * (α * δ.α); + } } - } + }) } } - fn apply_add>>(&self, y: &mut DVector, μ: I) { + fn apply_add>>(&self, y: &mut DVector, μ: I) { let grid = self.grid(); - for δ in μ.ref_instance() { - for &d in self.bt.iter_at(&δ.x) { - let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); - y[d] += sensor.apply(&δ.x) * δ.α; + μ.eval_ref(|μr| { + for δ in μr { + for &d in self.bt.iter_at(&δ.x) { + let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); + y[d] += sensor.apply(&δ.x) * δ.α; + } } - } + }) } } -impl BoundedLinear, Radon, L2, F> +impl BoundedLinear, Radon, L2, F> for SensorGrid where F: Float, - BT: SensorGridBT>, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, + BT: SensorGridBT, + S: Sensor, + P: Spread, + Convolution: Spread, { /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped /// with the Radon norm, and $ℝ^n$ with the Euclidean norm. - fn opnorm_bound(&self, _: Radon, _: L2) -> F { + fn opnorm_bound(&self, _: Radon, _: L2) -> DynResult { // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩ // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ @@ -271,22 +274,22 @@ // = |φ|_∞ √N_ψ |μ|_ℳ. // Hence let n = self.max_overlapping(); - self.base_sensor.bounds().uniform() * n.sqrt() + Ok(self.base_sensor.bounds().uniform() * n.sqrt()) } } -type SensorGridPreadjoint<'a, A, F, const N: usize> = PreadjointHelper<'a, A, RNDM>; +type SensorGridPreadjoint<'a, A, F, const N: usize> = PreadjointHelper<'a, A, RNDM>; -impl Preadjointable, DVector> +impl Preadjointable, DVector> for SensorGrid where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, { - type PreadjointCodomain = BTFN, BT, N>; + type PreadjointCodomain = BTFN, BT, N>; type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where @@ -303,10 +306,10 @@ for SensorGridPreadjoint<'a, SensorGrid, F, N> where F : Float, BT : SensorGridBT, - S : Sensor, - P : Spread, - Convolution : Spread + Lipschitz + DifferentiableMapping> + LocalAnalysis, - for<'b> as DifferentiableMapping>>::Differential<'b> : Lipschitz, + S : Sensor, + P : Spread, + Convolution : Spread + Lipschitz + DifferentiableMapping> + LocalAnalysis, + for<'b> as DifferentiableMapping>>::Differential<'b> : Lipschitz, { type FloatType = F; @@ -332,55 +335,49 @@ */ #[replace_float_literals(F::cast_from(literal))] -impl<'a, F, S, P, BT, const N: usize> BoundedCurvature for SensorGrid +impl<'a, F, S, P, BT, const N: usize> BasicCurvatureBoundEstimates for SensorGrid where - F: Float, + F: Float + ToNalgebraRealField, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + S: Sensor, + P: Spread, + DVector: Euclidean, + Convolution: Spread + Lipschitz - + DifferentiableMapping> + + DifferentiableMapping> + LocalAnalysis, - for<'b> as DifferentiableMapping>>::Differential<'b>: + for<'b> as DifferentiableMapping>>::Differential<'b>: Lipschitz, { - type FloatType = F; - - /// Returns factors $ℓ_F$ and $Θ²$ such that - /// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ)+F'(μ+Δ)|Δ⟩ ≤ Θ²|γ|(c_2)‖γ‖$, - /// where $Δ=(π_♯^1-π_♯^0)γ$. - /// - /// See Lemma 3.8, Lemma 5.10, Remark 5.14, and Example 5.15. - fn curvature_bound_components(&self) -> (Option, Option) { + fn basic_curvature_bound_components(&self) -> (DynResult, DynResult) { 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 ℓ_F = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l); - let θ2 = ψ_lip.map(|l| 4.0 * n_ψ * l.powi(2)); + let ℓ_F0 = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l); + let Θ2 = ψ_lip.map(|l| 4.0 * n_ψ * l.powi(2)); - (ℓ_F, θ2) + (ℓ_F0, Θ2) } } #[derive(Clone, Debug)] -pub struct SensorGridSupportGenerator +pub struct SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, + S: Sensor, + P: Spread, { base_sensor: Convolution, - grid: LinGrid, + grid: LinGrid, weights: DVector, } -impl SensorGridSupportGenerator +impl SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { #[inline] fn construct_sensor(&self, id: usize, w: F) -> Weighted, F> { @@ -397,12 +394,12 @@ } } -impl SupportGenerator for SensorGridSupportGenerator +impl SupportGenerator for SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { type Id = usize; type SupportType = Weighted, F>; @@ -434,14 +431,14 @@ } } -impl ForwardModel, F>, F> +impl ForwardModel, F>, F> for SensorGrid where F: Float + ToNalgebraRealField + nalgebra::RealField, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, { type Observable = DVector; @@ -456,17 +453,17 @@ } } -impl FindimQuadraticModel, F> for SensorGrid +impl FindimQuadraticModel, F> for SensorGrid where F: Float + ToNalgebraRealField + nalgebra::RealField, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, { fn findim_quadratic_model( &self, - μ: &DiscreteMeasure, F>, + μ: &DiscreteMeasure, F>, b: &Self::Observable, ) -> (DMatrix, DVector) { assert_eq!(b.len(), self.n_sensors()); @@ -483,29 +480,29 @@ } } -/// Implements the calculation a factor $L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model +/// Implements the calculation a factor $√L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model /// and $𝒟$ a seminorm of suitable form. /// /// **This assumes (but does not check) that the sensors are not overlapping.** #[replace_float_literals(F::cast_from(literal))] -impl AdjointProductBoundedBy, ConvolutionOp> - for SensorGrid +impl<'a, F, BT, S, P, K, const N: usize> + BoundedLinear, &'a ConvolutionOp, L2, F> for SensorGrid where - F: Float + nalgebra::RealField + ToNalgebraRealField, + F: Float + ToNalgebraRealField, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread, - K: SimpleConvolutionKernel, + S: Sensor, + P: Spread, + Convolution: Spread, + K: SimpleConvolutionKernel, AutoConvolution

: BoundedBy, + Weighted, F>: LocalAnalysis, { - type FloatType = F; - - fn adjoint_product_bound(&self, seminorm: &ConvolutionOp) -> Option { + fn opnorm_bound(&self, seminorm: &'a ConvolutionOp, _: L2) -> DynResult { // Sensors should not take on negative values to allow // A_*A to be upper bounded by a simple convolution of `spread`. + // TODO: Do we really need this restriction? if self.sensor.bounds().lower() < 0.0 { - return None; + return Err(anyhow!("Sensor not bounded from below by zero")); } // Calculate the factor $L_1$ for betwee $ℱ[ψ * ψ] ≤ L_1 ℱ[ρ]$ for $ψ$ the base spread @@ -517,33 +514,33 @@ let l0 = self.sensor.norm(Linfinity) * self.sensor.norm(L1); // The final transition factor is: - Some(l0 * l1) + Ok((l0 * l1).sqrt()) } } macro_rules! make_sensorgridsupportgenerator_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { impl std::ops::$trait_assign - for SensorGridSupportGenerator + for SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { fn $fn_assign(&mut self, t: F) { self.weights.$fn_assign(t); } } - impl std::ops::$trait for SensorGridSupportGenerator + impl std::ops::$trait for SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { - type Output = SensorGridSupportGenerator; + type Output = SensorGridSupportGenerator; fn $fn(mut self, t: F) -> Self::Output { std::ops::$trait_assign::$fn_assign(&mut self.weights, t); self @@ -551,14 +548,14 @@ } impl<'a, F, S, P, const N: usize> std::ops::$trait - for &'a SensorGridSupportGenerator + for &'a SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { - type Output = SensorGridSupportGenerator; + type Output = SensorGridSupportGenerator; fn $fn(self, t: F) -> Self::Output { SensorGridSupportGenerator { base_sensor: self.base_sensor.clone(), @@ -575,14 +572,14 @@ macro_rules! make_sensorgridsupportgenerator_unaryop { ($trait:ident, $fn:ident) => { - impl std::ops::$trait for SensorGridSupportGenerator + impl std::ops::$trait for SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { - type Output = SensorGridSupportGenerator; + type Output = SensorGridSupportGenerator; fn $fn(mut self) -> Self::Output { self.weights = self.weights.$fn(); self @@ -590,14 +587,14 @@ } impl<'a, F, S, P, const N: usize> std::ops::$trait - for &'a SensorGridSupportGenerator + for &'a SensorGridSupportGenerator where F: Float, - S: Sensor, - P: Spread, - Convolution: Spread, + S: Sensor, + P: Spread, + Convolution: Spread, { - type Output = SensorGridSupportGenerator; + type Output = SensorGridSupportGenerator; fn $fn(self) -> Self::Output { SensorGridSupportGenerator { base_sensor: self.base_sensor.clone(), @@ -612,13 +609,13 @@ make_sensorgridsupportgenerator_unaryop!(Neg, neg); impl<'a, F, S, P, BT, const N: usize> Mapping> - for PreadjointHelper<'a, SensorGrid, RNDM> + for PreadjointHelper<'a, SensorGrid, RNDM> where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, N>, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, N>, { type Codomain = SensorGridBTFN; @@ -634,12 +631,12 @@ } impl<'a, F, S, P, BT, const N: usize> Linear> - for PreadjointHelper<'a, SensorGrid, RNDM> + for PreadjointHelper<'a, SensorGrid, RNDM> where F: Float, BT: SensorGridBT, - S: Sensor, - P: Spread, - Convolution: Spread + LocalAnalysis, N>, + S: Sensor, + P: Spread, + Convolution: Spread + LocalAnalysis, N>, { }