Mercurial > repos > pointsource_algs / file comparison

comparison: src/seminorms.rs

src/seminorms.rs

branch
dev
changeset 61
4f468d35fa29
parent 54
b3312eee105c
equal deleted inserted replaced
60:9738b51d90d7 61:4f468d35fa29
4 the trait [`DiscreteMeasureOp`]. 4 the trait [`DiscreteMeasureOp`].
5 */ 5 */
6 6
7 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, SpikeIter, RNDM}; 7 use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, SpikeIter, RNDM};
8 use alg_tools::bisection_tree::*; 8 use alg_tools::bisection_tree::*;
9 use alg_tools::bounds::Bounded;
10 use alg_tools::error::DynResult;
9 use alg_tools::instance::Instance; 11 use alg_tools::instance::Instance;
10 use alg_tools::iter::{FilterMapX, Mappable}; 12 use alg_tools::iter::{FilterMapX, Mappable};
11 use alg_tools::linops::{BoundedLinear, Linear, Mapping}; 13 use alg_tools::linops::{BoundedLinear, Linear, Mapping};
12 use alg_tools::loc::Loc; 14 use alg_tools::loc::Loc;
13 use alg_tools::mapping::RealMapping; 15 use alg_tools::mapping::RealMapping;
14 use alg_tools::nalgebra_support::ToNalgebraRealField; 16 use alg_tools::nalgebra_support::ToNalgebraRealField;
15 use alg_tools::norms::Linfinity; 17 use alg_tools::norms::{Linfinity, Norm, NormExponent};
16 use alg_tools::sets::Cube; 18 use alg_tools::sets::Cube;
17 use alg_tools::types::*; 19 use alg_tools::types::*;
18 use itertools::Itertools; 20 use itertools::Itertools;
19 use nalgebra::DMatrix; 21 use nalgebra::DMatrix;
20 use std::iter::Zip; 22 use std::iter::Zip;
66 // 68 //
67 // Convolutions for discrete measures 69 // Convolutions for discrete measures
68 // 70 //
69 71
70 /// A trait alias for simple convolution kernels. 72 /// A trait alias for simple convolution kernels.
71 pub trait SimpleConvolutionKernel<F: Float, const N: usize>: 73 pub trait SimpleConvolutionKernel<const N: usize, F: Float = f64>:
72 RealMapping<F, N> + Support<F, N> + Bounded<F> + Clone + 'static 74 RealMapping<N, F> + Support<N, F> + Bounded<F> + Clone + 'static
73 { 75 {
74 } 76 }
75 77
76 impl<T, F: Float, const N: usize> SimpleConvolutionKernel<F, N> for T where 78 impl<T, F: Float, const N: usize> SimpleConvolutionKernel<N, F> for T where
77 T: RealMapping<F, N> + Support<F, N> + Bounded<F> + Clone + 'static 79 T: RealMapping<N, F> + Support<N, F> + Bounded<F> + Clone + 'static
78 { 80 {
79 } 81 }
80 82
81 /// [`SupportGenerator`] for [`ConvolutionOp`]. 83 /// [`SupportGenerator`] for [`ConvolutionOp`].
82 #[derive(Clone, Debug)] 84 #[derive(Clone, Debug)]
83 pub struct ConvolutionSupportGenerator<F: Float, K, const N: usize> 85 pub struct ConvolutionSupportGenerator<K, const N: usize, F: Float = f64>
84 where 86 where
85 K: SimpleConvolutionKernel<F, N>, 87 K: SimpleConvolutionKernel<N, F>,
86 { 88 {
87 kernel: K, 89 kernel: K,
88 centres: RNDM<F, N>, 90 centres: RNDM<N, F>,
89 } 91 }
90 92
91 impl<F: Float, K, const N: usize> ConvolutionSupportGenerator<F, K, N> 93 impl<F: Float, K, const N: usize> ConvolutionSupportGenerator<K, N, F>
92 where 94 where
93 K: SimpleConvolutionKernel<F, N>, 95 K: SimpleConvolutionKernel<N, F>,
94 { 96 {
95 /// Construct the convolution kernel corresponding to `δ`, i.e., one centered at `δ.x` and 97 /// Construct the convolution kernel corresponding to `δ`, i.e., one centered at `δ.x` and
96 /// weighted by `δ.α`. 98 /// weighted by `δ.α`.
97 #[inline] 99 #[inline]
98 fn construct_kernel<'a>( 100 fn construct_kernel<'a>(
99 &'a self, 101 &'a self,
100 δ: &'a DeltaMeasure<Loc<F, N>, F>, 102 δ: &'a DeltaMeasure<Loc<N, F>, F>,
101 ) -> Weighted<Shift<K, F, N>, F> { 103 ) -> Weighted<Shift<K, N, F>, F> {
102 self.kernel.clone().shift(δ.x).weigh(δ.α) 104 self.kernel.clone().shift(δ.x).weigh(δ.α)
103 } 105 }
104 106
105 /// This is a helper method for the implementation of [`ConvolutionSupportGenerator::all_data`]. 107 /// This is a helper method for the implementation of [`ConvolutionSupportGenerator::all_data`].
106 /// It filters out `δ` with zero weight, and otherwise returns the corresponding convolution 108 /// It filters out `δ` with zero weight, and otherwise returns the corresponding convolution
107 /// kernel. The `id` is passed through as-is. 109 /// kernel. The `id` is passed through as-is.
108 #[inline] 110 #[inline]
109 fn construct_kernel_and_id_filtered<'a>( 111 fn construct_kernel_and_id_filtered<'a>(
110 &'a self, 112 &'a self,
111 (id, δ): (usize, &'a DeltaMeasure<Loc<F, N>, F>), 113 (id, δ): (usize, &'a DeltaMeasure<Loc<N, F>, F>),
112 ) -> Option<(usize, Weighted<Shift<K, F, N>, F>)> { 114 ) -> Option<(usize, Weighted<Shift<K, N, F>, F>)> {
113 (δ.α != F::ZERO).then(|| (id.into(), self.construct_kernel(δ))) 115 (δ.α != F::ZERO).then(|| (id.into(), self.construct_kernel(δ)))
114 } 116 }
115 } 117 }
116 118
117 impl<F: Float, K, const N: usize> SupportGenerator<F, N> for ConvolutionSupportGenerator<F, K, N> 119 impl<F: Float, K, const N: usize> SupportGenerator<N, F> for ConvolutionSupportGenerator<K, N, F>
118 where 120 where
119 K: SimpleConvolutionKernel<F, N>, 121 K: SimpleConvolutionKernel<N, F>,
120 { 122 {
121 type Id = usize; 123 type Id = usize;
122 type SupportType = Weighted<Shift<K, F, N>, F>; 124 type SupportType = Weighted<Shift<K, N, F>, F>;
123 type AllDataIter<'a> = FilterMapX< 125 type AllDataIter<'a> = FilterMapX<
124 'a, 126 'a,
125 Zip<RangeFrom<usize>, SpikeIter<'a, Loc<F, N>, F>>, 127 Zip<RangeFrom<usize>, SpikeIter<'a, Loc<N, F>, F>>,
126 Self, 128 Self,
127 (Self::Id, Self::SupportType), 129 (Self::Id, Self::SupportType),
128 >; 130 >;
129 131
130 #[inline] 132 #[inline]
148 /// Representation of a convolution operator $𝒟$. 150 /// Representation of a convolution operator $𝒟$.
149 #[derive(Clone, Debug)] 151 #[derive(Clone, Debug)]
150 pub struct ConvolutionOp<F, K, BT, const N: usize> 152 pub struct ConvolutionOp<F, K, BT, const N: usize>
151 where 153 where
152 F: Float + ToNalgebraRealField, 154 F: Float + ToNalgebraRealField,
153 BT: BTImpl<F, N, Data = usize>, 155 BT: BTImpl<N, F, Data = usize>,
154 K: SimpleConvolutionKernel<F, N>, 156 K: SimpleConvolutionKernel<N, F>,
155 { 157 {
156 /// Depth of the [`BT`] bisection tree for the outputs [`Mapping::apply`]. 158 /// Depth of the [`BT`] bisection tree for the outputs [`Mapping::apply`].
157 depth: BT::Depth, 159 depth: BT::Depth,
158 /// Domain of the [`BT`] bisection tree for the outputs [`Mapping::apply`]. 160 /// Domain of the [`BT`] bisection tree for the outputs [`Mapping::apply`].
159 domain: Cube<F, N>, 161 domain: Cube<N, F>,
160 /// The convolution kernel 162 /// The convolution kernel
161 kernel: K, 163 kernel: K,
162 _phantoms: PhantomData<(F, BT)>, 164 _phantoms: PhantomData<(F, BT)>,
163 } 165 }
164 166
165 impl<F, K, BT, const N: usize> ConvolutionOp<F, K, BT, N> 167 impl<F, K, BT, const N: usize> ConvolutionOp<F, K, BT, N>
166 where 168 where
167 F: Float + ToNalgebraRealField, 169 F: Float + ToNalgebraRealField,
168 BT: BTImpl<F, N, Data = usize>, 170 BT: BTImpl<N, F, Data = usize>,
169 K: SimpleConvolutionKernel<F, N>, 171 K: SimpleConvolutionKernel<N, F>,
170 { 172 {
171 /// Creates a new convolution operator $𝒟$ with `kernel` on `domain`. 173 /// Creates a new convolution operator $𝒟$ with `kernel` on `domain`.
172 /// 174 ///
173 /// The output of [`Mapping::apply`] is a [`BT`] of given `depth`. 175 /// The output of [`Mapping::apply`] is a [`BT`] of given `depth`.
174 pub fn new(depth: BT::Depth, domain: Cube<F, N>, kernel: K) -> Self { 176 pub fn new(depth: BT::Depth, domain: Cube<N, F>, kernel: K) -> Self {
175 ConvolutionOp { 177 ConvolutionOp {
176 depth: depth, 178 depth: depth,
177 domain: domain, 179 domain: domain,
178 kernel: kernel, 180 kernel: kernel,
179 _phantoms: PhantomData, 181 _phantoms: PhantomData,
180 } 182 }
181 } 183 }
182 184
183 /// Returns the support generator for this convolution operator. 185 /// Returns the support generator for this convolution operator.
184 fn support_generator(&self, μ: RNDM<F, N>) -> ConvolutionSupportGenerator<F, K, N> { 186 fn support_generator(&self, μ: RNDM<N, F>) -> ConvolutionSupportGenerator<K, N, F> {
185 // TODO: can we avoid cloning μ? 187 // TODO: can we avoid cloning μ?
186 ConvolutionSupportGenerator { 188 ConvolutionSupportGenerator {
187 kernel: self.kernel.clone(), 189 kernel: self.kernel.clone(),
188 centres: μ, 190 centres: μ,
189 } 191 }
193 pub fn kernel(&self) -> &K { 195 pub fn kernel(&self) -> &K {
194 &self.kernel 196 &self.kernel
195 } 197 }
196 } 198 }
197 199
198 impl<F, K, BT, const N: usize> Mapping<RNDM<F, N>> for ConvolutionOp<F, K, BT, N> 200 impl<F, K, BT, const N: usize> Mapping<RNDM<N, F>> for ConvolutionOp<F, K, BT, N>
199 where 201 where
200 F: Float + ToNalgebraRealField, 202 F: Float + ToNalgebraRealField,
201 BT: BTImpl<F, N, Data = usize>, 203 BT: BTImpl<N, F, Data = usize>,
202 K: SimpleConvolutionKernel<F, N>, 204 K: SimpleConvolutionKernel<N, F>,
203 Weighted<Shift<K, F, N>, F>: LocalAnalysis<F, BT::Agg, N>, 205 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
204 { 206 {
205 type Codomain = BTFN<F, ConvolutionSupportGenerator<F, K, N>, BT, N>; 207 type Codomain = BTFN<F, ConvolutionSupportGenerator<K, N, F>, BT, N>;
206 208
207 fn apply<I>(&self, μ: I) -> Self::Codomain 209 fn apply<I>(&self, μ: I) -> Self::Codomain
208 where 210 where
209 I: Instance<RNDM<F, N>>, 211 I: Instance<RNDM<N, F>>,
210 { 212 {
211 let g = self.support_generator(μ.own()); 213 let g = self.support_generator(μ.own());
212 BTFN::construct(self.domain.clone(), self.depth, g) 214 BTFN::construct(self.domain.clone(), self.depth, g)
213 } 215 }
214 } 216 }
215 217
216 /// [`ConvolutionOp`]s as linear operators over [`DiscreteMeasure`]s. 218 /// [`ConvolutionOp`]s as linear operators over [`DiscreteMeasure`]s.
217 impl<F, K, BT, const N: usize> Linear<RNDM<F, N>> for ConvolutionOp<F, K, BT, N> 219 impl<F, K, BT, const N: usize> Linear<RNDM<N, F>> for ConvolutionOp<F, K, BT, N>
218 where 220 where
219 F: Float + ToNalgebraRealField, 221 F: Float + ToNalgebraRealField,
220 BT: BTImpl<F, N, Data = usize>, 222 BT: BTImpl<N, F, Data = usize>,
221 K: SimpleConvolutionKernel<F, N>, 223 K: SimpleConvolutionKernel<N, F>,
222 Weighted<Shift<K, F, N>, F>: LocalAnalysis<F, BT::Agg, N>, 224 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
223 { 225 {
224 } 226 }
225 227
226 impl<F, K, BT, const N: usize> BoundedLinear<RNDM<F, N>, Radon, Linfinity, F> 228 impl<F, K, BT, const N: usize> BoundedLinear<RNDM<N, F>, Radon, Linfinity, F>
227 for ConvolutionOp<F, K, BT, N> 229 for ConvolutionOp<F, K, BT, N>
228 where 230 where
229 F: Float + ToNalgebraRealField, 231 F: Float + ToNalgebraRealField,
230 BT: BTImpl<F, N, Data = usize>, 232 BT: BTImpl<N, F, Data = usize>,
231 K: SimpleConvolutionKernel<F, N>, 233 K: SimpleConvolutionKernel<N, F>,
232 Weighted<Shift<K, F, N>, F>: LocalAnalysis<F, BT::Agg, N>, 234 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
233 { 235 {
234 fn opnorm_bound(&self, _: Radon, _: Linfinity) -> F { 236 fn opnorm_bound(&self, _: Radon, _: Linfinity) -> DynResult<F> {
235 // With μ = ∑_i α_i δ_{x_i}, we have 237 // With μ = ∑_i α_i δ_{x_i}, we have
236 // |𝒟μ|_∞ 238 // |𝒟μ|_∞
237 // = sup_z |∑_i α_i φ(z - x_i)| 239 // = sup_z |∑_i α_i φ(z - x_i)|
238 // ≤ sup_z ∑_i |α_i| |φ(z - x_i)| 240 // ≤ sup_z ∑_i |α_i| |φ(z - x_i)|
239 // ≤ ∑_i |α_i| |φ|_∞ 241 // ≤ ∑_i |α_i| |φ|_∞
240 // = |μ|_ℳ |φ|_∞ 242 // = |μ|_ℳ |φ|_∞
241 self.kernel.bounds().uniform() 243 Ok(self.kernel.bounds().uniform())
242 } 244 }
243 } 245 }
244 246
245 impl<F, K, BT, const N: usize> DiscreteMeasureOp<Loc<F, N>, F> for ConvolutionOp<F, K, BT, N> 247 impl<'a, F, K, BT, const N: usize> NormExponent for &'a ConvolutionOp<F, K, BT, N>
246 where 248 where
247 F: Float + ToNalgebraRealField, 249 F: Float + ToNalgebraRealField,
248 BT: BTImpl<F, N, Data = usize>, 250 BT: BTImpl<N, F, Data = usize>,
249 K: SimpleConvolutionKernel<F, N>, 251 K: SimpleConvolutionKernel<N, F>,
250 Weighted<Shift<K, F, N>, F>: LocalAnalysis<F, BT::Agg, N>, 252 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
251 { 253 {
252 type PreCodomain = PreBTFN<F, ConvolutionSupportGenerator<F, K, N>, N>; 254 }
255
256 impl<'a, F, K, BT, const N: usize> Norm<&'a ConvolutionOp<F, K, BT, N>, F> for RNDM<N, F>
257 where
258 F: Float + ToNalgebraRealField,
259 BT: BTImpl<N, F, Data = usize>,
260 K: SimpleConvolutionKernel<N, F>,
261 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
262 {
263 fn norm(&self, op𝒟: &'a ConvolutionOp<F, K, BT, N>) -> F {
264 self.apply(op𝒟.apply(self)).sqrt()
265 }
266 }
267
268 impl<F, K, BT, const N: usize> DiscreteMeasureOp<Loc<N, F>, F> for ConvolutionOp<F, K, BT, N>
269 where
270 F: Float + ToNalgebraRealField,
271 BT: BTImpl<N, F, Data = usize>,
272 K: SimpleConvolutionKernel<N, F>,
273 Weighted<Shift<K, N, F>, F>: LocalAnalysis<F, BT::Agg, N>,
274 {
275 type PreCodomain = PreBTFN<F, ConvolutionSupportGenerator<K, N, F>, N>;
253 276
254 fn findim_matrix<'a, I>(&self, points: I) -> DMatrix<F::MixedType> 277 fn findim_matrix<'a, I>(&self, points: I) -> DMatrix<F::MixedType>
255 where 278 where
256 I: ExactSizeIterator<Item = &'a Loc<F, N>> + Clone, 279 I: ExactSizeIterator<Item = &'a Loc<N, F>> + Clone,
257 { 280 {
258 // TODO: Preliminary implementation. It be best to use sparse matrices or 281 // TODO: Preliminary implementation. It be best to use sparse matrices or
259 // possibly explicit operators without matrices 282 // possibly explicit operators without matrices
260 let n = points.len(); 283 let n = points.len();
261 let points_clone = points.clone(); 284 let points_clone = points.clone();
266 } 289 }
267 290
268 /// A version of [`Mapping::apply`] that does not instantiate the [`BTFN`] codomain with 291 /// A version of [`Mapping::apply`] that does not instantiate the [`BTFN`] codomain with
269 /// a bisection tree, instead returning a [`PreBTFN`]. This can improve performance when 292 /// a bisection tree, instead returning a [`PreBTFN`]. This can improve performance when
270 /// the output is to be added as the right-hand-side operand to a proper BTFN. 293 /// the output is to be added as the right-hand-side operand to a proper BTFN.
271 fn preapply(&self, μ: RNDM<F, N>) -> Self::PreCodomain { 294 fn preapply(&self, μ: RNDM<N, F>) -> Self::PreCodomain {
272 BTFN::new_pre(self.support_generator(μ)) 295 BTFN::new_pre(self.support_generator(μ))
273 } 296 }
274 } 297 }
275 298
276 /// Generates an scalar operation (e.g. [`std::ops::Mul`], [`std::ops::Div`]) 299 /// Generates an scalar operation (e.g. [`std::ops::Mul`], [`std::ops::Div`])
277 /// for [`ConvolutionSupportGenerator`]. 300 /// for [`ConvolutionSupportGenerator`].
278 macro_rules! make_convolutionsupportgenerator_scalarop_rhs { 301 macro_rules! make_convolutionsupportgenerator_scalarop_rhs {
279 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { 302 ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => {
280 impl<F: Float, K: SimpleConvolutionKernel<F, N>, const N: usize> std::ops::$trait_assign<F> 303 impl<F: Float, K: SimpleConvolutionKernel<N, F>, const N: usize> std::ops::$trait_assign<F>
281 for ConvolutionSupportGenerator<F, K, N> 304 for ConvolutionSupportGenerator<K, N, F>
282 { 305 {
283 fn $fn_assign(&mut self, t: F) { 306 fn $fn_assign(&mut self, t: F) {
284 self.centres.$fn_assign(t); 307 self.centres.$fn_assign(t);
285 } 308 }
286 } 309 }
287 310
288 impl<F: Float, K: SimpleConvolutionKernel<F, N>, const N: usize> std::ops::$trait<F> 311 impl<F: Float, K: SimpleConvolutionKernel<N, F>, const N: usize> std::ops::$trait<F>
289 for ConvolutionSupportGenerator<F, K, N> 312 for ConvolutionSupportGenerator<K, N, F>
290 { 313 {
291 type Output = ConvolutionSupportGenerator<F, K, N>; 314 type Output = ConvolutionSupportGenerator<K, N, F>;
292 fn $fn(mut self, t: F) -> Self::Output { 315 fn $fn(mut self, t: F) -> Self::Output {
293 std::ops::$trait_assign::$fn_assign(&mut self.centres, t); 316 std::ops::$trait_assign::$fn_assign(&mut self.centres, t);
294 self 317 self
295 } 318 }
296 } 319 }
297 impl<'a, F: Float, K: SimpleConvolutionKernel<F, N>, const N: usize> std::ops::$trait<F> 320 impl<'a, F: Float, K: SimpleConvolutionKernel<N, F>, const N: usize> std::ops::$trait<F>
298 for &'a ConvolutionSupportGenerator<F, K, N> 321 for &'a ConvolutionSupportGenerator<K, N, F>
299 { 322 {
300 type Output = ConvolutionSupportGenerator<F, K, N>; 323 type Output = ConvolutionSupportGenerator<K, N, F>;
301 fn $fn(self, t: F) -> Self::Output { 324 fn $fn(self, t: F) -> Self::Output {
302 ConvolutionSupportGenerator { 325 ConvolutionSupportGenerator {
303 kernel: self.kernel.clone(), 326 kernel: self.kernel.clone(),
304 centres: (&self.centres).$fn(t), 327 centres: (&self.centres).$fn(t),
305 } 328 }
312 make_convolutionsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign); 335 make_convolutionsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign);
313 336
314 /// Generates an unary operation (e.g. [`std::ops::Neg`]) for [`ConvolutionSupportGenerator`]. 337 /// Generates an unary operation (e.g. [`std::ops::Neg`]) for [`ConvolutionSupportGenerator`].
315 macro_rules! make_convolutionsupportgenerator_unaryop { 338 macro_rules! make_convolutionsupportgenerator_unaryop {
316 ($trait:ident, $fn:ident) => { 339 ($trait:ident, $fn:ident) => {
317 impl<F: Float, K: SimpleConvolutionKernel<F, N>, const N: usize> std::ops::$trait 340 impl<F: Float, K: SimpleConvolutionKernel<N, F>, const N: usize> std::ops::$trait
318 for ConvolutionSupportGenerator<F, K, N> 341 for ConvolutionSupportGenerator<K, N, F>
319 { 342 {
320 type Output = ConvolutionSupportGenerator<F, K, N>; 343 type Output = ConvolutionSupportGenerator<K, N, F>;
321 fn $fn(mut self) -> Self::Output { 344 fn $fn(mut self) -> Self::Output {
322 self.centres = self.centres.$fn(); 345 self.centres = self.centres.$fn();
323 self 346 self
324 } 347 }
325 } 348 }
326 349
327 impl<'a, F: Float, K: SimpleConvolutionKernel<F, N>, const N: usize> std::ops::$trait 350 impl<'a, F: Float, K: SimpleConvolutionKernel<N, F>, const N: usize> std::ops::$trait
328 for &'a ConvolutionSupportGenerator<F, K, N> 351 for &'a ConvolutionSupportGenerator<K, N, F>
329 { 352 {
330 type Output = ConvolutionSupportGenerator<F, K, N>; 353 type Output = ConvolutionSupportGenerator<K, N, F>;
331 fn $fn(self) -> Self::Output { 354 fn $fn(self) -> Self::Output {
332 ConvolutionSupportGenerator { 355 ConvolutionSupportGenerator {
333 kernel: self.kernel.clone(), 356 kernel: self.kernel.clone(),
334 centres: (&self.centres).$fn(), 357 centres: (&self.centres).$fn(),
335 } 358 }

Mercurial Repository: pointsource_algs

mercurial