Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model.rs

src/forward_model.rs

branch
dev
changeset 35
b087e3eab191
parent 34
efa60bc4f743
child 44
03251c546744
equal deleted inserted replaced
34:efa60bc4f743 35:b087e3eab191
1 /*! 1 /*!
2 Forward models from discrete measures to observations. 2 Forward models from discrete measures to observations.
3 */ 3 */
4 4
5 use numeric_literals::replace_float_literals;
6 use nalgebra::base::{
7 DMatrix,
8 DVector
9 };
10 use std::iter::Zip;
11 use std::ops::RangeFrom;
12 use std::marker::PhantomData;
13
14 pub use alg_tools::linops::*; 5 pub use alg_tools::linops::*;
15 use alg_tools::euclidean::Euclidean; 6 use alg_tools::euclidean::Euclidean;
16 use alg_tools::norms::{
17 L1, Linfinity, L2, Norm
18 };
19 use alg_tools::bisection_tree::*;
20 use alg_tools::mapping::RealMapping;
21 use alg_tools::lingrid::*;
22 use alg_tools::iter::{MapX, Mappable};
23 use alg_tools::nalgebra_support::ToNalgebraRealField;
24 use alg_tools::tabledump::write_csv;
25 use alg_tools::error::DynError; 7 use alg_tools::error::DynError;
26 use alg_tools::maputil::map2; 8 use alg_tools::instance::Instance;
9 use alg_tools::norms::{NormExponent, L2, Norm};
27 10
28 use crate::types::*; 11 use crate::types::*;
29 use crate::measures::*; 12 use crate::measures::Radon;
30 use crate::seminorms::{ 13 pub mod sensor_grid;
31 ConvolutionOp, 14 pub mod bias;
32 SimpleConvolutionKernel,
33 };
34 use crate::kernels::{
35 Convolution,
36 AutoConvolution,
37 BoundedBy,
38 };
39 use crate::types::L2Squared;
40 use crate::transport::TransportLipschitz;
41
42 pub type RNDM<F, const N : usize> = DiscreteMeasure<Loc<F,N>, F>;
43 15
44 /// `ForwardeModel`s are bounded preadjointable linear operators $A ∈ 𝕃(𝒵(Ω); E)$ 16 /// `ForwardeModel`s are bounded preadjointable linear operators $A ∈ 𝕃(𝒵(Ω); E)$
45 /// where $𝒵(Ω) ⊂ ℳ(Ω)$ is the space of sums of delta measures, presented by 17 /// where $𝒵(Ω) ⊂ ℳ(Ω)$ is the space of sums of delta measures, presented by
46 /// [`DiscreteMeasure`], and $E$ is a [`Euclidean`] space. 18 /// [`crate::measures::DiscreteMeasure`], and $E$ is a [`Euclidean`] space.
47 pub trait ForwardModel<Domain, F : Float + ToNalgebraRealField> 19 pub trait ForwardModel<Domain : Space, F : Float = f64, E : NormExponent = Radon>
48 : BoundedLinear<DiscreteMeasure<Domain, F>, Codomain=Self::Observable, FloatType=F> 20 : BoundedLinear<Domain, E, L2, F, Codomain=Self::Observable>
49 + GEMV<F, DiscreteMeasure<Domain, F>, Self::Observable> 21 + GEMV<F, Domain, Self::Observable>
50 + Linear<DeltaMeasure<Domain, F>, Codomain=Self::Observable> 22 + Preadjointable<Domain, Self::Observable>
51 + Preadjointable<DiscreteMeasure<Domain, F>, Self::Observable> { 23 where
24 for<'a> Self::Observable : Instance<Self::Observable>,
25 Domain : Norm<F, E>,
26 {
52 /// The codomain or value space (of “observables”) for this operator. 27 /// The codomain or value space (of “observables”) for this operator.
53 /// It is assumed to be a [`Euclidean`] space, and therefore also (identified with) 28 /// It is assumed to be a [`Euclidean`] space, and therefore also (identified with)
54 /// the domain of the preadjoint. 29 /// the domain of the preadjoint.
55 type Observable : Euclidean<F, Output=Self::Observable> 30 type Observable : Euclidean<F, Output=Self::Observable>
56 + AXPY<F> 31 + AXPY<F>
32 + Space
57 + Clone; 33 + Clone;
58
59 /// Return A_*A and A_* b
60 fn findim_quadratic_model(
61 &self,
62 μ : &DiscreteMeasure<Domain, F>,
63 b : &Self::Observable
64 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>);
65 34
66 /// Write an observable into a file. 35 /// Write an observable into a file.
67 fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError; 36 fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError;
68 37
69 /// Returns a zero observable 38 /// Returns a zero observable
70 fn zero_observable(&self) -> Self::Observable; 39 fn zero_observable(&self) -> Self::Observable;
71
72 /// Returns an empty (uninitialised) observable.
73 ///
74 /// This is used as a placeholder for temporary [`std::mem::replace`] move operations.
75 fn empty_observable(&self) -> Self::Observable;
76 } 40 }
77 41
78 pub type ShiftedSensor<F, S, P, const N : usize> = Shift<Convolution<S, P>, F, N>; 42 /// Trait for operators $A$ for which $A_*A$ is bounded by some other operator.
79 43 pub trait AdjointProductBoundedBy<Domain : Space, D> : Linear<Domain> {
80 /// Trait for physical convolution models. Has blanket implementation for all cases. 44 type FloatType : Float;
81 pub trait Spread<F : Float, const N : usize> 45 /// Return $L$ such that $A_*A ≤ LD$.
82 : 'static + Clone + Support<F, N> + RealMapping<F, N> + Bounded<F> {} 46 fn adjoint_product_bound(&self, other : &D) -> Option<Self::FloatType>;
83
84 impl<F, T, const N : usize> Spread<F, N> for T
85 where F : Float,
86 T : 'static + Clone + Support<F, N> + Bounded<F> + RealMapping<F, N> {}
87
88 /// Trait for compactly supported sensors. Has blanket implementation for all cases.
89 pub trait Sensor<F : Float, const N : usize> : Spread<F, N> + Norm<F, L1> + Norm<F, Linfinity> {}
90
91 impl<F, T, const N : usize> Sensor<F, N> for T
92 where F : Float,
93 T : Spread<F, N> + Norm<F, L1> + Norm<F, Linfinity> {}
94
95
96 pub trait SensorGridBT<F, S, P, const N : usize> :
97 Clone + BTImpl<F, N, Data=usize, Agg=Bounds<F>>
98 where F : Float,
99 S : Sensor<F, N>,
100 P : Spread<F, N> {}
101
102 impl<F, S, P, T, const N : usize>
103 SensorGridBT<F, S, P, N>
104 for T
105 where T : Clone + BTImpl<F, N, Data=usize, Agg=Bounds<F>>,
106 F : Float,
107 S : Sensor<F, N>,
108 P : Spread<F, N> {}
109
110 // We need type alias bounds to access associated types
111 #[allow(type_alias_bounds)]
112 type SensorGridBTFN<F, S, P, BT : SensorGridBT<F, S, P, N>, const N : usize>
113 = BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>;
114
115 /// Sensor grid forward model
116 #[derive(Clone)]
117 pub struct SensorGrid<F, S, P, BT, const N : usize>
118 where F : Float,
119 S : Sensor<F, N>,
120 P : Spread<F, N>,
121 Convolution<S, P> : Spread<F, N>,
122 BT : SensorGridBT<F, S, P, N>, {
123 domain : Cube<F, N>,
124 sensor_count : [usize; N],
125 sensor : S,
126 spread : P,
127 base_sensor : Convolution<S, P>,
128 bt : BT,
129 } 47 }
130 48
131 impl<F, S, P, BT, const N : usize> SensorGrid<F, S, P, BT, N> 49 /// Trait for operators $A$ for which $A_*A$ is bounded by a diagonal operator.
132 where F : Float, 50 pub trait AdjointProductPairBoundedBy<Domain : Space, D1, D2> : Linear<Domain> {
133 BT : SensorGridBT<F, S, P, N>, 51 type FloatType : Float;
134 S : Sensor<F, N>, 52 /// Return $(L, L_z)$ such that $A_*A ≤ (L_1 D_1, L_2 D_2)$.
135 P : Spread<F, N>, 53 fn adjoint_product_pair_bound(&self, other1 : &D1, other_2 : &D2)
136 Convolution<S, P> : Spread<F, N>, 54 -> Option<(Self::FloatType, Self::FloatType)>;
137 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> { 55 }
138 56
139 pub fn new( 57 /// Trait for [`ForwardModel`]s whose preadjoint has Lipschitz values.
140 domain : Cube<F, N>, 58 pub trait LipschitzValues {
141 sensor_count : [usize; N], 59 type FloatType : Float;
142 sensor : S, 60 /// Return (if one exists) a factor $L$ such that $A_*z$ is $L$-Lipschitz for all
143 spread : P, 61 /// $z$ in the unit ball.
144 depth : BT::Depth 62 fn value_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
145 ) -> Self { 63 None
146 let base_sensor = Convolution(sensor.clone(), spread.clone());
147 let bt = BT::new(domain, depth);
148 let mut sensorgrid = SensorGrid {
149 domain,
150 sensor_count,
151 sensor,
152 spread,
153 base_sensor,
154 bt,
155 };
156
157 for (x, id) in sensorgrid.grid().into_iter().zip(0usize..) {
158 let s = sensorgrid.shifted_sensor(x);
159 sensorgrid.bt.insert(id, &s);
160 }
161
162 sensorgrid
163 } 64 }
164 65
165 pub fn grid(&self) -> LinGrid<F, N> { 66 /// Return (if one exists) a factor $L$ such that $∇A_*z$ is $L$-Lipschitz for all
166 lingrid_centered(&self.domain, &self.sensor_count) 67 /// $z$ in the unit ball.
167 } 68 fn value_diff_unit_lipschitz_factor(&self) -> Option<Self::FloatType> {
168 69 None
169 pub fn n_sensors(&self) -> usize {
170 self.sensor_count.iter().product()
171 }
172
173 #[inline]
174 fn shifted_sensor(&self, x : Loc<F, N>) -> ShiftedSensor<F, S, P, N> {
175 self.base_sensor.clone().shift(x)
176 }
177
178 #[inline]
179 fn _zero_observable(&self) -> DVector<F> {
180 DVector::zeros(self.n_sensors())
181 } 70 }
182 } 71 }
183 72
184 impl<F, S, P, BT, const N : usize> Apply<RNDM<F, N>> for SensorGrid<F, S, P, BT, N>
185 where F : Float,
186 BT : SensorGridBT<F, S, P, N>,
187 S : Sensor<F, N>,
188 P : Spread<F, N>,
189 Convolution<S, P> : Spread<F, N>,
190 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
191
192 type Output = DVector<F>;
193
194 #[inline]
195 fn apply(&self, μ : RNDM<F, N>) -> DVector<F> {
196 self.apply(&μ)
197 }
198 }
199
200 impl<'a, F, S, P, BT, const N : usize> Apply<&'a RNDM<F, N>> for SensorGrid<F, S, P, BT, N>
201 where F : Float,
202 BT : SensorGridBT<F, S, P, N>,
203 S : Sensor<F, N>,
204 P : Spread<F, N>,
205 Convolution<S, P> : Spread<F, N>,
206 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
207
208 type Output = DVector<F>;
209
210 fn apply(&self, μ : &'a RNDM<F, N>) -> DVector<F> {
211 let mut res = self._zero_observable();
212 self.apply_add(&mut res, μ);
213 res
214 }
215 }
216
217 impl<F, S, P, BT, const N : usize> Linear<RNDM<F, N>> for SensorGrid<F, S, P, BT, N>
218 where F : Float,
219 BT : SensorGridBT<F, S, P, N>,
220 S : Sensor<F, N>,
221 P : Spread<F, N>,
222 Convolution<S, P> : Spread<F, N>,
223 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
224 type Codomain = DVector<F>;
225 }
226
227
228 #[replace_float_literals(F::cast_from(literal))]
229 impl<F, S, P, BT, const N : usize> GEMV<F, RNDM<F, N>, DVector<F>> for SensorGrid<F, S, P, BT, N>
230 where F : Float,
231 BT : SensorGridBT<F, S, P, N>,
232 S : Sensor<F, N>,
233 P : Spread<F, N>,
234 Convolution<S, P> : Spread<F, N>,
235 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
236
237 fn gemv(&self, y : &mut DVector<F>, α : F, μ : &RNDM<F, N>, β : F) {
238 let grid = self.grid();
239 if β == 0.0 {
240 y.fill(0.0)
241 } else if β != 1.0 {
242 *y *= β; // Need to multiply first, as we have to be able to add to y.
243 }
244 if α == 1.0 {
245 self.apply_add(y, μ)
246 } else {
247 for δ in μ.iter_spikes() {
248 for &d in self.bt.iter_at(&δ.x) {
249 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
250 y[d] += sensor.apply(&δ.x) * (α * δ.α);
251 }
252 }
253 }
254 }
255
256 fn apply_add(&self, y : &mut DVector<F>, μ : &RNDM<F, N>) {
257 let grid = self.grid();
258 for δ in μ.iter_spikes() {
259 for &d in self.bt.iter_at(&δ.x) {
260 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
261 y[d] += sensor.apply(&δ.x) * δ.α;
262 }
263 }
264 }
265
266 }
267
268 impl<F, S, P, BT, const N : usize> Apply<DeltaMeasure<Loc<F, N>, F>>
269 for SensorGrid<F, S, P, BT, N>
270 where F : Float,
271 BT : SensorGridBT<F, S, P, N>,
272 S : Sensor<F, N>,
273 P : Spread<F, N>,
274 Convolution<S, P> : Spread<F, N>,
275 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
276
277 type Output = DVector<F>;
278
279 #[inline]
280 fn apply(&self, δ : DeltaMeasure<Loc<F, N>, F>) -> DVector<F> {
281 self.apply(&δ)
282 }
283 }
284
285 impl<'a, F, S, P, BT, const N : usize> Apply<&'a DeltaMeasure<Loc<F, N>, F>>
286 for SensorGrid<F, S, P, BT, N>
287 where F : Float,
288 BT : SensorGridBT<F, S, P, N>,
289 S : Sensor<F, N>,
290 P : Spread<F, N>,
291 Convolution<S, P> : Spread<F, N>,
292 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
293
294 type Output = DVector<F>;
295
296 fn apply(&self, δ : &DeltaMeasure<Loc<F, N>, F>) -> DVector<F> {
297 let mut res = DVector::zeros(self.n_sensors());
298 let grid = self.grid();
299 for &d in self.bt.iter_at(&δ.x) {
300 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
301 res[d] += sensor.apply(&δ.x) * δ.α;
302 }
303 res
304 }
305 }
306
307 impl<F, S, P, BT, const N : usize> Linear<DeltaMeasure<Loc<F, N>, F>> for SensorGrid<F, S, P, BT, N>
308 where F : Float,
309 BT : SensorGridBT<F, S, P, N>,
310 S : Sensor<F, N>,
311 P : Spread<F, N>,
312 Convolution<S, P> : Spread<F, N>,
313 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
314 type Codomain = DVector<F>;
315 }
316
317 impl<F, S, P, BT, const N : usize> BoundedLinear<RNDM<F, N>> for SensorGrid<F, S, P, BT, N>
318 where F : Float,
319 BT : SensorGridBT<F, S, P, N, Agg=Bounds<F>>,
320 S : Sensor<F, N>,
321 P : Spread<F, N>,
322 Convolution<S, P> : Spread<F, N>,
323 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> {
324 type FloatType = F;
325
326 /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped
327 /// with the Radon norm, and $ℝ^n$ with the Euclidean norm.
328 fn opnorm_bound(&self) -> F {
329 // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have
330 // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩
331 // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ
332 // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ
333 // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ ∑ |z_i| |μ|_ℳ
334 // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √n |z|_2 |μ|_ℳ
335 // = |φ|_∞ √n |μ|_ℳ.
336 // Hence
337 let n = F::cast_from(self.n_sensors());
338 self.base_sensor.bounds().uniform() * n.sqrt()
339 }
340 }
341
342 type SensorGridPreadjoint<'a, A, F, const N : usize> = PreadjointHelper<'a, A, RNDM<F,N>>;
343
344
345 impl<F, S, P, BT, const N : usize>
346 Preadjointable<RNDM<F, N>, DVector<F>>
347 for SensorGrid<F, S, P, BT, N>
348 where F : Float,
349 BT : SensorGridBT<F, S, P, N>,
350 S : Sensor<F, N>,
351 P : Spread<F, N>,
352 Convolution<S, P> : Spread<F, N>,
353 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
354 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
355 type PreadjointCodomain = BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>;
356 type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where Self : 'a;
357
358 fn preadjoint(&self) -> Self::Preadjoint<'_> {
359 PreadjointHelper::new(self)
360 }
361 }
362
363 #[derive(Clone,Debug)]
364 pub struct SensorGridSupportGenerator<F, S, P, const N : usize>
365 where F : Float,
366 S : Sensor<F, N>,
367 P : Spread<F, N> {
368 base_sensor : Convolution<S, P>,
369 grid : LinGrid<F, N>,
370 weights : DVector<F>
371 }
372
373 impl<F, S, P, const N : usize> SensorGridSupportGenerator<F, S, P, N>
374 where F : Float,
375 S : Sensor<F, N>,
376 P : Spread<F, N>,
377 Convolution<S, P> : Spread<F, N> {
378
379 #[inline]
380 fn construct_sensor(&self, id : usize, w : F) -> Weighted<ShiftedSensor<F, S, P, N>, F> {
381 let x = self.grid.entry_linear_unchecked(id);
382 self.base_sensor.clone().shift(x).weigh(w)
383 }
384
385 #[inline]
386 fn construct_sensor_and_id<'a>(&'a self, (id, w) : (usize, &'a F))
387 -> (usize, Weighted<ShiftedSensor<F, S, P, N>, F>) {
388 (id.into(), self.construct_sensor(id, *w))
389 }
390 }
391
392 impl<F, S, P, const N : usize> SupportGenerator<F, N>
393 for SensorGridSupportGenerator<F, S, P, N>
394 where F : Float,
395 S : Sensor<F, N>,
396 P : Spread<F, N>,
397 Convolution<S, P> : Spread<F, N> {
398 type Id = usize;
399 type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>;
400 type AllDataIter<'a> = MapX<'a, Zip<RangeFrom<usize>,
401 std::slice::Iter<'a, F>>,
402 Self,
403 (Self::Id, Self::SupportType)>
404 where Self : 'a;
405
406 #[inline]
407 fn support_for(&self, d : Self::Id) -> Self::SupportType {
408 self.construct_sensor(d, self.weights[d])
409 }
410
411 #[inline]
412 fn support_count(&self) -> usize {
413 self.weights.len()
414 }
415
416 #[inline]
417 fn all_data(&self) -> Self::AllDataIter<'_> {
418 (0..).zip(self.weights.as_slice().iter()).mapX(self, Self::construct_sensor_and_id)
419 }
420 }
421
422 /// Helper structure for constructing preadjoints of `S` where `S : Linear<X>`.
423 /// [`Linear`] needs to be implemented for each instance, but [`Adjointable`]
424 /// and [`BoundedLinear`] have blanket implementations.
425 #[derive(Clone,Debug)]
426 pub struct PreadjointHelper<'a, S : 'a, X> {
427 forward_op : &'a S,
428 _domain : PhantomData<X>
429 }
430
431 impl<'a, S : 'a, X> PreadjointHelper<'a, S, X> {
432 pub fn new(forward_op : &'a S) -> Self {
433 PreadjointHelper { forward_op, _domain: PhantomData }
434 }
435 }
436
437 impl<'a, X, Ypre, S> Adjointable<Ypre, X>
438 for PreadjointHelper<'a, S, X>
439 where Self : Linear<Ypre>,
440 S : Clone + Linear<X> {
441 type AdjointCodomain = S::Codomain;
442 type Adjoint<'b> = S where Self : 'b;
443 fn adjoint(&self) -> Self::Adjoint<'_> {
444 self.forward_op.clone()
445 }
446 }
447
448 impl<'a, X, Ypre, S> BoundedLinear<Ypre>
449 for PreadjointHelper<'a, S, X>
450 where Self : Linear<Ypre>,
451 S : 'a + Clone + BoundedLinear<X> {
452 type FloatType = S::FloatType;
453 fn opnorm_bound(&self) -> Self::FloatType {
454 self.forward_op.opnorm_bound()
455 }
456 }
457
458
459 impl<'a, 'b, F, S, P, BT, const N : usize> Apply<&'b DVector<F>>
460 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<F,N>>
461 where F : Float,
462 BT : SensorGridBT<F, S, P, N>,
463 S : Sensor<F, N>,
464 P : Spread<F, N>,
465 Convolution<S, P> : Spread<F, N>,
466 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
467 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
468
469 type Output = SensorGridBTFN<F, S, P, BT, N>;
470
471 fn apply(&self, x : &'b DVector<F>) -> Self::Output {
472 self.apply(x.clone())
473 }
474 }
475
476 impl<'a, F, S, P, BT, const N : usize> Apply<DVector<F>>
477 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<F,N>>
478 where F : Float,
479 BT : SensorGridBT<F, S, P, N>,
480 S : Sensor<F, N>,
481 P : Spread<F, N>,
482 Convolution<S, P> : Spread<F, N>,
483 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
484 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
485
486 type Output = SensorGridBTFN<F, S, P, BT, N>;
487
488 fn apply(&self, x : DVector<F>) -> Self::Output {
489 let fwd = &self.forward_op;
490 let generator = SensorGridSupportGenerator{
491 base_sensor : fwd.base_sensor.clone(),
492 grid : fwd.grid(),
493 weights : x
494 };
495 BTFN::new_refresh(&fwd.bt, generator)
496 }
497 }
498
499 impl<'a, F, S, P, BT, const N : usize> Linear<DVector<F>>
500 for PreadjointHelper<'a, SensorGrid<F, S, P, BT, N>, RNDM<F,N>>
501 where F : Float,
502 BT : SensorGridBT<F, S, P, N>,
503 S : Sensor<F, N>,
504 P : Spread<F, N>,
505 Convolution<S, P> : Spread<F, N>,
506 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
507 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
508
509 type Codomain = SensorGridBTFN<F, S, P, BT, N>;
510 }
511
512 impl<F, S, P, BT, const N : usize> ForwardModel<Loc<F, N>, F>
513 for SensorGrid<F, S, P, BT, N>
514 where F : Float + ToNalgebraRealField<MixedType=F> + nalgebra::RealField,
515 BT : SensorGridBT<F, S, P, N>,
516 S : Sensor<F, N>,
517 P : Spread<F, N>,
518 Convolution<S, P> : Spread<F, N>,
519 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
520 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
521 type Observable = DVector<F>;
522
523 fn findim_quadratic_model(
524 &self,
525 μ : &DiscreteMeasure<Loc<F, N>, F>,
526 b : &Self::Observable
527 ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>) {
528 assert_eq!(b.len(), self.n_sensors());
529 let mut mA = DMatrix::zeros(self.n_sensors(), μ.len());
530 let grid = self.grid();
531 for (mut mAcol, δ) in mA.column_iter_mut().zip(μ.iter_spikes()) {
532 for &d in self.bt.iter_at(&δ.x) {
533 let sensor = self.shifted_sensor(grid.entry_linear_unchecked(d));
534 mAcol[d] += sensor.apply(&δ.x);
535 }
536 }
537 let mAt = mA.transpose();
538 (&mAt * mA, &mAt * b)
539 }
540
541 fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError {
542 let it = self.grid().into_iter().zip(b.iter()).map(|(x, &v)| (x, v));
543 write_csv(it, prefix + ".txt")
544 }
545
546 #[inline]
547 fn zero_observable(&self) -> Self::Observable {
548 self._zero_observable()
549 }
550
551 #[inline]
552 fn empty_observable(&self) -> Self::Observable {
553 DVector::zeros(0)
554 }
555
556 }
557
558 /// Implements the calculation a factor $L$ such that $A_*A ≤ L 𝒟$ for $A$ the forward model
559 /// and $𝒟$ a seminorm of suitable form.
560 ///
561 /// **This assumes (but does not check) that the sensors are not overlapping.**
562 #[replace_float_literals(F::cast_from(literal))]
563 impl<'a, F, BT, S, P, K, const N : usize> Lipschitz<&'a ConvolutionOp<F, K, BT, N>>
564 for SensorGrid<F, S, P, BT, N>
565 where F : Float + nalgebra::RealField + ToNalgebraRealField,
566 BT : SensorGridBT<F, S, P, N>,
567 S : Sensor<F, N>,
568 P : Spread<F, N>,
569 Convolution<S, P> : Spread<F, N>,
570 K : SimpleConvolutionKernel<F, N>,
571 AutoConvolution<P> : BoundedBy<F, K> {
572
573 type FloatType = F;
574
575 fn lipschitz_factor(&self, seminorm : &'a ConvolutionOp<F, K, BT, N>) -> Option<F> {
576 // Sensors should not take on negative values to allow
577 // A_*A to be upper bounded by a simple convolution of `spread`.
578 if self.sensor.bounds().lower() < 0.0 {
579 return None
580 }
581
582 // Calculate the factor $L_1$ for betwee $ℱ[ψ * ψ] ≤ L_1 ℱ[ρ]$ for $ψ$ the base spread
583 // and $ρ$ the kernel of the seminorm.
584 let l1 = AutoConvolution(self.spread.clone()).bounding_factor(seminorm.kernel())?;
585
586 // Calculate the factor for transitioning from $A_*A$ to `AutoConvolution<P>`, where A
587 // consists of several `Convolution<S, P>` for the physical model `P` and the sensor `S`.
588 let l0 = self.sensor.norm(Linfinity) * self.sensor.norm(L1);
589
590 // The final transition factor is:
591 Some(l0 * l1)
592 }
593 }
594
595 #[replace_float_literals(F::cast_from(literal))]
596 impl<F, BT, S, P, const N : usize> TransportLipschitz<L2Squared>
597 for SensorGrid<F, S, P, BT, N>
598 where F : Float + ToNalgebraRealField,
599 BT : SensorGridBT<F, S, P, N>,
600 S : Sensor<F, N>,
601 P : Spread<F, N>,
602 Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F> {
603 type FloatType = F;
604
605 fn transport_lipschitz_factor(&self, L2Squared : L2Squared) -> Self::FloatType {
606 // We estimate the factor by N_ψL^2, where L is the 2-norm Lipschitz factor of
607 // the base sensor (sensor * base_spread), and N_ψ the maximum overlap.
608 let l = self.base_sensor.lipschitz_factor(L2).unwrap();
609 let w = self.base_sensor.support_hint().width();
610 let d = map2(self.domain.width(), &self.sensor_count, |wi, &i| wi/F::cast_from(i));
611 let n = w.iter()
612 .zip(d.iter())
613 .map(|(&wi, &di)| (wi/di).ceil())
614 .reduce(F::mul)
615 .unwrap();
616 2.0 * n * l.powi(2)
617 }
618 }
619
620
621 macro_rules! make_sensorgridsupportgenerator_scalarop_rhs {
622 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
623 impl<F, S, P, const N : usize>
624 std::ops::$trait_assign<F>
625 for SensorGridSupportGenerator<F, S, P, N>
626 where F : Float,
627 S : Sensor<F, N>,
628 P : Spread<F, N>,
629 Convolution<S, P> : Spread<F, N> {
630 fn $fn_assign(&mut self, t : F) {
631 self.weights.$fn_assign(t);
632 }
633 }
634
635 impl<F, S, P, const N : usize>
636 std::ops::$trait<F>
637 for SensorGridSupportGenerator<F, S, P, N>
638 where F : Float,
639 S : Sensor<F, N>,
640 P : Spread<F, N>,
641 Convolution<S, P> : Spread<F, N> {
642 type Output = SensorGridSupportGenerator<F, S, P, N>;
643 fn $fn(mut self, t : F) -> Self::Output {
644 std::ops::$trait_assign::$fn_assign(&mut self.weights, t);
645 self
646 }
647 }
648
649 impl<'a, F, S, P, const N : usize>
650 std::ops::$trait<F>
651 for &'a SensorGridSupportGenerator<F, S, P, N>
652 where F : Float,
653 S : Sensor<F, N>,
654 P : Spread<F, N>,
655 Convolution<S, P> : Spread<F, N> {
656 type Output = SensorGridSupportGenerator<F, S, P, N>;
657 fn $fn(self, t : F) -> Self::Output {
658 SensorGridSupportGenerator{
659 base_sensor : self.base_sensor.clone(),
660 grid : self.grid,
661 weights : (&self.weights).$fn(t)
662 }
663 }
664 }
665 }
666 }
667
668 make_sensorgridsupportgenerator_scalarop_rhs!(Mul, mul, MulAssign, mul_assign);
669 make_sensorgridsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign);
670
671 macro_rules! make_sensorgridsupportgenerator_unaryop {
672 ($trait:ident, $fn:ident) => {
673 impl<F, S, P, const N : usize>
674 std::ops::$trait
675 for SensorGridSupportGenerator<F, S, P, N>
676 where F : Float,
677 S : Sensor<F, N>,
678 P : Spread<F, N>,
679 Convolution<S, P> : Spread<F, N> {
680 type Output = SensorGridSupportGenerator<F, S, P, N>;
681 fn $fn(mut self) -> Self::Output {
682 self.weights = self.weights.$fn();
683 self
684 }
685 }
686
687 impl<'a, F, S, P, const N : usize>
688 std::ops::$trait
689 for &'a SensorGridSupportGenerator<F, S, P, N>
690 where F : Float,
691 S : Sensor<F, N>,
692 P : Spread<F, N>,
693 Convolution<S, P> : Spread<F, N> {
694 type Output = SensorGridSupportGenerator<F, S, P, N>;
695 fn $fn(self) -> Self::Output {
696 SensorGridSupportGenerator{
697 base_sensor : self.base_sensor.clone(),
698 grid : self.grid,
699 weights : (&self.weights).$fn()
700 }
701 }
702 }
703 }
704 }
705
706 make_sensorgridsupportgenerator_unaryop!(Neg, neg);

Mercurial Repository: pointsource_algs

mercurial