ground_motion_lib/
vectorized.rs

1//! # Vectorized Ground Motion Calculations
2//!
3//! This module provides parallelized routines for ground motion prediction calculations
4//! and post-processing over collections of site points. It leverages [`Rayon`] for efficient
5//! multi-threaded computation on large datasets.
6//!
7//! ## Features
8//!
9//! - **Parallel GMPE prediction:** Compute ground motion values for large sets of site points in parallel.
10//! - **Summary statistics computation:** Derive key statistical descriptors (mean, standard deviation, min, max, median)
11//!   from ground motion prediction results, with parallelism applied to key operations.
12//!
13//! ## Primary Types and Functions
14//!
15//! - [`calc_gmpe_vec`]: Perform parallel ground motion prediction for a vector of [`Vs30Point`] instances.
16//! - [`compute_stats`]: Calculate summary statistics over a collection of predicted [`GmpePoint`] values.
17//! - [`Stats`]: Struct representing the computed statistical summary.
18//!
19//! ## Parallelism
20//!
21//! This module uses [`Rayon`](https://docs.rs/rayon/latest/rayon/) for thread-safe, data-parallel operations:
22//!
23//! - `par_iter()` for distributing GMPE calculations and statistical reductions across threads.
24//! - Number of threads is controlled by the `RAYON_NUM_THREADS` environment variable or defaults
25//!   to the number of logical CPU cores.
26//!
27//! ## Usage Example
28//!
29//! ```rust
30//! use ground_motion_lib::gmm::{Vs30Point, Earthquake, Magnitude};
31//! use ground_motion_lib::configs::get_mf2013_lib_configs;
32//! use ground_motion_lib::vectorized::{calc_gmpe_vec, compute_stats};
33//!
34//! let points = vec![
35//!     Vs30Point::new(142.5, 50.0, 400., Some(200.), Some(0)),
36//!     Vs30Point::new(142.6, 50.1, 350., Some(150.), Some(1)),
37//! ];
38//!
39//! let eq = Earthquake {
40//!     lon: 142.4,
41//!     lat: 50.0,
42//!     depth: 10.0,
43//!     magnitude: 6.5,
44//!     magnitude_kind: Magnitude::Mw,
45//! };
46//!
47//! let gmpe_ref = get_mf2013_lib_configs().get("config_mf2013_crustal_pga").unwrap();
48//! let results = calc_gmpe_vec(&points, gmpe_ref, &eq);
49//! let stats = compute_stats(&results);
50//!
51//! println!("Stats: {stats:?}");
52//! ```
53//!
54//! ## See Also
55//!
56//! - [`crate::gmm::GmpePoint`]
57//! - [`crate::gmm::Vs30Point`]
58//! - [`crate::gmm::GroundMotionModeling`]
59//! - [`crate::mf2013::MF2013`]
60//!
61//! ## Thread Safety
62//!
63//! All operations in this module are thread-safe and make use of [`Rayon`] for concurrency.
64
65use crate::gmm::{Earthquake, GmpePoint, GroundMotionModeling, Vs30Point};
66use rayon::prelude::*;
67
68/// Calculate ground motion predictions for a set of site points in parallel.
69///
70/// This function takes a slice of `Vs30Point` site points, a reference to a ground motion prediction
71/// equation (GMPE) implementation, and an earthquake definition, and computes ground motion values
72/// (`GmpePoint`) for each site point using the provided GMPE model.
73///
74/// The calculation is performed in parallel using Rayon to improve performance for large datasets.
75///
76/// # Type Parameters
77///
78/// * `T` - A type implementing the `GroundMotionModeling` trait.
79///   Must also implement `Sync` to allow safe parallel access across threads.
80///
81/// # Arguments
82///
83/// * `points` - A slice of `Vs30Point` instances representing the site points for which
84///   ground motion predictions will be calculated.
85/// * `gmpe` - A reference to a type implementing the `GroundMotionModeling` trait, representing
86///   the GMPE model to be used for the calculations.
87/// * `eq` - A reference to the `Earthquake` instance describing the earthquake event.
88///
89/// # Returns
90///
91/// A `Vec<GmpePoint>` containing the calculated ground motion values for each input site point.
92///
93/// # Examples
94///
95/// ```rust
96/// use ground_motion_lib::gmm::{Vs30Point, Earthquake, Magnitude, GroundMotionModeling};
97/// use ground_motion_lib::mf2013::MF2013;
98/// use ground_motion_lib::configs::get_mf2013_lib_configs;
99/// use ground_motion_lib::vectorized::calc_gmpe_vec;
100///
101/// let points = vec![
102///     Vs30Point::new(142.5, 50.0, 400., Some(200.), Some(0)),
103///     Vs30Point::new(142.6, 50.1, 350., Some(150.), Some(1)),
104/// ];
105///
106/// let eq = Earthquake {
107///     lon: 142.4,
108///     lat: 50.0,
109///     depth: 10.0,
110///     magnitude: 6.5,
111///     magnitude_kind: Magnitude::Mw,
112/// };
113///
114/// let gmpe_ref = get_mf2013_lib_configs().get("config_mf2013_crustal_pga").unwrap();
115///
116/// let results = calc_gmpe_vec(&points, gmpe_ref, &eq);
117/// println!("{results:?}");
118/// ```
119///
120/// # Parallelism
121///
122/// This function uses Rayon’s `par_iter()` to distribute the workload of computing ground motion
123/// values across multiple threads. The number of threads is typically controlled by the
124/// `RAYON_NUM_THREADS` environment variable, or by default equals the number of logical CPU cores.
125/// - [`Rayon`](https://docs.rs/rayon/latest/rayon/)
126///
127/// # See Also
128///
129/// - [`GmpePoint`](crate::gmm::GmpePoint)
130/// - [`Vs30Point`](crate::gmm::Vs30Point)
131/// - [`GroundMotionModeling`](crate::gmm::GroundMotionModeling)
132///
133pub fn calc_gmpe_vec<T: GroundMotionModeling + Sync>(
134    points: &[Vs30Point],
135    gmpe: &T,
136    eq: &Earthquake,
137) -> Vec<GmpePoint> {
138    points
139        .par_iter()
140        .map(|point| point.get_gm(gmpe, eq))
141        .collect()
142}
143
144/// Struct for computed summary statistics
145#[derive(Debug, PartialEq)]
146pub struct Stats {
147    pub mean: f64,
148    pub std_dev: f64,
149    pub min: f64,
150    pub max: f64,
151    pub median: f64,
152}
153
154/// Compute summary statistics (mean, standard deviation, minimum, maximum, and median)
155/// for a list of `GmpePoint` values.
156///
157/// This function processes the `value` field from a slice of `GmpePoint` instances,
158/// computing key statistical measures.
159/// Most operations are parallelized using Rayon where possible (e.g., sum, variance, min, max)
160/// to improve performance on larger datasets.
161///
162/// # Arguments
163///
164/// * `points` - A slice of `GmpePoint` instances to analyze.
165///
166/// # Returns
167///
168/// A `Stats` struct containing:
169/// - `mean` — the arithmetic mean
170/// - `std_dev` — the sample standard deviation
171/// - `min` — the minimum value
172/// - `max` — the maximum value
173/// - `median` — the median value (sorted centrally)
174///
175/// # Example
176///
177/// ```rust
178/// use ground_motion_lib::gmm::GmpePoint;
179/// use ground_motion_lib::vectorized::{compute_stats, Stats};
180///
181/// let points = vec![
182///     GmpePoint::new_pga(147.1, 50.1, 1.1),
183///     GmpePoint::new_pga(147.2, 50.2, 1.2),
184///     GmpePoint::new_pga(147.3, 50.3, 1.3),
185/// ];
186///
187/// let stats = compute_stats(&points);
188///
189/// println!("Mean: {}", stats.mean);
190/// println!("Min: {}", stats.min);
191/// println!("Max: {}", stats.max);
192/// println!("Median: {}", stats.median);
193/// println!("Std Dev: {}", stats.std_dev);
194/// ```
195///
196/// # Parallelism
197///
198/// - Sum, variance, min, and max calculations use `Rayon`’s parallel iterators.
199/// - Median is computed single-threaded via an in-place sort since sorting isn’t parallelized here.
200///
201/// # Panics
202///
203/// This function will panic if called with an empty slice.
204///
205/// # See Also
206///
207/// - [`Rayon`](https://docs.rs/rayon/latest/rayon/)
208/// - [`GmpePoint`](crate::gmm::GmpePoint)
209/// - [`Stats`](crate::vectorized::Stats)
210///
211pub fn compute_stats(points: &[GmpePoint]) -> Stats {
212    let n = points.len() as f64;
213
214    // Extract values into a Vec<f64> to operate on
215    let mut values: Vec<f64> = points.iter().map(|p| p.value).collect();
216
217    // Compute sum in parallel, then mean
218    let sum: f64 = values.par_iter().sum();
219    let mean = sum / n;
220
221    // Compute variance (sample, denominator is n-1)
222    let variance: f64 = values
223        .par_iter()
224        .map(|v| {
225            let diff = v - mean;
226            diff * diff
227        })
228        .sum::<f64>()
229        / (n - 1.0);
230    let std_dev = variance.sqrt();
231
232    // Compute min and max via parallel reduction
233    let min = values
234        .par_iter()
235        .cloned()
236        .reduce(|| f64::INFINITY, f64::min);
237
238    let max = values
239        .par_iter()
240        .cloned()
241        .reduce(|| f64::NEG_INFINITY, f64::max);
242
243    // Compute median by sorting values locally (single-threaded)
244    values.sort_by(|a, b| a.partial_cmp(b).unwrap());
245    let median = if values.len() % 2 == 0 {
246        let mid = values.len() / 2;
247        (values[mid - 1] + values[mid]) / 2.0
248    } else {
249        values[values.len() / 2]
250    };
251
252    Stats {
253        mean,
254        std_dev,
255        min,
256        max,
257        median,
258    }
259}
260
261#[cfg(test)]
262mod tests {
263    use super::*;
264    use crate::gmm::GmpePointKind;
265
266    #[test]
267    fn test_compute_stats() {
268        let points = vec![
269            GmpePoint {
270                lon: 0.0,
271                lat: 0.0,
272                value: 1.0,
273                kind: GmpePointKind::Pga,
274            },
275            GmpePoint {
276                lon: 0.0,
277                lat: 0.0,
278                value: 2.0,
279                kind: GmpePointKind::Pga,
280            },
281            GmpePoint {
282                lon: 0.0,
283                lat: 0.0,
284                value: 3.0,
285                kind: GmpePointKind::Pga,
286            },
287            GmpePoint {
288                lon: 0.0,
289                lat: 0.0,
290                value: 4.0,
291                kind: GmpePointKind::Pga,
292            },
293            GmpePoint {
294                lon: 0.0,
295                lat: 0.0,
296                value: 5.0,
297                kind: GmpePointKind::Pga,
298            },
299        ];
300
301        let stats = compute_stats(&points);
302
303        // Expected values calculated by hand
304        let expected = Stats {
305            mean: 3.0,
306            std_dev: 1.5811388300841898, // sqrt(2.5)
307            min: 1.0,
308            max: 5.0,
309            median: 3.0,
310        };
311
312        assert!((stats.mean - expected.mean).abs() < 1e-10);
313        assert!((stats.std_dev - expected.std_dev).abs() < 1e-10);
314        assert_eq!(stats.min, expected.min);
315        assert_eq!(stats.max, expected.max);
316        assert_eq!(stats.median, expected.median);
317    }
318}
ground_motion_lib/vectorized.rs

ground_motion_lib/
vectorized.rs
