Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model/sensor_grid.rs

src/forward_model/sensor_grid.rs

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

Mercurial Repository: pointsource_algs

mercurial