Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model.rs

src/forward_model.rs

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

Mercurial Repository: pointsource_algs

mercurial