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