Mercurial > repos > pointsource_algs / file comparison

comparison: src/kernels/hat_convolution.rs

src/kernels/hat_convolution.rs

changeset 72
e9a460a0e638
parent 61
4f468d35fa29
equal deleted inserted replaced
71:e2953ffd4e0b 72:e9a460a0e638
95 type Codomain = S::Type; 95 type Codomain = S::Type;
96 96
97 #[inline] 97 #[inline]
98 fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain { 98 fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain {
99 let σ = self.radius(); 99 let σ = self.radius();
100 y.decompose().product_map(|x| self.value_1d_σ1(x / σ) / σ) 100 y.decompose().product_map(|x| self.value_1d_σ1(x / σ)) / σ.powi(N as i32)
101 } 101 }
102 } 102 }
103 103
104 #[replace_float_literals(S::Type::cast_from(literal))] 104 #[replace_float_literals(S::Type::cast_from(literal))]
105 impl<S, const N: usize> Lipschitz<L1> for HatConv<S, N> 105 impl<S, const N: usize> Lipschitz<L2> for HatConv<S, N>
106 where 106 where
107 S: Constant, 107 S: Constant,
108 { 108 {
109 type FloatType = S::Type; 109 type FloatType = S::Type;
110 #[inline] 110 #[inline]
111 fn lipschitz_factor(&self, L1: L1) -> DynResult<Self::FloatType> { 111 fn lipschitz_factor(&self, _L2: L2) -> DynResult<Self::FloatType> {
112 // For any ψ_i, we have 112 // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have
113 // ∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i) 113 // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j < i} ψ_j(y_j) ∏_{j > i} ψ_j(x_j)
114 // = [ψ_1(x_1)-ψ_1(y_1)] ∏_{i=2}^N ψ_i(x_i) 114 // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have
115 // + ψ_1(y_1)[ ∏_{i=2}^N ψ_i(x_i) - ∏_{i=2}^N ψ_i(y_i)] 115 // |ψ(x) - ψ(y)| ≤ ∑_i L_g|x_i-y_i| M_g^{n-1} ≤ √n L_g M_g^{n-1} |x-y|₂
116 // = ∑_{j=1}^N [ψ_j(x_j)-ψ_j(y_j)]∏_{i > j} ψ_i(x_i) ∏_{i < j} ψ_i(y_i) 116 // where L_g is the Lipschitz factor of g, and M_g a bound on it.
117 // Thus
118 // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)|
119 // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_j |ψ_j|
120 let σ = self.radius(); 117 let σ = self.radius();
121 let l1d = self.lipschitz_1d_σ1() / (σ * σ); 118 let l1d = self.lipschitz_1d_σ1();
122 let m1d = self.value_1d_σ1(0.0) / σ; 119 let m1d = self.value_1d_σ1(0.0);
123 Ok(l1d * m1d.powi(N as i32 - 1)) 120 Ok(
124 } 121 S::Type::cast_from(N).sqrt() * l1d * m1d.powi((N - 1) as i32)
125 } 122 / σ.powi((N - 1 + 2) as i32),
126 123 )
127 impl<S, const N: usize> Lipschitz<L2> for HatConv<S, N>
128 where
129 S: Constant,
130 {
131 type FloatType = S::Type;
132 #[inline]
133 fn lipschitz_factor(&self, L2: L2) -> DynResult<Self::FloatType> {
134 self.lipschitz_factor(L1)
135 .map(|l1| l1 * <S::Type>::cast_from(N).sqrt())
136 } 124 }
137 } 125 }
138 126
139 impl<'a, S, const N: usize> DifferentiableImpl<Loc<N, S::Type>> for HatConv<S, N> 127 impl<'a, S, const N: usize> DifferentiableImpl<Loc<N, S::Type>> for HatConv<S, N>
140 where 128 where
144 132
145 #[inline] 133 #[inline]
146 fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative { 134 fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative {
147 let y = y0.decompose(); 135 let y = y0.decompose();
148 let σ = self.radius(); 136 let σ = self.radius();
149 let σ2 = σ * σ; 137 let tmp = y.map(|x| x / σ);
150 let vs = y.map(|x| self.value_1d_σ1(x / σ) / σ); 138 let vs = tmp.map(|x| self.value_1d_σ1(x));
151 product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ) / σ2) 139 product_differential(&tmp, &vs, |x| self.diff_1d_σ1(x)) / σ.powi((N + 1) as i32)
152 } 140 }
153 } 141 }
154 142
155 #[replace_float_literals(S::Type::cast_from(literal))] 143 #[replace_float_literals(S::Type::cast_from(literal))]
156 impl<'a, F: Float, S, const N: usize> LipschitzDifferentiableImpl<Loc<N, S::Type>, L2> 144 impl<'a, F: Float, S, const N: usize> LipschitzDifferentiableImpl<Loc<N, S::Type>, L2>
162 150
163 #[inline] 151 #[inline]
164 fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> { 152 fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> {
165 let σ = self.radius(); 153 let σ = self.radius();
166 product_differential_lipschitz_factor::<F, N>( 154 product_differential_lipschitz_factor::<F, N>(
167 self.value_1d_σ1(0.0) / σ, 155 self.value_1d_σ1(0.0),
168 self.lipschitz_1d_σ1() / (σ * σ), 156 self.lipschitz_1d_σ1(),
169 self.maxabsdiff_1d_σ1() / (σ * σ), 157 self.maxabsdiff_1d_σ1(),
170 self.lipschitz_diff_1d_σ1() / (σ * σ), 158 self.lipschitz_diff_1d_σ1(),
171 ) 159 )
160 .map(|l| l / σ.powi((N + 2) as i32))
161 // We have f(x)=f₀(x/σ)/σ^N, where f₀(z)=∏_{i=1}^N f₁(z_i).
162 // Thus |f'(x)-f'(y)| = |f₀'(x/σ)-f₀'(y/σ)|/σ^{N+1} ≤ L_{f₀'}/σ^{N+2}.
172 } 163 }
173 } 164 }
174 165
175 #[replace_float_literals(S::Type::cast_from(literal))] 166 #[replace_float_literals(S::Type::cast_from(literal))]
176 impl<'a, F: Float, S, const N: usize> HatConv<S, N> 167 impl<'a, F: Float, S, const N: usize> HatConv<S, N>
193 } 184 }
194 185
195 /// Computes the differential of the kernel for $n=1$ with $σ=1$. 186 /// Computes the differential of the kernel for $n=1$ with $σ=1$.
196 #[inline] 187 #[inline]
197 fn diff_1d_σ1(&self, x: F) -> F { 188 fn diff_1d_σ1(&self, x: F) -> F {
198 let y = x.abs(); 189 if x >= 1.0 || x <= -1.0 {
199 if y >= 1.0 {
200 0.0 190 0.0
201 } else if y > 0.5 { 191 } else if x > 0.5 {
202 -8.0 * (y - 1.0).powi(2) 192 -8.0 * (x - 1.0).powi(2)
203 } else 193 } else if x < -0.5 {
204 /* 0 ≤ y ≤ 0.5 */ 194 8.0 * (x + 1.0).powi(2)
205 { 195 } else {
206 (24.0 * y - 16.0) * y 196 /* 0 ≤ y ≤ 0.5 */
197 (24.0 * x.abs() - 16.0) * x
207 } 198 }
208 } 199 }
209 200
210 /// Computes the Lipschitz factor of the kernel for $n=1$ with $σ=1$. 201 /// Computes the Lipschitz factor of the kernel for $n=1$ with $σ=1$.
211 #[inline] 202 #[inline]
217 /// Computes the maximum absolute differential of the kernel for $n=1$ with $σ=1$. 208 /// Computes the maximum absolute differential of the kernel for $n=1$ with $σ=1$.
218 #[inline] 209 #[inline]
219 fn maxabsdiff_1d_σ1(&self) -> F { 210 fn maxabsdiff_1d_σ1(&self) -> F {
220 // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis 211 // Maximal absolute differential achieved at ±0.5 by diff_1d_σ1 analysis
221 2.0 212 2.0
222 }
223
224 /// Computes the second differential of the kernel for $n=1$ with $σ=1$.
225 #[inline]
226 #[allow(dead_code)]
227 fn diff2_1d_σ1(&self, x: F) -> F {
228 let y = x.abs();
229 if y >= 1.0 {
230 0.0
231 } else if y > 0.5 {
232 -16.0 * (y - 1.0)
233 } else
234 /* 0 ≤ y ≤ 0.5 */
235 {
236 48.0 * y - 16.0
237 }
238 } 213 }
239 214
240 /// Computes the differential of the kernel for $n=1$ with $σ=1$. 215 /// Computes the differential of the kernel for $n=1$ with $σ=1$.
241 #[inline] 216 #[inline]
242 fn lipschitz_diff_1d_σ1(&self) -> F { 217 fn lipschitz_diff_1d_σ1(&self) -> F {
330 // $$ 305 // $$
331 // [χ_{-β,β} * u_σ](x) 306 // [χ_{-β,β} * u_σ](x)
332 // = ∫_{x-β}^{x+β} u_σ(z) d z 307 // = ∫_{x-β}^{x+β} u_σ(z) d z
333 // = (1/σ)∫_{x-β}^{x+β} u_1(z/σ) d z 308 // = (1/σ)∫_{x-β}^{x+β} u_1(z/σ) d z
334 // = ∫_{(x-β)/σ}^{(x+β)/σ} u_1(z) d z 309 // = ∫_{(x-β)/σ}^{(x+β)/σ} u_1(z) d z
335 // = [χ_{-β/σ, β/σ} * u_1](x/σ) 310 // = (χ_{[-β/σ, β/σ]} * u_1)(x/σ)
336 // $$ 311 // $$
337 self.value_1d_σ1(x / σ, β / σ) 312 self.value_1d_σ1(x / σ, β / σ)
338 }) 313 })
339 } 314 }
340 } 315 }
352 fn differential_impl<I: Instance<Loc<N, F>>>(&self, y0: I) -> Loc<N, F> { 327 fn differential_impl<I: Instance<Loc<N, F>>>(&self, y0: I) -> Loc<N, F> {
353 let y = y0.decompose(); 328 let y = y0.decompose();
354 let Convolution(ref ind, ref hatconv) = self; 329 let Convolution(ref ind, ref hatconv) = self;
355 let β = ind.r.value(); 330 let β = ind.r.value();
356 let σ = hatconv.radius(); 331 let σ = hatconv.radius();
357 let σ2 = σ * σ;
358 332
359 let vs = y.map(|x| self.value_1d_σ1(x / σ, β / σ)); 333 let vs = y.map(|x| self.value_1d_σ1(x / σ, β / σ));
360 product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ) / σ2) 334 product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ)) / σ
361 } 335 }
362 } 336 }
363 337
364 /// Integrate $f$, whose support is $[c, d]$, on $[a, b]$. 338 /// Integrate $f'$, whose support is $[c, d]$, on $[a, b]$.
339 /// The value $f$ is given, not the derivative.
365 /// If $b > d$, add $g()$ to the result. 340 /// If $b > d$, add $g()$ to the result.
366 #[inline] 341 #[inline]
367 #[replace_float_literals(F::cast_from(literal))] 342 #[replace_float_literals(F::cast_from(literal))]
368 fn i<F: Float>(a: F, b: F, c: F, d: F, f: impl Fn(F) -> F, g: impl Fn() -> F) -> F { 343 fn i<F: Float>(a: F, b: F, c: F, d: F, f: impl Fn(F) -> F, g: impl Fn() -> F) -> F {
369 if b < c { 344 if b < c {
401 pub fn value_1d_σ1(&self, x: F, β: F) -> F { 376 pub fn value_1d_σ1(&self, x: F, β: F) -> F {
402 // The integration interval 377 // The integration interval
403 let a = x - β; 378 let a = x - β;
404 let b = x + β; 379 let b = x + β;
405 380
406 #[inline]
407 fn pow4<F: Float>(x: F) -> F {
408 let y = x * x;
409 y * y
410 }
411
412 // Observe the factor 1/6 at the front from the antiderivatives below. 381 // Observe the factor 1/6 at the front from the antiderivatives below.
413 // The factor 4 is from normalisation of the original function. 382 // The factor 4 is from normalisation of the original function.
414 (4.0 / 6.0) 383 (4.0 / 6.0)
415 * i( 384 * i(
416 a, 385 a,
417 b, 386 b,
418 -1.0, 387 -1.0,
419 -0.5, 388 -0.5,
420 // (2/3) (y+1)^3 on -1 < y ≤ -1/2 389 // (2/3) (y+1)^3 on -1 < y ≤ -1/2
421 // The antiderivative is (2/12)(y+1)^4 = (1/6)(y+1)^4 390 // The antiderivative is (2/12)(y+1)^4 = (1/6)(y+1)^4
422 |y| pow4(y + 1.0), 391 |y| (y + 1.0).powi(4),
423 || { 392 || {
424 i( 393 i(
425 a, 394 a,
426 b, 395 b,
427 -0.5, 396 -0.5,
428 0.0, 397 0.5,
429 // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0 398 // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0
430 // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y 399 // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y
431 |y| y * (-y * y * (y * 3.0 + 4.0) + 2.0), 400 // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
401 // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y
402 |y| y * (y * y * (y.abs() * 3.0 - 4.0) + 2.0),
432 || { 403 || {
433 i( 404 i(
434 a, 405 a,
435 b, 406 b,
436 0.0,
437 0.5, 407 0.5,
438 // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 408 1.0,
439 // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y 409 // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
440 |y| y * (y * y * (y * 3.0 - 4.0) + 2.0), 410 // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4
441 || { 411 |y| -(y - 1.0).powi(4),
442 i( 412 || 0.0,
443 a,
444 b,
445 0.5,
446 1.0,
447 // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
448 // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4
449 |y| -pow4(y - 1.0),
450 || 0.0,
451 )
452 },
453 ) 413 )
454 }, 414 },
455 ) 415 )
456 }, 416 },
457 ) 417 )
464 pub fn diff_1d_σ1(&self, x: F, β: F) -> F { 424 pub fn diff_1d_σ1(&self, x: F, β: F) -> F {
465 // The integration interval 425 // The integration interval
466 let a = x - β; 426 let a = x - β;
467 let b = x + β; 427 let b = x + β;
468 428
469 // The factor 4 is from normalisation of the original function. 429 i(
470 4.0 * i(
471 a, 430 a,
472 b, 431 b,
473 -1.0, 432 -1.0,
474 -0.5, 433 -0.5,
475 // (2/3) (y+1)^3 on -1 < y ≤ -1/2 434 // (2/3) (y+1)^3 on -1 < y ≤ -1/2
476 |y| (2.0 / 3.0) * (y + 1.0).powi(3), 435 |y| (8.0 / 3.0) * (y + 1.0).powi(3),
477 || { 436 || {
478 i( 437 i(
479 a, 438 a,
480 b, 439 b,
481 -0.5, 440 -0.5,
482 0.0, 441 0.5,
483 // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0 442 // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2
484 |y| -2.0 * (y + 1.0) * y * y + (1.0 / 3.0), 443 |y| 8.0 * (y.abs() - 1.0) * y * y + (4.0 / 3.0),
485 || { 444 || {
486 i( 445 i(
487 a, 446 a,
488 b, 447 b,
489 0.0,
490 0.5, 448 0.5,
491 // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 449 1.0,
492 |y| 2.0 * (y - 1.0) * y * y + (1.0 / 3.0), 450 // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
493 || { 451 |y| -(8.0 / 3.0) * (y - 1.0).powi(3),
494 i( 452 || 0.0,
495 a,
496 b,
497 0.5,
498 1.0,
499 // -(2/3) (y-1)^3 on 1/2 < y ≤ 1
500 |y| -(2.0 / 3.0) * (y - 1.0).powi(3),
501 || 0.0,
502 )
503 },
504 ) 453 )
505 }, 454 },
506 ) 455 )
507 }, 456 },
508 ) 457 )
523 None 472 None
524 } 473 }
525 } 474 }
526 */ 475 */
527 476
477 #[replace_float_literals(F::cast_from(literal))]
528 impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>> 478 impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>>
529 where 479 where
530 R: Constant<Type = F>, 480 R: Constant<Type = F>,
531 C: Constant<Type = F>, 481 C: Constant<Type = F>,
532 { 482 {

Mercurial Repository: pointsource_algs

mercurial