gstools_core/field.rs
1//! Compute the randomization methods for random field generations.
2
3use ndarray::{Array1, Array2, ArrayView1, ArrayView2, Axis, Zip};
4use rayon::prelude::*;
5
6use crate::short_vec::ShortVec;
7
8/// The randomization method for scalar fields.
9///
10/// Computes the isotropic spatial random field $u(x)$ by the randomization method according to
11/// $$
12/// u(x) = \sum_{i=1}^N (z_{1,i} \cdot \cos(\langle k_i, x \rangle) +
13/// z_{2,i} \cdot \sin(\langle k_i, x \rangle))
14/// $$
15/// with
16/// * $N$ being the number of Fourier modes
17/// * $z_1, z_2$ being independent samples from a standard normal distribution
18/// * $k$ being the samples from the spectral density distribution of the covariance model
19/// and are given by the argument `cov_samples`.
20///
21/// # Arguments
22///
23/// * `cov_samples` - the samples from the spectral density distribution of the covariance model
24/// <br> dim = (spatial dim. of field $d$, Fourier modes $N$)
25/// * `z1` - independent samples from a standard normal distribution
26/// <br> dim = Fourier modes $N$
27/// * `z2` - independent samples from a standard normal distribution
28/// <br> dim = Fourier modes $N$
29/// * `pos` - the position $x$ where the spatial random field is calculated
30/// <br> dim = (spatial dim. of field $d$, no. of spatial points where field is calculated $j$)
31/// * `num_threads` - the number of parallel threads used, if None, use rayon's default
32///
33/// # Returns
34///
35/// * `summed_modes` - the isotropic spatial field
36/// <br> dim = no. of spatial points where field is calculated $j$
37pub fn summator(
38 cov_samples: ArrayView2<'_, f64>,
39 z1: ArrayView1<'_, f64>,
40 z2: ArrayView1<'_, f64>,
41 pos: ArrayView2<'_, f64>,
42 num_threads: Option<usize>,
43) -> Array1<f64> {
44 assert_eq!(cov_samples.dim().0, pos.dim().0);
45 assert_eq!(cov_samples.dim().1, z1.dim());
46 assert_eq!(cov_samples.dim().1, z2.dim());
47
48 rayon::ThreadPoolBuilder::new()
49 .num_threads(num_threads.unwrap_or(rayon::current_num_threads()))
50 .build()
51 .unwrap()
52 .install(|| {
53 Zip::from(pos.columns()).par_map_collect(|pos| {
54 Zip::from(cov_samples.columns()).and(z1).and(z2).fold(
55 0.0,
56 |sum, sample, &z1, &z2| {
57 let phase = sample.dot(&pos);
58 let z12 = z1 * phase.cos() + z2 * phase.sin();
59
60 sum + z12
61 },
62 )
63 })
64 })
65}
66
67/// The randomization method for vector fields.
68///
69/// Computes the isotropic incompressible spatial random field $u(x)$ by the randomization method according to
70/// $$
71/// u(x)\_i = \sum_{j=1}^N p_i(k_j) (z_{1,j} \cdot \cos(\langle k_j, x \rangle) +
72/// z_{2,j} \cdot \sin(\langle k_j, x \rangle))
73/// $$
74/// with
75/// * $N$ being the number of Fourier modes
76/// * $z_1, z_2$ being independent samples from a standard normal distribution
77/// * $k$ being the samples from the spectral density distribution of the covariance model
78/// and are given by the argument `cov_samples`.
79/// * $p_i(k_j) = e_1 - \frac{k_ik_1}{|k|^2}$ being the projector ensuring the incompressibility
80///
81/// # Arguments
82///
83/// * `cov_samples` - the samples from the spectral density distribution of the covariance model
84/// <br> dim = (spatial dim. of field $d$, Fourier modes $N$)
85/// * `z1` - independent samples from a standard normal distribution
86/// <br> dim = Fourier modes $N$
87/// * `z2` - independent samples from a standard normal distribution
88/// <br> dim = Fourier modes $N$
89/// * `pos` - the position $x$ where the spatial random field is calculated
90/// <br> dim = (spatial dim. of field $d$, no. of spatial points where field is calculated $j$)
91/// * `num_threads` - the number of parallel threads used, if None, use rayon's default
92///
93/// # Returns
94///
95/// * `summed_modes` - the isotropic incompressible spatial field
96/// <br> dim = (spatial dim. of field $d$, no. of spatial points where field is calculated $j$)
97pub fn summator_incompr(
98 cov_samples: ArrayView2<'_, f64>,
99 z1: ArrayView1<'_, f64>,
100 z2: ArrayView1<'_, f64>,
101 pos: ArrayView2<'_, f64>,
102 num_threads: Option<usize>,
103) -> Array2<f64> {
104 assert_eq!(cov_samples.dim().0, pos.dim().0);
105 assert_eq!(cov_samples.dim().1, z1.dim());
106 assert_eq!(cov_samples.dim().1, z2.dim());
107
108 fn inner<const N: usize>(
109 cov_samples: ArrayView2<'_, f64>,
110 z1: ArrayView1<'_, f64>,
111 z2: ArrayView1<'_, f64>,
112 pos: ArrayView2<'_, f64>,
113 num_threads: Option<usize>,
114 ) -> Array2<f64> {
115 let cov_samples = cov_samples
116 .axis_iter(Axis(1))
117 .map(ShortVec::<N>::from_array)
118 .collect::<Vec<_>>();
119
120 let pos = pos
121 .axis_iter(Axis(1))
122 .map(ShortVec::<N>::from_array)
123 .collect::<Vec<_>>();
124
125 let summed_modes = rayon::ThreadPoolBuilder::new()
126 .num_threads(num_threads.unwrap_or(rayon::current_num_threads()))
127 .build()
128 .unwrap()
129 .install(|| {
130 cov_samples
131 .par_iter()
132 .zip(z1.axis_iter(Axis(0)))
133 .zip(z2.axis_iter(Axis(0)))
134 .with_min_len(100)
135 .fold(
136 || vec![ShortVec::<N>::zero(); pos.len()],
137 |mut summed_modes, ((cov_samples, z1), z2)| {
138 let k_2 = cov_samples[0] / cov_samples.dot(cov_samples);
139 let z1 = z1.into_scalar();
140 let z2 = z2.into_scalar();
141
142 pos.par_iter()
143 .zip(&mut summed_modes)
144 .for_each(|(pos, sum)| {
145 let phase = cov_samples.dot(pos);
146 let z12 = z1 * phase.cos() + z2 * phase.sin();
147
148 sum[0] += (1.0 - cov_samples[0] * k_2) * z12;
149
150 (1..N).for_each(|idx| {
151 sum[idx] -= cov_samples[idx] * k_2 * z12;
152 });
153 });
154
155 summed_modes
156 },
157 )
158 .reduce_with(|mut lhs, rhs| {
159 lhs.iter_mut().zip(&rhs).for_each(|(lhs, rhs)| lhs.add(rhs));
160
161 lhs
162 })
163 .unwrap()
164 });
165
166 Array2::<f64>::from_shape_vec(
167 (pos.len(), N),
168 summed_modes
169 .into_iter()
170 .flat_map(ShortVec::<N>::into_iter)
171 .collect(),
172 )
173 .unwrap()
174 .reversed_axes()
175 }
176
177 match pos.dim().0 {
178 2 => inner::<2>(cov_samples, z1, z2, pos, num_threads),
179 3 => inner::<3>(cov_samples, z1, z2, pos, num_threads),
180 _ => panic!("Only two- and three-dimensional problems are supported."),
181 }
182}
183
184/// The Fourier method for scalar fields.
185///
186/// Computes the periodic, isotropic spatial random field $u(x)$ by the Fourier
187/// method according to
188/// $$
189/// u(x) = \sum_{i=1}^N \sqrt{2S(k_i\Delta k}\left(
190/// z_{1,i} \cdot \cos(\langle k_i, x \rangle) +
191/// z_{2,i} \cdot \sin(\langle k_i, x \rangle)\right)
192/// $$
193/// with
194/// * $S$ being the spectrum of the covariance model
195/// * $N$ being the number of Fourier modes
196/// * $z_1, z_2$ being independent samples from a standard normal distribution
197/// * $k$ being the equidistant Fourier grid
198/// and are given by the argument `modes`.
199/// * $\Delta k$ being the cell size of the Fourier grid
200///
201/// # Arguments
202///
203/// * `spectrum_factor` - the pre-calculated factor $\sqrt{2S(k_i\Delta k)}
204/// <br> dim = Fourier modes $N$
205/// * `modes` - equidistant Fourier grid, $k$ in Eq. above
206/// <br> dim = (spatial dim. of field $d$, Fourier modes $N$)
207/// * `z1` - independent samples from a standard normal distribution
208/// <br> dim = Fourier modes $N$
209/// * `z2` - independent samples from a standard normal distribution
210/// <br> dim = Fourier modes $N$
211/// * `pos` - the position $x$ where the spatial random field is calculated
212/// <br> dim = (spatial dim. of field $d$, no. of spatial points where field is calculated $j$)
213/// * `num_threads` - the number of parallel threads used, if None, use rayon's default
214///
215/// # Returns
216///
217/// * `summed_modes` - the isotropic spatial field
218/// <br> dim = no. of spatial points where field is calculated $j$
219pub fn summator_fourier(
220 spectrum_factor: ArrayView1<'_, f64>,
221 modes: ArrayView2<'_, f64>,
222 z1: ArrayView1<'_, f64>,
223 z2: ArrayView1<'_, f64>,
224 pos: ArrayView2<'_, f64>,
225 num_threads: Option<usize>,
226) -> Array1<f64> {
227 assert_eq!(modes.dim().0, pos.dim().0);
228 assert_eq!(modes.dim().1, z1.dim());
229 assert_eq!(modes.dim().1, z2.dim());
230
231 rayon::ThreadPoolBuilder::new()
232 .num_threads(num_threads.unwrap_or(rayon::current_num_threads()))
233 .build()
234 .unwrap()
235 .install(|| {
236 Zip::from(pos.columns()).par_map_collect(|pos| {
237 Zip::from(modes.columns())
238 .and(spectrum_factor)
239 .and(z1)
240 .and(z2)
241 .fold(0.0, |sum, k, &spectrum_factor, &z1, &z2| {
242 let phase = k.dot(&pos);
243 let z12 = spectrum_factor * (z1 * phase.cos() + z2 * phase.sin());
244
245 sum + z12
246 })
247 })
248 })
249}
250
251#[cfg(test)]
252mod tests {
253 use super::*;
254
255 use approx::assert_ulps_eq;
256 use ndarray::{arr1, arr2};
257
258 struct Setup {
259 cov_samples: Array2<f64>,
260 z_1: Array1<f64>,
261 z_2: Array1<f64>,
262 pos: Array2<f64>,
263 }
264
265 impl Setup {
266 fn new() -> Self {
267 Self {
268 cov_samples: arr2(&[
269 [
270 -2.15, 1.04, 0.69, -1.09, -1.54, -2.32, -1.81, -2.78, 1.57, -3.44,
271 ],
272 [
273 0.19, -1.24, -2.10, -2.86, -0.63, -0.51, -1.68, -0.07, 0.29, -0.007,
274 ],
275 [
276 0.98, -2.83, -0.10, 3.23, 0.51, 0.13, -1.03, 1.53, -0.51, 2.82,
277 ],
278 ]),
279 z_1: arr1(&[
280 -1.93, 0.46, 0.66, 0.02, -0.10, 1.29, 0.93, -1.14, 1.81, 1.47,
281 ]),
282 z_2: arr1(&[
283 -0.26, 0.98, -1.30, 0.66, 0.57, -0.25, -0.31, -0.29, 0.69, 1.14,
284 ]),
285 pos: arr2(&[
286 [0.00, 1.43, 2.86, 4.29, 5.71, 7.14, 9.57, 10.00],
287 [-5.00, -3.57, -2.14, -0.71, 0.71, 2.14, 3.57, 5.00],
288 [-6.00, -4.00, -2.00, 0.00, 2.00, 4.00, 6.00, 8.00],
289 ]),
290 }
291 }
292 }
293
294 struct SetupFourier {
295 spectrum_factor: Array1<f64>,
296 modes: Array2<f64>,
297 z_1: Array1<f64>,
298 z_2: Array1<f64>,
299 pos: Array2<f64>,
300 }
301
302 impl SetupFourier {
303 fn new() -> Self {
304 Self {
305 spectrum_factor: arr1(&[
306 -2.15, 1.04, 0.69, -1.09, -1.54, -2.32, -1.81, -2.78, 1.57, -3.44,
307 ]),
308 modes: arr2(&[
309 [
310 -2.15, 1.04, 0.69, -1.09, -1.54, -2.32, -1.81, -2.78, 1.57, -3.44,
311 ],
312 [
313 0.19, -1.24, -2.10, -2.86, -0.63, -0.51, -1.68, -0.07, 0.29, -0.007,
314 ],
315 [
316 0.98, -2.83, -0.10, 3.23, 0.51, 0.13, -1.03, 1.53, -0.51, 2.82,
317 ],
318 ]),
319 z_1: arr1(&[
320 -1.93, 0.46, 0.66, 0.02, -0.10, 1.29, 0.93, -1.14, 1.81, 1.47,
321 ]),
322 z_2: arr1(&[
323 -0.26, 0.98, -1.30, 0.66, 0.57, -0.25, -0.31, -0.29, 0.69, 1.14,
324 ]),
325 pos: arr2(&[
326 [0.00, 1.43, 2.86, 4.29, 5.71, 7.14, 9.57, 10.00],
327 [-5.00, -3.57, -2.14, -0.71, 0.71, 2.14, 3.57, 5.00],
328 [-6.00, -4.00, -2.00, 0.00, 2.00, 4.00, 6.00, 8.00],
329 ]),
330 }
331 }
332 }
333
334 #[test]
335 fn test_summate_fourier_2d() {
336 let setup = SetupFourier::new();
337 assert_eq!(
338 summator_fourier(
339 setup.spectrum_factor.view(),
340 setup.modes.view(),
341 setup.z_1.view(),
342 setup.z_2.view(),
343 setup.pos.view(),
344 None,
345 ),
346 arr1(&[
347 1.0666558330143816,
348 -3.5855143411414883,
349 -2.70208228699285,
350 9.808554698975039,
351 0.01634921830347258,
352 -2.2356422006860663,
353 14.730786907708966,
354 -2.851408419726332,
355 ])
356 );
357 }
358
359 #[test]
360 fn test_summate_2d() {
361 let setup = Setup::new();
362
363 assert_eq!(
364 summator(
365 setup.cov_samples.view(),
366 setup.z_1.view(),
367 setup.z_2.view(),
368 setup.pos.view(),
369 None,
370 ),
371 arr1(&[
372 0.3773130601113641,
373 -4.298994445846448,
374 0.9285578931297425,
375 0.893013192171638,
376 -1.4956409956178418,
377 -1.488542499264307,
378 0.19211668257573278,
379 2.3427520079106143
380 ])
381 );
382 }
383
384 #[test]
385 fn test_summate_incompr_2d() {
386 let setup = Setup::new();
387
388 assert_ulps_eq!(
389 summator_incompr(
390 setup.cov_samples.view(),
391 setup.z_1.view(),
392 setup.z_2.view(),
393 setup.pos.view(),
394 None,
395 ),
396 arr2(&[
397 [
398 0.7026540940472319,
399 -1.9323916721330978,
400 -0.4166102970790725,
401 0.27803989953742114,
402 -2.0809691290114567,
403 0.20148641078244162,
404 0.7758364517737109,
405 0.12811415623445488
406 ],
407 [
408 0.3498241912898348,
409 -0.07775049450238455,
410 -0.5970579726508763,
411 0.03011066817308309,
412 -0.6406632397415202,
413 0.4669548537557405,
414 0.908893008714896,
415 -0.5120295866263118
416 ],
417 [
418 0.2838955719581232,
419 -0.9042103150526011,
420 -0.6494289973178196,
421 -0.5654019280252776,
422 -0.8386683161758316,
423 -0.4648269322196026,
424 -0.0656185245433833,
425 1.6593799470196355
426 ]
427 ]),
428 max_ulps = 6
429 );
430 }
431}