Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model/sensor_grid.rs

src/forward_model/sensor_grid.rs

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

Mercurial Repository: pointsource_algs

mercurial