Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model/sensor_grid.rs

src/forward_model/sensor_grid.rs

branch
dev
changeset 61
4f468d35fa29
parent 49
6b0db7251ebe
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
1 /*! 1 /*!
2 Sensor grid forward model 2 Sensor grid forward model
3 */ 3 */
4 4
5 use nalgebra::base::{DMatrix, DVector}; 5 use super::{BasicCurvatureBoundEstimates, ForwardModel};
6 use numeric_literals::replace_float_literals; 6 use crate::frank_wolfe::FindimQuadraticModel;
7 use std::iter::Zip; 7 use crate::kernels::{AutoConvolution, BoundedBy, Convolution};
8 use std::ops::RangeFrom; 8 use crate::measures::{DiscreteMeasure, Radon, RNDM};
9 9 use crate::preadjoint_helper::PreadjointHelper;
10 use crate::seminorms::{ConvolutionOp, SimpleConvolutionKernel};
11 use crate::types::*;
10 use alg_tools::bisection_tree::*; 12 use alg_tools::bisection_tree::*;
11 use alg_tools::error::DynError; 13 use alg_tools::bounds::Bounded;
14 use alg_tools::error::{DynError, DynResult};
15 use alg_tools::euclidean::Euclidean;
12 use alg_tools::instance::Instance; 16 use alg_tools::instance::Instance;
13 use alg_tools::iter::{MapX, Mappable}; 17 use alg_tools::iter::{MapX, Mappable};
14 use alg_tools::lingrid::*; 18 use alg_tools::lingrid::*;
15 pub use alg_tools::linops::*; 19 pub use alg_tools::linops::*;
16 use alg_tools::mapping::{DifferentiableMapping, RealMapping}; 20 use alg_tools::mapping::{DifferentiableMapping, RealMapping};
17 use alg_tools::maputil::map2; 21 use alg_tools::maputil::map2;
18 use alg_tools::nalgebra_support::ToNalgebraRealField; 22 use alg_tools::nalgebra_support::ToNalgebraRealField;
19 use alg_tools::norms::{Linfinity, Norm, L1, L2}; 23 use alg_tools::norms::{Linfinity, Norm, L1, L2};
20 use alg_tools::tabledump::write_csv; 24 use alg_tools::tabledump::write_csv;
21 25 use anyhow::anyhow;
22 use super::{AdjointProductBoundedBy, BoundedCurvature, ForwardModel}; 26 use nalgebra::base::{DMatrix, DVector};
23 use crate::frank_wolfe::FindimQuadraticModel; 27 use numeric_literals::replace_float_literals;
24 use crate::kernels::{AutoConvolution, BoundedBy, Convolution}; 28 use std::iter::Zip;
25 use crate::measures::{DiscreteMeasure, Radon}; 29 use std::ops::RangeFrom;
26 use crate::preadjoint_helper::PreadjointHelper; 30
27 use crate::seminorms::{ConvolutionOp, SimpleConvolutionKernel}; 31 pub type ShiftedSensor<F, S, P, const N: usize> = Shift<Convolution<S, P>, N, F>;
28 use crate::types::*;
29
30 type RNDM<F, const N: usize> = DiscreteMeasure<Loc<F, N>, F>;
31
32 pub type ShiftedSensor<F, S, P, const N: usize> = Shift<Convolution<S, P>, F, N>;
33 32
34 /// Trait for physical convolution models. Has blanket implementation for all cases. 33 /// Trait for physical convolution models. Has blanket implementation for all cases.
35 pub trait Spread<F: Float, const N: usize>: 34 pub trait Spread<const N: usize, F: Float = f64>:
36 'static + Clone + Support<F, N> + RealMapping<F, N> + Bounded<F> 35 'static + Clone + Support<N, F> + RealMapping<N, F> + Bounded<F>
37 { 36 {
38 } 37 }
39 38
40 impl<F, T, const N: usize> Spread<F, N> for T 39 impl<F, T, const N: usize> Spread<N, F> for T
41 where 40 where
42 F: Float, 41 F: Float,
43 T: 'static + Clone + Support<F, N> + Bounded<F> + RealMapping<F, N>, 42 T: 'static + Clone + Support<N, F> + Bounded<F> + RealMapping<N, F>,
44 { 43 {
45 } 44 }
46 45
47 /// Trait for compactly supported sensors. Has blanket implementation for all cases. 46 /// Trait for compactly supported sensors. Has blanket implementation for all cases.
48 pub trait Sensor<F: Float, const N: usize>: 47 pub trait Sensor<const N: usize, F: Float = f64>:
49 Spread<F, N> + Norm<F, L1> + Norm<F, Linfinity> 48 Spread<N, F> + Norm<L1, F> + Norm<Linfinity, F>
50 { 49 {
51 } 50 }
52 51
53 impl<F, T, const N: usize> Sensor<F, N> for T 52 impl<F, T, const N: usize> Sensor<N, F> for T
54 where 53 where
55 F: Float, 54 F: Float,
56 T: Spread<F, N> + Norm<F, L1> + Norm<F, Linfinity>, 55 T: Spread<N, F> + Norm<L1, F> + Norm<Linfinity, F>,
57 { 56 {
58 } 57 }
59 58
60 pub trait SensorGridBT<F, S, P, const N: usize>: 59 pub trait SensorGridBT<F, S, P, const N: usize>:
61 Clone + BTImpl<F, N, Data = usize, Agg = Bounds<F>> 60 Clone + BTImpl<N, F, Data = usize, Agg = Bounds<F>>
62 where 61 where
63 F: Float, 62 F: Float,
64 S: Sensor<F, N>, 63 S: Sensor<N, F>,
65 P: Spread<F, N>, 64 P: Spread<N, F>,
66 { 65 {
67 } 66 }
68 67
69 impl<F, S, P, T, const N: usize> SensorGridBT<F, S, P, N> for T 68 impl<F, S, P, T, const N: usize> SensorGridBT<F, S, P, N> for T
70 where 69 where
71 T: Clone + BTImpl<F, N, Data = usize, Agg = Bounds<F>>, 70 T: Clone + BTImpl<N, F, Data = usize, Agg = Bounds<F>>,
72 F: Float, 71 F: Float,
73 S: Sensor<F, N>, 72 S: Sensor<N, F>,
74 P: Spread<F, N>, 73 P: Spread<N, F>,
75 { 74 {
76 } 75 }
77 76
78 // We need type alias bounds to access associated types 77 // We need type alias bounds to access associated types
79 #[allow(type_alias_bounds)] 78 #[allow(type_alias_bounds)]
80 pub type SensorGridBTFN<F, S, P, BT: SensorGridBT<F, S, P, N>, const N: usize> = 79 pub type SensorGridBTFN<F, S, P, BT: SensorGridBT<F, S, P, N>, const N: usize> =
81 BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>; 80 BTFN<F, SensorGridSupportGenerator<S, F, P, N>, BT, N>;
82 81
83 /// Sensor grid forward model 82 /// Sensor grid forward model
84 #[derive(Clone)] 83 #[derive(Clone)]
85 pub struct SensorGrid<F, S, P, BT, const N: usize> 84 pub struct SensorGrid<F, S, P, BT, const N: usize>
86 where 85 where
87 F: Float, 86 F: Float,
88 S: Sensor<F, N>, 87 S: Sensor<N, F>,
89 P: Spread<F, N>, 88 P: Spread<N, F>,
90 Convolution<S, P>: Spread<F, N>, 89 Convolution<S, P>: Spread<N, F>,
91 BT: SensorGridBT<F, S, P, N>, 90 BT: SensorGridBT<F, S, P, N>,
92 { 91 {
93 domain: Cube<F, N>, 92 domain: Cube<N, F>,
94 sensor_count: [usize; N], 93 sensor_count: [usize; N],
95 sensor: S, 94 sensor: S,
96 spread: P, 95 spread: P,
97 base_sensor: Convolution<S, P>, 96 base_sensor: Convolution<S, P>,
98 bt: BT, 97 bt: BT,
100 99
101 impl<F, S, P, BT, const N: usize> SensorGrid<F, S, P, BT, N> 100 impl<F, S, P, BT, const N: usize> SensorGrid<F, S, P, BT, N>
102 where 101 where
103 F: Float, 102 F: Float,
104 BT: SensorGridBT<F, S, P, N>, 103 BT: SensorGridBT<F, S, P, N>,
105 S: Sensor<F, N>, 104 S: Sensor<N, F>,
106 P: Spread<F, N>, 105 P: Spread<N, F>,
107 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 106 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, BT::Agg, N>,
108 { 107 {
109 /// Create a new sensor grid. 108 /// Create a new sensor grid.
110 /// 109 ///
111 /// The parameter `depth` indicates the search depth of the created [`BT`]s 110 /// The parameter `depth` indicates the search depth of the created [`BT`]s
112 /// for the adjoint values. 111 /// for the adjoint values.
113 pub fn new( 112 pub fn new(
114 domain: Cube<F, N>, 113 domain: Cube<N, F>,
115 sensor_count: [usize; N], 114 sensor_count: [usize; N],
116 sensor: S, 115 sensor: S,
117 spread: P, 116 spread: P,
118 depth: BT::Depth, 117 depth: BT::Depth,
119 ) -> Self { 118 ) -> Self {
139 138
140 impl<F, S, P, BT, const N: usize> SensorGrid<F, S, P, BT, N> 139 impl<F, S, P, BT, const N: usize> SensorGrid<F, S, P, BT, N>
141 where 140 where
142 F: Float, 141 F: Float,
143 BT: SensorGridBT<F, S, P, N>, 142 BT: SensorGridBT<F, S, P, N>,
144 S: Sensor<F, N>, 143 S: Sensor<N, F>,
145 P: Spread<F, N>, 144 P: Spread<N, F>,
146 Convolution<S, P>: Spread<F, N>, 145 Convolution<S, P>: Spread<N, F>,
147 { 146 {
148 /// Return the grid of sensor locations. 147 /// Return the grid of sensor locations.
149 pub fn grid(&self) -> LinGrid<F, N> { 148 pub fn grid(&self) -> LinGrid<N, F> {
150 lingrid_centered(&self.domain, &self.sensor_count) 149 lingrid_centered(&self.domain, &self.sensor_count)
151 } 150 }
152 151
153 /// Returns the number of sensors (number of grid points) 152 /// Returns the number of sensors (number of grid points)
154 pub fn n_sensors(&self) -> usize { 153 pub fn n_sensors(&self) -> usize {
155 self.sensor_count.iter().product() 154 self.sensor_count.iter().product()
156 } 155 }
157 156
158 /// Constructs a sensor shifted by `x`. 157 /// Constructs a sensor shifted by `x`.
159 #[inline] 158 #[inline]
160 fn shifted_sensor(&self, x: Loc<F, N>) -> ShiftedSensor<F, S, P, N> { 159 fn shifted_sensor(&self, x: Loc<N, F>) -> ShiftedSensor<F, S, P, N> {
161 self.base_sensor.clone().shift(x) 160 self.base_sensor.clone().shift(x)
162 } 161 }
163 162
164 #[inline] 163 #[inline]
165 fn _zero_observable(&self) -> DVector<F> { 164 fn _zero_observable(&self) -> DVector<F> {
178 .reduce(F::mul) 177 .reduce(F::mul)
179 .unwrap() 178 .unwrap()
180 } 179 }
181 } 180 }
182 181
183 impl<F, S, P, BT, const N: usize> Mapping<RNDM<F, N>> for SensorGrid<F, S, P, BT, N> 182 impl<F, S, P, BT, const N: usize> Mapping<RNDM<N, F>> for SensorGrid<F, S, P, BT, N>
184 where 183 where
185 F: Float, 184 F: Float,
186 BT: SensorGridBT<F, S, P, N>, 185 BT: SensorGridBT<F, S, P, N>,
187 S: Sensor<F, N>, 186 S: Sensor<N, F>,
188 P: Spread<F, N>, 187 P: Spread<N, F>,
189 Convolution<S, P>: Spread<F, N>, 188 Convolution<S, P>: Spread<N, F>,
190 { 189 {
191 type Codomain = DVector<F>; 190 type Codomain = DVector<F>;
192 191
193 #[inline] 192 #[inline]
194 fn apply<I: Instance<RNDM<F, N>>>(&self, μ: I) -> DVector<F> { 193 fn apply<I: Instance<RNDM<N, F>>>(&self, μ: I) -> DVector<F> {
195 let mut y = self._zero_observable(); 194 let mut y = self._zero_observable();
196 self.apply_add(&mut y, μ); 195 self.apply_add(&mut y, μ);
197 y 196 y
198 } 197 }
199 } 198 }
200 199
201 impl<F, S, P, BT, const N: usize> Linear<RNDM<F, N>> for SensorGrid<F, S, P, BT, N> 200 impl<F, S, P, BT, const N: usize> Linear<RNDM<N, F>> for SensorGrid<F, S, P, BT, N>
202 where 201 where
203 F: Float, 202 F: Float,
204 BT: SensorGridBT<F, S, P, N>, 203 BT: SensorGridBT<F, S, P, N>,
205 S: Sensor<F, N>, 204 S: Sensor<N, F>,
206 P: Spread<F, N>, 205 P: Spread<N, F>,
207 Convolution<S, P>: Spread<F, N>, 206 Convolution<S, P>: Spread<N, F>,
208 { 207 {
209 } 208 }
210 209
211 #[replace_float_literals(F::cast_from(literal))] 210 #[replace_float_literals(F::cast_from(literal))]
212 impl<F, S, P, BT, const N: usize> GEMV<F, RNDM<F, N>, DVector<F>> for SensorGrid<F, S, P, BT, N> 211 impl<F, S, P, BT, const N: usize> GEMV<F, RNDM<N, F>, DVector<F>> for SensorGrid<F, S, P, BT, N>
213 where 212 where
214 F: Float, 213 F: Float,
215 BT: SensorGridBT<F, S, P, N>, 214 BT: SensorGridBT<F, S, P, N>,
216 S: Sensor<F, N>, 215 S: Sensor<N, F>,
217 P: Spread<F, N>, 216 P: Spread<N, F>,
218 Convolution<S, P>: Spread<F, N>, 217 Convolution<S, P>: Spread<N, F>,
219 { 218 {
220 fn gemv<I: Instance<RNDM<F, N>>>(&self, y: &mut DVector<F>, α: F, μ: I, β: F) { 219 fn gemv<I: Instance<RNDM<N, F>>>(&self, y: &mut DVector<F>, α: F, μ: I, β: F) {
221 let grid = self.grid(); 220 let grid = self.grid();
222 if β == 0.0 { 221 if β == 0.0 {
223 y.fill(0.0) 222 y.fill(0.0)
224 } else if β != 1.0 { 223 } else if β != 1.0 {
225 *y *= β; // Need to multiply first, as we have to be able to add to y. 224 *y *= β; // Need to multiply first, as we have to be able to add to y.
226 } 225 }
227 if α == 1.0 { 226 if α == 1.0 {
228 self.apply_add(y, μ) 227 self.apply_add(y, μ)
229 } else { 228 } else {
230 for δ in μ.ref_instance() { 229 μ.eval_ref(|μr| {
230 for δ in μr {
231 for &d in self.bt.iter_at(&δ.x) {
232 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
233 y[d] += sensor.apply(&δ.x) * (α * δ.α);
234 }
235 }
236 })
237 }
238 }
239
240 fn apply_add<I: Instance<RNDM<N, F>>>(&self, y: &mut DVector<F>, μ: I) {
241 let grid = self.grid();
242 μ.eval_ref(|μr| {
243 for δ in μr {
231 for &d in self.bt.iter_at(&δ.x) { 244 for &d in self.bt.iter_at(&δ.x) {
232 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d)); 245 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
233 y[d] += sensor.apply(&δ.x) * (α * δ.α); 246 y[d] += sensor.apply(&δ.x) * δ.α;
234 } 247 }
235 } 248 }
236 } 249 })
237 } 250 }
238 251 }
239 fn apply_add<I: Instance<RNDM<F, N>>>(&self, y: &mut DVector<F>, μ: I) { 252
240 let grid = self.grid(); 253 impl<F, S, P, BT, const N: usize> BoundedLinear<RNDM<N, F>, Radon, L2, F>
241 for δ in μ.ref_instance() {
242 for &d in self.bt.iter_at(&δ.x) {
243 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
244 y[d] += sensor.apply(&δ.x) * δ.α;
245 }
246 }
247 }
248 }
249
250 impl<F, S, P, BT, const N: usize> BoundedLinear<RNDM<F, N>, Radon, L2, F>
251 for SensorGrid<F, S, P, BT, N> 254 for SensorGrid<F, S, P, BT, N>
252 where 255 where
253 F: Float, 256 F: Float,
254 BT: SensorGridBT<F, S, P, N, Agg = Bounds<F>>, 257 BT: SensorGridBT<F, S, P, N>,
255 S: Sensor<F, N>, 258 S: Sensor<N, F>,
256 P: Spread<F, N>, 259 P: Spread<N, F>,
257 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 260 Convolution<S, P>: Spread<N, F>,
258 { 261 {
259 /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped 262 /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped
260 /// with the Radon norm, and $ℝ^n$ with the Euclidean norm. 263 /// with the Radon norm, and $ℝ^n$ with the Euclidean norm.
261 fn opnorm_bound(&self, _: Radon, _: L2) -> F { 264 fn opnorm_bound(&self, _: Radon, _: L2) -> DynResult<F> {
262 // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have 265 // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have
263 // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩ 266 // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩
264 // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ 267 // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ
265 // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ 268 // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ
266 // ≤ sup_{|z|_2 ≤ 1} |φ(y)| ∑_{i:th sensor active at y}|z_i| |μ|_ℳ 269 // ≤ sup_{|z|_2 ≤ 1} |φ(y)| ∑_{i:th sensor active at y}|z_i| |μ|_ℳ
269 // where N_ψ is the maximum number of sensors that overlap, and 272 // where N_ψ is the maximum number of sensors that overlap, and
270 // |z|_2 is restricted to the active sensors. 273 // |z|_2 is restricted to the active sensors.
271 // = |φ|_∞ √N_ψ |μ|_ℳ. 274 // = |φ|_∞ √N_ψ |μ|_ℳ.
272 // Hence 275 // Hence
273 let n = self.max_overlapping(); 276 let n = self.max_overlapping();
274 self.base_sensor.bounds().uniform() * n.sqrt() 277 Ok(self.base_sensor.bounds().uniform() * n.sqrt())
275 } 278 }
276 } 279 }
277 280
278 type SensorGridPreadjoint<'a, A, F, const N: usize> = PreadjointHelper<'a, A, RNDM<F, N>>; 281 type SensorGridPreadjoint<'a, A, F, const N: usize> = PreadjointHelper<'a, A, RNDM<N, F>>;
279 282
280 impl<F, S, P, BT, const N: usize> Preadjointable<RNDM<F, N>, DVector<F>> 283 impl<F, S, P, BT, const N: usize> Preadjointable<RNDM<N, F>, DVector<F>>
281 for SensorGrid<F, S, P, BT, N> 284 for SensorGrid<F, S, P, BT, N>
282 where 285 where
283 F: Float, 286 F: Float,
284 BT: SensorGridBT<F, S, P, N>, 287 BT: SensorGridBT<F, S, P, N>,
285 S: Sensor<F, N>, 288 S: Sensor<N, F>,
286 P: Spread<F, N>, 289 P: Spread<N, F>,
287 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 290 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, BT::Agg, N>,
288 { 291 {
289 type PreadjointCodomain = BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>; 292 type PreadjointCodomain = BTFN<F, SensorGridSupportGenerator<S, F, P, N>, BT, N>;
290 type Preadjoint<'a> 293 type Preadjoint<'a>
291 = SensorGridPreadjoint<'a, Self, F, N> 294 = SensorGridPreadjoint<'a, Self, F, N>
292 where 295 where
293 Self: 'a; 296 Self: 'a;
294 297
301 #[replace_float_literals(F::cast_from(literal))] 304 #[replace_float_literals(F::cast_from(literal))]
302 impl<'a, F, S, P, BT, const N : usize> LipschitzValues 305 impl<'a, F, S, P, BT, const N : usize> LipschitzValues
303 for SensorGridPreadjoint<'a, SensorGrid<F, S, P, BT, N>, F, N> 306 for SensorGridPreadjoint<'a, SensorGrid<F, S, P, BT, N>, F, N>
304 where F : Float, 307 where F : Float,
305 BT : SensorGridBT<F, S, P, N>, 308 BT : SensorGridBT<F, S, P, N>,
306 S : Sensor<F, N>, 309 S : Sensor<N, F>,
307 P : Spread<F, N>, 310 P : Spread<N, F>,
308 Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType=F> + DifferentiableMapping<Loc<F,N>> + LocalAnalysis<F, BT::Agg, N>, 311 Convolution<S, P> : Spread<N, F> + Lipschitz<L2, FloatType=F> + DifferentiableMapping<Loc<N, F>> + LocalAnalysis<F, BT::Agg, N>,
309 for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<F,N>>>::Differential<'b> : Lipschitz<L2, FloatType=F>, 312 for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<N, F>>>::Differential<'b> : Lipschitz<L2, FloatType=F>,
310 { 313 {
311 314
312 type FloatType = F; 315 type FloatType = F;
313 316
314 fn value_unit_lipschitz_factor(&self) -> Option<F> { 317 fn value_unit_lipschitz_factor(&self) -> Option<F> {
330 } 333 }
331 } 334 }
332 */ 335 */
333 336
334 #[replace_float_literals(F::cast_from(literal))] 337 #[replace_float_literals(F::cast_from(literal))]
335 impl<'a, F, S, P, BT, const N: usize> BoundedCurvature for SensorGrid<F, S, P, BT, N> 338 impl<'a, F, S, P, BT, const N: usize> BasicCurvatureBoundEstimates<F> for SensorGrid<F, S, P, BT, N>
336 where 339 where
337 F: Float, 340 F: Float + ToNalgebraRealField,
338 BT: SensorGridBT<F, S, P, N>, 341 BT: SensorGridBT<F, S, P, N>,
339 S: Sensor<F, N>, 342 S: Sensor<N, F>,
340 P: Spread<F, N>, 343 P: Spread<N, F>,
341 Convolution<S, P>: Spread<F, N> 344 DVector<F>: Euclidean<F>,
345 Convolution<S, P>: Spread<N, F>
342 + Lipschitz<L2, FloatType = F> 346 + Lipschitz<L2, FloatType = F>
343 + DifferentiableMapping<Loc<F, N>> 347 + DifferentiableMapping<Loc<N, F>>
344 + LocalAnalysis<F, BT::Agg, N>, 348 + LocalAnalysis<F, BT::Agg, N>,
345 for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<F, N>>>::Differential<'b>: 349 for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<N, F>>>::Differential<'b>:
346 Lipschitz<L2, FloatType = F>, 350 Lipschitz<L2, FloatType = F>,
347 { 351 {
348 type FloatType = F; 352 fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>) {
349
350 /// Returns factors $ℓ_F$ and $Θ²$ such that
351 /// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ)+F'(μ+Δ)|Δ⟩ ≤ Θ²|γ|(c_2)‖γ‖$,
352 /// where $Δ=(π_♯^1-π_♯^0)γ$.
353 ///
354 /// See Lemma 3.8, Lemma 5.10, Remark 5.14, and Example 5.15.
355 fn curvature_bound_components(&self) -> (Option<Self::FloatType>, Option<Self::FloatType>) {
356 let n_ψ = self.max_overlapping(); 353 let n_ψ = self.max_overlapping();
357 let ψ_diff_lip = self.base_sensor.diff_ref().lipschitz_factor(L2); 354 let ψ_diff_lip = self.base_sensor.diff_ref().lipschitz_factor(L2);
358 let ψ_lip = self.base_sensor.lipschitz_factor(L2); 355 let ψ_lip = self.base_sensor.lipschitz_factor(L2);
359 let ℓ_F = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l); 356 let ℓ_F0 = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l);
360 let θ2 = ψ_lip.map(|l| 4.0 * n_ψ * l.powi(2)); 357 let Θ2 = ψ_lip.map(|l| 4.0 * n_ψ * l.powi(2));
361 358
362 (ℓ_F, θ2) 359 (ℓ_F0, Θ2)
363 } 360 }
364 } 361 }
365 362
366 #[derive(Clone, Debug)] 363 #[derive(Clone, Debug)]
367 pub struct SensorGridSupportGenerator<F, S, P, const N: usize> 364 pub struct SensorGridSupportGenerator<S, F, P, const N: usize>
368 where 365 where
369 F: Float, 366 F: Float,
370 S: Sensor<F, N>, 367 S: Sensor<N, F>,
371 P: Spread<F, N>, 368 P: Spread<N, F>,
372 { 369 {
373 base_sensor: Convolution<S, P>, 370 base_sensor: Convolution<S, P>,
374 grid: LinGrid<F, N>, 371 grid: LinGrid<N, F>,
375 weights: DVector<F>, 372 weights: DVector<F>,
376 } 373 }
377 374
378 impl<F, S, P, const N: usize> SensorGridSupportGenerator<F, S, P, N> 375 impl<F, S, P, const N: usize> SensorGridSupportGenerator<S, F, P, N>
379 where 376 where
380 F: Float, 377 F: Float,
381 S: Sensor<F, N>, 378 S: Sensor<N, F>,
382 P: Spread<F, N>, 379 P: Spread<N, F>,
383 Convolution<S, P>: Spread<F, N>, 380 Convolution<S, P>: Spread<N, F>,
384 { 381 {
385 #[inline] 382 #[inline]
386 fn construct_sensor(&self, id: usize, w: F) -> Weighted<ShiftedSensor<F, S, P, N>, F> { 383 fn construct_sensor(&self, id: usize, w: F) -> Weighted<ShiftedSensor<F, S, P, N>, F> {
387 let x = self.grid.entry_linear_unchecked(id); 384 let x = self.grid.entry_linear_unchecked(id);
388 self.base_sensor.clone().shift(x).weigh(w) 385 self.base_sensor.clone().shift(x).weigh(w)
395 ) -> (usize, Weighted<ShiftedSensor<F, S, P, N>, F>) { 392 ) -> (usize, Weighted<ShiftedSensor<F, S, P, N>, F>) {
396 (id.into(), self.construct_sensor(id, *w)) 393 (id.into(), self.construct_sensor(id, *w))
397 } 394 }
398 } 395 }
399 396
400 impl<F, S, P, const N: usize> SupportGenerator<F, N> for SensorGridSupportGenerator<F, S, P, N> 397 impl<F, S, P, const N: usize> SupportGenerator<N, F> for SensorGridSupportGenerator<S, F, P, N>
401 where 398 where
402 F: Float, 399 F: Float,
403 S: Sensor<F, N>, 400 S: Sensor<N, F>,
404 P: Spread<F, N>, 401 P: Spread<N, F>,
405 Convolution<S, P>: Spread<F, N>, 402 Convolution<S, P>: Spread<N, F>,
406 { 403 {
407 type Id = usize; 404 type Id = usize;
408 type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>; 405 type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>;
409 type AllDataIter<'a> 406 type AllDataIter<'a>
410 = MapX< 407 = MapX<
432 .zip(self.weights.as_slice().iter()) 429 .zip(self.weights.as_slice().iter())
433 .mapX(self, Self::construct_sensor_and_id) 430 .mapX(self, Self::construct_sensor_and_id)
434 } 431 }
435 } 432 }
436 433
437 impl<F, S, P, BT, const N: usize> ForwardModel<DiscreteMeasure<Loc<F, N>, F>, F> 434 impl<F, S, P, BT, const N: usize> ForwardModel<DiscreteMeasure<Loc<N, F>, F>, F>
438 for SensorGrid<F, S, P, BT, N> 435 for SensorGrid<F, S, P, BT, N>
439 where 436 where
440 F: Float + ToNalgebraRealField<MixedType = F> + nalgebra::RealField, 437 F: Float + ToNalgebraRealField<MixedType = F> + nalgebra::RealField,
441 BT: SensorGridBT<F, S, P, N>, 438 BT: SensorGridBT<F, S, P, N>,
442 S: Sensor<F, N>, 439 S: Sensor<N, F>,
443 P: Spread<F, N>, 440 P: Spread<N, F>,
444 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 441 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, BT::Agg, N>,
445 { 442 {
446 type Observable = DVector<F>; 443 type Observable = DVector<F>;
447 444
448 fn write_observable(&self, b: &Self::Observable, prefix: String) -> DynError { 445 fn write_observable(&self, b: &Self::Observable, prefix: String) -> DynError {
449 let it = self.grid().into_iter().zip(b.iter()).map(|(x, &v)| (x, v)); 446 let it = self.grid().into_iter().zip(b.iter()).map(|(x, &v)| (x, v));
454 fn zero_observable(&self) -> Self::Observable { 451 fn zero_observable(&self) -> Self::Observable {
455 self._zero_observable() 452 self._zero_observable()
456 } 453 }
457 } 454 }
458 455
459 impl<F, S, P, BT, const N: usize> FindimQuadraticModel<Loc<F, N>, F> for SensorGrid<F, S, P, BT, N> 456 impl<F, S, P, BT, const N: usize> FindimQuadraticModel<Loc<N, F>, F> for SensorGrid<F, S, P, BT, N>
460 where 457 where
461 F: Float + ToNalgebraRealField<MixedType = F> + nalgebra::RealField, 458 F: Float + ToNalgebraRealField<MixedType = F> + nalgebra::RealField,
462 BT: SensorGridBT<F, S, P, N>, 459 BT: SensorGridBT<F, S, P, N>,
463 S: Sensor<F, N>, 460 S: Sensor<N, F>,
464 P: Spread<F, N>, 461 P: Spread<N, F>,
465 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 462 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, BT::Agg, N>,
466 { 463 {
467 fn findim_quadratic_model( 464 fn findim_quadratic_model(
468 &self, 465 &self,
469 μ: &DiscreteMeasure<Loc<F, N>, F>, 466 μ: &DiscreteMeasure<Loc<N, F>, F>,
470 b: &Self::Observable, 467 b: &Self::Observable,
471 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>) { 468 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>) {
472 assert_eq!(b.len(), self.n_sensors()); 469 assert_eq!(b.len(), self.n_sensors());
473 let mut mA = DMatrix::zeros(self.n_sensors(), μ.len()); 470 let mut mA = DMatrix::zeros(self.n_sensors(), μ.len());
474 let grid = self.grid(); 471 let grid = self.grid();
481 let mAt = mA.transpose(); 478 let mAt = mA.transpose();
482 (&mAt * mA, &mAt * b) 479 (&mAt * mA, &mAt * b)
483 } 480 }
484 } 481 }
485 482
486 /// Implements the calculation a factor $L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model 483 /// Implements the calculation a factor $√L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model
487 /// and $𝒟$ a seminorm of suitable form. 484 /// and $𝒟$ a seminorm of suitable form.
488 /// 485 ///
489 /// **This assumes (but does not check) that the sensors are not overlapping.** 486 /// **This assumes (but does not check) that the sensors are not overlapping.**
490 #[replace_float_literals(F::cast_from(literal))] 487 #[replace_float_literals(F::cast_from(literal))]
491 impl<F, BT, S, P, K, const N: usize> AdjointProductBoundedBy<RNDM<F, N>, ConvolutionOp<F, K, BT, N>> 488 impl<'a, F, BT, S, P, K, const N: usize>
492 for SensorGrid<F, S, P, BT, N> 489 BoundedLinear<RNDM<N, F>, &'a ConvolutionOp<F, K, BT, N>, L2, F> for SensorGrid<F, S, P, BT, N>
493 where 490 where
494 F: Float + nalgebra::RealField + ToNalgebraRealField, 491 F: Float + ToNalgebraRealField,
495 BT: SensorGridBT<F, S, P, N>, 492 BT: SensorGridBT<F, S, P, N>,
496 S: Sensor<F, N>, 493 S: Sensor<N, F>,
497 P: Spread<F, N>, 494 P: Spread<N, F>,
498 Convolution<S, P>: Spread<F, N>, 495 Convolution<S, P>: Spread<N, F>,
499 K: SimpleConvolutionKernel<F, N>, 496 K: SimpleConvolutionKernel<N, F>,
500 AutoConvolution<P>: BoundedBy<F, K>, 497 AutoConvolution<P>: BoundedBy<F, K>,
501 { 498 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
502 type FloatType = F; 499 {
503 500 fn opnorm_bound(&self, seminorm: &'a ConvolutionOp<F, K, BT, N>, _: L2) -> DynResult<F> {
504 fn adjoint_product_bound(&self, seminorm: &ConvolutionOp<F, K, BT, N>) -> Option<F> {
505 // Sensors should not take on negative values to allow 501 // Sensors should not take on negative values to allow
506 // A_*A to be upper bounded by a simple convolution of `spread`. 502 // A_*A to be upper bounded by a simple convolution of `spread`.
503 // TODO: Do we really need this restriction?
507 if self.sensor.bounds().lower() < 0.0 { 504 if self.sensor.bounds().lower() < 0.0 {
508 return None; 505 return Err(anyhow!("Sensor not bounded from below by zero"));
509 } 506 }
510 507
511 // Calculate the factor $L_1$ for betwee $ℱ[ψ * ψ] ≤ L_1 ℱ[ρ]$ for $ψ$ the base spread 508 // Calculate the factor $L_1$ for betwee $ℱ[ψ * ψ] ≤ L_1 ℱ[ρ]$ for $ψ$ the base spread
512 // and $ρ$ the kernel of the seminorm. 509 // and $ρ$ the kernel of the seminorm.
513 let l1 = AutoConvolution(self.spread.clone()).bounding_factor(seminorm.kernel())?; 510 let l1 = AutoConvolution(self.spread.clone()).bounding_factor(seminorm.kernel())?;
515 // Calculate the factor for transitioning from $A_*A$ to `AutoConvolution<P>`, where A 512 // Calculate the factor for transitioning from $A_*A$ to `AutoConvolution<P>`, where A
516 // consists of several `Convolution<S, P>` for the physical model `P` and the sensor `S`. 513 // consists of several `Convolution<S, P>` for the physical model `P` and the sensor `S`.
517 let l0 = self.sensor.norm(Linfinity) * self.sensor.norm(L1); 514 let l0 = self.sensor.norm(Linfinity) * self.sensor.norm(L1);
518 515
519 // The final transition factor is: 516 // The final transition factor is:
520 Some(l0 * l1) 517 Ok((l0 * l1).sqrt())
521 } 518 }
522 } 519 }
523 520
524 macro_rules! make_sensorgridsupportgenerator_scalarop_rhs { 521 macro_rules! make_sensorgridsupportgenerator_scalarop_rhs {
525 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { 522 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
526 impl<F, S, P, const N: usize> std::ops::$trait_assign<F> 523 impl<F, S, P, const N: usize> std::ops::$trait_assign<F>
527 for SensorGridSupportGenerator<F, S, P, N> 524 for SensorGridSupportGenerator<S, F, P, N>
528 where 525 where
529 F: Float, 526 F: Float,
530 S: Sensor<F, N>, 527 S: Sensor<N, F>,
531 P: Spread<F, N>, 528 P: Spread<N, F>,
532 Convolution<S, P>: Spread<F, N>, 529 Convolution<S, P>: Spread<N, F>,
533 { 530 {
534 fn $fn_assign(&mut self, t: F) { 531 fn $fn_assign(&mut self, t: F) {
535 self.weights.$fn_assign(t); 532 self.weights.$fn_assign(t);
536 } 533 }
537 } 534 }
538 535
539 impl<F, S, P, const N: usize> std::ops::$trait<F> for SensorGridSupportGenerator<F, S, P, N> 536 impl<F, S, P, const N: usize> std::ops::$trait<F> for SensorGridSupportGenerator<S, F, P, N>
540 where 537 where
541 F: Float, 538 F: Float,
542 S: Sensor<F, N>, 539 S: Sensor<N, F>,
543 P: Spread<F, N>, 540 P: Spread<N, F>,
544 Convolution<S, P>: Spread<F, N>, 541 Convolution<S, P>: Spread<N, F>,
545 { 542 {
546 type Output = SensorGridSupportGenerator<F, S, P, N>; 543 type Output = SensorGridSupportGenerator<S, F, P, N>;
547 fn $fn(mut self, t: F) -> Self::Output { 544 fn $fn(mut self, t: F) -> Self::Output {
548 std::ops::$trait_assign::$fn_assign(&mut self.weights, t); 545 std::ops::$trait_assign::$fn_assign(&mut self.weights, t);
549 self 546 self
550 } 547 }
551 } 548 }
552 549
553 impl<'a, F, S, P, const N: usize> std::ops::$trait<F> 550 impl<'a, F, S, P, const N: usize> std::ops::$trait<F>
554 for &'a SensorGridSupportGenerator<F, S, P, N> 551 for &'a SensorGridSupportGenerator<S, F, P, N>
555 where 552 where
556 F: Float, 553 F: Float,
557 S: Sensor<F, N>, 554 S: Sensor<N, F>,
558 P: Spread<F, N>, 555 P: Spread<N, F>,
559 Convolution<S, P>: Spread<F, N>, 556 Convolution<S, P>: Spread<N, F>,
560 { 557 {
561 type Output = SensorGridSupportGenerator<F, S, P, N>; 558 type Output = SensorGridSupportGenerator<S, F, P, N>;
562 fn $fn(self, t: F) -> Self::Output { 559 fn $fn(self, t: F) -> Self::Output {
563 SensorGridSupportGenerator { 560 SensorGridSupportGenerator {
564 base_sensor: self.base_sensor.clone(), 561 base_sensor: self.base_sensor.clone(),
565 grid: self.grid, 562 grid: self.grid,
566 weights: (&self.weights).$fn(t), 563 weights: (&self.weights).$fn(t),
573 make_sensorgridsupportgenerator_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); 570 make_sensorgridsupportgenerator_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
574 make_sensorgridsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign); 571 make_sensorgridsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign);
575 572
576 macro_rules! make_sensorgridsupportgenerator_unaryop { 573 macro_rules! make_sensorgridsupportgenerator_unaryop {
577 ($trait:ident, $fn:ident) => { 574 ($trait:ident, $fn:ident) => {
578 impl<F, S, P, const N: usize> std::ops::$trait for SensorGridSupportGenerator<F, S, P, N> 575 impl<F, S, P, const N: usize> std::ops::$trait for SensorGridSupportGenerator<S, F, P, N>
579 where 576 where
580 F: Float, 577 F: Float,
581 S: Sensor<F, N>, 578 S: Sensor<N, F>,
582 P: Spread<F, N>, 579 P: Spread<N, F>,
583 Convolution<S, P>: Spread<F, N>, 580 Convolution<S, P>: Spread<N, F>,
584 { 581 {
585 type Output = SensorGridSupportGenerator<F, S, P, N>; 582 type Output = SensorGridSupportGenerator<S, F, P, N>;
586 fn $fn(mut self) -> Self::Output { 583 fn $fn(mut self) -> Self::Output {
587 self.weights = self.weights.$fn(); 584 self.weights = self.weights.$fn();
588 self 585 self
589 } 586 }
590 } 587 }
591 588
592 impl<'a, F, S, P, const N: usize> std::ops::$trait 589 impl<'a, F, S, P, const N: usize> std::ops::$trait
593 for &'a SensorGridSupportGenerator<F, S, P, N> 590 for &'a SensorGridSupportGenerator<S, F, P, N>
594 where 591 where
595 F: Float, 592 F: Float,
596 S: Sensor<F, N>, 593 S: Sensor<N, F>,
597 P: Spread<F, N>, 594 P: Spread<N, F>,
598 Convolution<S, P>: Spread<F, N>, 595 Convolution<S, P>: Spread<N, F>,
599 { 596 {
600 type Output = SensorGridSupportGenerator<F, S, P, N>; 597 type Output = SensorGridSupportGenerator<S, F, P, N>;
601 fn $fn(self) -> Self::Output { 598 fn $fn(self) -> Self::Output {
602 SensorGridSupportGenerator { 599 SensorGridSupportGenerator {
603 base_sensor: self.base_sensor.clone(), 600 base_sensor: self.base_sensor.clone(),
604 grid: self.grid, 601 grid: self.grid,
605 weights: (&self.weights).$fn(), 602 weights: (&self.weights).$fn(),
610 } 607 }
611 608
612 make_sensorgridsupportgenerator_unaryop!(Neg, neg); 609 make_sensorgridsupportgenerator_unaryop!(Neg, neg);
613 610
614 impl<'a, F, S, P, BT, const N: usize> Mapping<DVector<F>> 611 impl<'a, F, S, P, BT, const N: usize> Mapping<DVector<F>>
615 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<F, N>> 612 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<N, F>>
616 where 613 where
617 F: Float, 614 F: Float,
618 BT: SensorGridBT<F, S, P, N>, 615 BT: SensorGridBT<F, S, P, N>,
619 S: Sensor<F, N>, 616 S: Sensor<N, F>,
620 P: Spread<F, N>, 617 P: Spread<N, F>,
621 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, 618 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, Bounds<F>, N>,
622 { 619 {
623 type Codomain = SensorGridBTFN<F, S, P, BT, N>; 620 type Codomain = SensorGridBTFN<F, S, P, BT, N>;
624 621
625 fn apply<I: Instance<DVector<F>>>(&self, x: I) -> Self::Codomain { 622 fn apply<I: Instance<DVector<F>>>(&self, x: I) -> Self::Codomain {
626 let fwd = &self.forward_op; 623 let fwd = &self.forward_op;
632 BTFN::new_refresh(&fwd.bt, generator) 629 BTFN::new_refresh(&fwd.bt, generator)
633 } 630 }
634 } 631 }
635 632
636 impl<'a, F, S, P, BT, const N: usize> Linear<DVector<F>> 633 impl<'a, F, S, P, BT, const N: usize> Linear<DVector<F>>
637 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<F, N>> 634 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<N, F>>
638 where 635 where
639 F: Float, 636 F: Float,
640 BT: SensorGridBT<F, S, P, N>, 637 BT: SensorGridBT<F, S, P, N>,
641 S: Sensor<F, N>, 638 S: Sensor<N, F>,
642 P: Spread<F, N>, 639 P: Spread<N, F>,
643 Convolution<S, P>: Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, 640 Convolution<S, P>: Spread<N, F> + LocalAnalysis<F, Bounds<F>, N>,
644 { 641 {
645 } 642 }

Mercurial Repository: pointsource_algs

mercurial