Mercurial > repos > pointsource_algs / file comparison

comparison: src/forward_model/sensor_grid.rs

src/forward_model/sensor_grid.rs

branch
dev
changeset 39
6316d68b58af
parent 35
b087e3eab191
child 44
03251c546744
equal deleted inserted replaced
37:c5d8bd1a7728 39:6316d68b58af
92 pub struct SensorGrid<F, S, P, BT, const N : usize> 92 pub struct SensorGrid<F, S, P, BT, const N : usize>
93 where F : Float, 93 where F : Float,
94 S : Sensor<F, N>, 94 S : Sensor<F, N>,
95 P : Spread<F, N>, 95 P : Spread<F, N>,
96 Convolution<S, P> : Spread<F, N>, 96 Convolution<S, P> : Spread<F, N>,
97 BT : SensorGridBT<F, S, P, N>, { 97 BT : SensorGridBT<F, S, P, N>,
98 {
98 domain : Cube<F, N>, 99 domain : Cube<F, N>,
99 sensor_count : [usize; N], 100 sensor_count : [usize; N],
100 sensor : S, 101 sensor : S,
101 spread : P, 102 spread : P,
102 base_sensor : Convolution<S, P>, 103 base_sensor : Convolution<S, P>,
107 where F : Float, 108 where F : Float,
108 BT : SensorGridBT<F, S, P, N>, 109 BT : SensorGridBT<F, S, P, N>,
109 S : Sensor<F, N>, 110 S : Sensor<F, N>,
110 P : Spread<F, N>, 111 P : Spread<F, N>,
111 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 112 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>,
112 /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>*/ { 113 {
113 114
114 /// Create a new sensor grid. 115 /// Create a new sensor grid.
115 /// 116 ///
116 /// The parameter `depth` indicates the search depth of the created [`BT`]s 117 /// The parameter `depth` indicates the search depth of the created [`BT`]s
117 /// for the adjoint values. 118 /// for the adjoint values.
146 impl<F, S, P, BT, const N : usize> SensorGrid<F, S, P, BT, N> 147 impl<F, S, P, BT, const N : usize> SensorGrid<F, S, P, BT, N>
147 where F : Float, 148 where F : Float,
148 BT : SensorGridBT<F, S, P, N>, 149 BT : SensorGridBT<F, S, P, N>,
149 S : Sensor<F, N>, 150 S : Sensor<F, N>,
150 P : Spread<F, N>, 151 P : Spread<F, N>,
151 Convolution<S, P> : Spread<F, N> { 152 Convolution<S, P> : Spread<F, N>
153 {
152 154
153 /// Return the grid of sensor locations. 155 /// Return the grid of sensor locations.
154 pub fn grid(&self) -> LinGrid<F, N> { 156 pub fn grid(&self) -> LinGrid<F, N> {
155 lingrid_centered(&self.domain, &self.sensor_count) 157 lingrid_centered(&self.domain, &self.sensor_count)
156 } 158 }
188 F : Float, 190 F : Float,
189 BT : SensorGridBT<F, S, P, N>, 191 BT : SensorGridBT<F, S, P, N>,
190 S : Sensor<F, N>, 192 S : Sensor<F, N>,
191 P : Spread<F, N>, 193 P : Spread<F, N>,
192 Convolution<S, P> : Spread<F, N>, 194 Convolution<S, P> : Spread<F, N>,
193 //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
194 { 195 {
195 196
196 type Codomain = DVector<F>; 197 type Codomain = DVector<F>;
197 198
198 #[inline] 199 #[inline]
209 F : Float, 210 F : Float,
210 BT : SensorGridBT<F, S, P, N>, 211 BT : SensorGridBT<F, S, P, N>,
211 S : Sensor<F, N>, 212 S : Sensor<F, N>,
212 P : Spread<F, N>, 213 P : Spread<F, N>,
213 Convolution<S, P> : Spread<F, N>, 214 Convolution<S, P> : Spread<F, N>,
214 //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>
215 { } 215 { }
216 216
217 217
218 #[replace_float_literals(F::cast_from(literal))] 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> 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, 220 where F : Float,
221 BT : SensorGridBT<F, S, P, N>, 221 BT : SensorGridBT<F, S, P, N>,
222 S : Sensor<F, N>, 222 S : Sensor<F, N>,
223 P : Spread<F, N>, 223 P : Spread<F, N>,
224 Convolution<S, P> : Spread<F, N>, 224 Convolution<S, P> : Spread<F, N>,
225 //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>
226 { 225 {
227 226
228 fn gemv<I : Instance<RNDM<F, N>>>( 227 fn gemv<I : Instance<RNDM<F, N>>>(
229 &self, y : &mut DVector<F>, α : F, μ : I, β : F 228 &self, y : &mut DVector<F>, α : F, μ : I, β : F
230 ) { 229 ) {
266 for SensorGrid<F, S, P, BT, N> 265 for SensorGrid<F, S, P, BT, N>
267 where F : Float, 266 where F : Float,
268 BT : SensorGridBT<F, S, P, N, Agg=Bounds<F>>, 267 BT : SensorGridBT<F, S, P, N, Agg=Bounds<F>>,
269 S : Sensor<F, N>, 268 S : Sensor<F, N>,
270 P : Spread<F, N>, 269 P : Spread<F, N>,
271 Convolution<S, P> : Spread<F, N>, 270 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>
272 ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N> { 271 {
273 272
274 /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped 273 /// An estimate on the operator norm in $𝕃(ℳ(Ω); ℝ^n)$ with $ℳ(Ω)$ equipped
275 /// with the Radon norm, and $ℝ^n$ with the Euclidean norm. 274 /// with the Radon norm, and $ℝ^n$ with the Euclidean norm.
276 fn opnorm_bound(&self, _ : Radon, _ : L2) -> F { 275 fn opnorm_bound(&self, _ : Radon, _ : L2) -> F {
277 // With {x_i}_{i=1}^n the grid centres and φ the kernel, we have 276 // 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|μ⟩ 277 // |Aμ|_2 = sup_{|z|_2 ≤ 1} ⟨z,Αμ⟩ = sup_{|z|_2 ≤ 1} ⟨A^*z|μ⟩
279 // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ 278 // ≤ sup_{|z|_2 ≤ 1} |A^*z|_∞ |μ|_ℳ
280 // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ 279 // = sup_{|z|_2 ≤ 1} |∑ φ(· - x_i)z_i|_∞ |μ|_ℳ
281 // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ ∑ |z_i| |μ|_ℳ 280 // ≤ sup_{|z|_2 ≤ 1} |φ(y)| ∑_{i:th sensor active at y}|z_i| |μ|_ℳ
282 // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √n |z|_2 |μ|_ℳ 281 // where the supremum of |∑ φ(· - x_i)z_i|_∞ is reached at y
283 // = |φ|_∞ √n |μ|_ℳ. 282 // ≤ sup_{|z|_2 ≤ 1} |φ|_∞ √N_ψ |z|_2 |μ|_ℳ
283 // where N_ψ is the maximum number of sensors that overlap, and
284 // |z|_2 is restricted to the active sensors.
285 // = |φ|_∞ √N_ψ |μ|_ℳ.
284 // Hence 286 // Hence
285 let n = F::cast_from(self.n_sensors()); 287 let n = self.max_overlapping();
286 self.base_sensor.bounds().uniform() * n.sqrt() 288 self.base_sensor.bounds().uniform() * n.sqrt()
287 } 289 }
288 } 290 }
289 291
290 type SensorGridPreadjoint<'a, A, F, const N : usize> = PreadjointHelper<'a, A, RNDM<F,N>>; 292 type SensorGridPreadjoint<'a, A, F, const N : usize> = PreadjointHelper<'a, A, RNDM<F,N>>;
295 for SensorGrid<F, S, P, BT, N> 297 for SensorGrid<F, S, P, BT, N>
296 where F : Float, 298 where F : Float,
297 BT : SensorGridBT<F, S, P, N>, 299 BT : SensorGridBT<F, S, P, N>,
298 S : Sensor<F, N>, 300 S : Sensor<F, N>,
299 P : Spread<F, N>, 301 P : Spread<F, N>,
300 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 302 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>
301 /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, 303 {
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 PreadjointCodomain = BTFN<F, SensorGridSupportGenerator<F, S, P, N>, BT, N>;
304 type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where Self : 'a; 305 type Preadjoint<'a> = SensorGridPreadjoint<'a, Self, F, N> where Self : 'a;
305 306
306 fn preadjoint(&self) -> Self::Preadjoint<'_> { 307 fn preadjoint(&self) -> Self::Preadjoint<'_> {
307 PreadjointHelper::new(self) 308 PreadjointHelper::new(self)
315 BT : SensorGridBT<F, S, P, N>, 316 BT : SensorGridBT<F, S, P, N>,
316 S : Sensor<F, N>, 317 S : Sensor<F, N>,
317 P : Spread<F, N>, 318 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 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 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 {
321 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ {
322 322
323 type FloatType = F; 323 type FloatType = F;
324 324
325 fn value_unit_lipschitz_factor(&self) -> Option<F> { 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 326 // The Lipschitz factor of the sensors has to be scaled by the square root of twice
343 343
344 #[derive(Clone,Debug)] 344 #[derive(Clone,Debug)]
345 pub struct SensorGridSupportGenerator<F, S, P, const N : usize> 345 pub struct SensorGridSupportGenerator<F, S, P, const N : usize>
346 where F : Float, 346 where F : Float,
347 S : Sensor<F, N>, 347 S : Sensor<F, N>,
348 P : Spread<F, N> { 348 P : Spread<F, N>
349 {
349 base_sensor : Convolution<S, P>, 350 base_sensor : Convolution<S, P>,
350 grid : LinGrid<F, N>, 351 grid : LinGrid<F, N>,
351 weights : DVector<F> 352 weights : DVector<F>
352 } 353 }
353 354
354 impl<F, S, P, const N : usize> SensorGridSupportGenerator<F, S, P, N> 355 impl<F, S, P, const N : usize> SensorGridSupportGenerator<F, S, P, N>
355 where F : Float, 356 where F : Float,
356 S : Sensor<F, N>, 357 S : Sensor<F, N>,
357 P : Spread<F, N>, 358 P : Spread<F, N>,
358 Convolution<S, P> : Spread<F, N> { 359 Convolution<S, P> : Spread<F, N>
360 {
359 361
360 #[inline] 362 #[inline]
361 fn construct_sensor(&self, id : usize, w : F) -> Weighted<ShiftedSensor<F, S, P, N>, F> { 363 fn construct_sensor(&self, id : usize, w : F) -> Weighted<ShiftedSensor<F, S, P, N>, F> {
362 let x = self.grid.entry_linear_unchecked(id); 364 let x = self.grid.entry_linear_unchecked(id);
363 self.base_sensor.clone().shift(x).weigh(w) 365 self.base_sensor.clone().shift(x).weigh(w)
373 impl<F, S, P, const N : usize> SupportGenerator<F, N> 375 impl<F, S, P, const N : usize> SupportGenerator<F, N>
374 for SensorGridSupportGenerator<F, S, P, N> 376 for SensorGridSupportGenerator<F, S, P, N>
375 where F : Float, 377 where F : Float,
376 S : Sensor<F, N>, 378 S : Sensor<F, N>,
377 P : Spread<F, N>, 379 P : Spread<F, N>,
378 Convolution<S, P> : Spread<F, N> { 380 Convolution<S, P> : Spread<F, N>
381 {
379 type Id = usize; 382 type Id = usize;
380 type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>; 383 type SupportType = Weighted<ShiftedSensor<F, S, P, N>, F>;
381 type AllDataIter<'a> = MapX<'a, Zip<RangeFrom<usize>, 384 type AllDataIter<'a> = MapX<'a, Zip<RangeFrom<usize>,
382 std::slice::Iter<'a, F>>, 385 std::slice::Iter<'a, F>>,
383 Self, 386 Self,
405 where F : Float + ToNalgebraRealField<MixedType=F> + nalgebra::RealField, 408 where F : Float + ToNalgebraRealField<MixedType=F> + nalgebra::RealField,
406 BT : SensorGridBT<F, S, P, N>, 409 BT : SensorGridBT<F, S, P, N>,
407 S : Sensor<F, N>, 410 S : Sensor<F, N>,
408 P : Spread<F, N>, 411 P : Spread<F, N>,
409 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 412 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>,
410 /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, 413 {
411 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ {
412 type Observable = DVector<F>; 414 type Observable = DVector<F>;
413 415
414 fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError { 416 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)); 417 let it = self.grid().into_iter().zip(b.iter()).map(|(x, &v)| (x, v));
416 write_csv(it, prefix + ".txt") 418 write_csv(it, prefix + ".txt")
426 for SensorGrid<F, S, P, BT, N> 428 for SensorGrid<F, S, P, BT, N>
427 where F : Float + ToNalgebraRealField<MixedType=F> + nalgebra::RealField, 429 where F : Float + ToNalgebraRealField<MixedType=F> + nalgebra::RealField,
428 BT : SensorGridBT<F, S, P, N>, 430 BT : SensorGridBT<F, S, P, N>,
429 S : Sensor<F, N>, 431 S : Sensor<F, N>,
430 P : Spread<F, N>, 432 P : Spread<F, N>,
431 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>, 433 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, BT::Agg, N>
432 /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, 434 {
433 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ {
434 435
435 fn findim_quadratic_model( 436 fn findim_quadratic_model(
436 &self, 437 &self,
437 μ : &DiscreteMeasure<Loc<F, N>, F>, 438 μ : &DiscreteMeasure<Loc<F, N>, F>,
438 b : &Self::Observable 439 b : &Self::Observable
463 BT : SensorGridBT<F, S, P, N>, 464 BT : SensorGridBT<F, S, P, N>,
464 S : Sensor<F, N>, 465 S : Sensor<F, N>,
465 P : Spread<F, N>, 466 P : Spread<F, N>,
466 Convolution<S, P> : Spread<F, N>, 467 Convolution<S, P> : Spread<F, N>,
467 K : SimpleConvolutionKernel<F, N>, 468 K : SimpleConvolutionKernel<F, N>,
468 AutoConvolution<P> : BoundedBy<F, K> { 469 AutoConvolution<P> : BoundedBy<F, K>
470 {
469 471
470 type FloatType = F; 472 type FloatType = F;
471 473
472 fn adjoint_product_bound(&self, seminorm : &ConvolutionOp<F, K, BT, N>) -> Option<F> { 474 fn adjoint_product_bound(&self, seminorm : &ConvolutionOp<F, K, BT, N>) -> Option<F> {
473 // Sensors should not take on negative values to allow 475 // Sensors should not take on negative values to allow
494 for SensorGrid<F, S, P, BT, N> 496 for SensorGrid<F, S, P, BT, N>
495 where F : Float + ToNalgebraRealField, 497 where F : Float + ToNalgebraRealField,
496 BT : SensorGridBT<F, S, P, N>, 498 BT : SensorGridBT<F, S, P, N>,
497 S : Sensor<F, N>, 499 S : Sensor<F, N>,
498 P : Spread<F, N>, 500 P : Spread<F, N>,
499 Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F> { 501 Convolution<S, P> : Spread<F, N> + Lipschitz<L2, FloatType = F>
502 {
500 type FloatType = F; 503 type FloatType = F;
501 504
502 fn transport_lipschitz_factor(&self, L2Squared : L2Squared) -> Self::FloatType { 505 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 506 // 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. 507 // the base sensor (sensor * base_spread), and N_ψ the maximum overlap.
602 where F : Float, 605 where F : Float,
603 BT : SensorGridBT<F, S, P, N>, 606 BT : SensorGridBT<F, S, P, N>,
604 S : Sensor<F, N>, 607 S : Sensor<F, N>,
605 P : Spread<F, N>, 608 P : Spread<F, N>,
606 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, 609 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>,
607 //ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, 610 {
608 /*Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ {
609 611
610 type Codomain = SensorGridBTFN<F, S, P, BT, N>; 612 type Codomain = SensorGridBTFN<F, S, P, BT, N>;
611 613
612 fn apply<I : Instance<DVector<F>>>(&self, x : I) -> Self::Codomain { 614 fn apply<I : Instance<DVector<F>>>(&self, x : I) -> Self::Codomain {
613 let fwd = &self.forward_op; 615 let fwd = &self.forward_op;
625 where F : Float, 627 where F : Float,
626 BT : SensorGridBT<F, S, P, N>, 628 BT : SensorGridBT<F, S, P, N>,
627 S : Sensor<F, N>, 629 S : Sensor<F, N>,
628 P : Spread<F, N>, 630 P : Spread<F, N>,
629 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>, 631 Convolution<S, P> : Spread<F, N> + LocalAnalysis<F, Bounds<F>, N>,
630 /*ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>, 632 { }
631 Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N>*/ {
632
633 }
634

Mercurial Repository: pointsource_algs

mercurial