hpr_analysis/sensitivity/
benchmark.rs

1//! Test functions whose Sobol' indices are known in closed form, to check an analysis against.
2//!
3//! **Guide:** [Sensitivity analysis][guide]'s *Checked against* section shows both methods on
4//! them.
5//!
6//! [guide]: https://nrdptel.github.io/hpr-sim/sensitivity.html#checked-against
7//!
8//! - [`Ishigami`]: `y = sin x₁ + a sin² x₂ + b x₃⁴ sin x₁` on `[−π, π]³`. It is nonlinear and
9//!   not additive: `x₃` does nothing alone (`S₃ = 0`) but a quarter of the variance with `x₁`.
10//! - [`SobolG`]: `y = Πᵢ (|4xᵢ − 2| + aᵢ)/(1 + aᵢ)` on `[0, 1]^k`. Each `aᵢ ≥ 0` sets how much
11//!   factor `i` matters, from most at `aᵢ = 0` to almost nothing at `aᵢ = 99`.
12//!
13//! Each evaluates a point given as a slice, the form [`super::morris`] and [`super::sobol`] pass,
14//! and gives `NaN` for a point of the wrong length, which an analysis then refuses.
15
16use std::f64::consts::PI;
17
18use serde::{Deserialize, Serialize};
19
20use super::Factor;
21use crate::error::AnalysisError;
22
23/// Ishigami and Homma's function, `y = sin x₁ + a sin² x₂ + b x₃⁴ sin x₁`, each `xᵢ` uniform on
24/// `[−π, π]`.
25///
26/// Its mean is `a/2` and its variance and the parts of it are (I. M. Sobol' and Y. L. Levitan,
27/// "On the use of variance reducing multipliers in Monte Carlo computations of a global
28/// sensitivity index", *Computer Physics Communications* 117, 52–61, 1999,
29/// <https://doi.org/10.1016/S0010-4655(98)00156-8>, p. 57, §5; A. Saltelli and others, *Global
30/// Sensitivity Analysis: The Primer*, Wiley, 2008, pp. 179–182)
31///
32/// - `V = 1/2 + a²/8 + bπ⁴/5 + b²π⁸/18`,
33/// - `V₁ = (1 + bπ⁴/5)²/2`, `V₂ = a²/8`, `V₃ = 0`,
34/// - `V₁₃ = (1/18 − 1/50) b²π⁸ = 8b²π⁸/225`, every other part zero,
35///
36/// so `Sᵢ = Vᵢ/V`, `S_T₁ = (V₁ + V₁₃)/V`, `S_T₂ = V₂/V` and `S_T₃ = V₁₃/V`. The usual
37/// parameters, [`Ishigami::STANDARD`], are `a = 7`, `b = 0.1`. T. Ishigami and T. Homma, "An
38/// importance quantification technique in uncertainty analysis for computer models", *Proceedings
39/// of ISUMA '90*, 398–403, 1990, <https://doi.org/10.1109/ISUMA.1990.151285>, introduced it.
40#[derive(Debug, Clone, Copy, PartialEq, Serialize, Deserialize)]
41pub struct Ishigami {
42    /// The weight of `sin² x₂`.
43    pub a: f64,
44    /// The weight of `x₃⁴ sin x₁`.
45    pub b: f64,
46}
47
48impl Ishigami {
49    /// `a = 7`, `b = 0.1`.
50    pub const STANDARD: Self = Self { a: 7.0, b: 0.1 };
51
52    /// Its three factors, `x1` to `x3`, each on `[−π, π]`.
53    ///
54    /// # Errors
55    ///
56    /// None in practice; [`Factor::new`]'s, as the type requires.
57    pub fn factors() -> Result<Vec<Factor>, AnalysisError> {
58        ["x1", "x2", "x3"]
59            .into_iter()
60            .map(|name| Factor::new(name, -PI, PI))
61            .collect()
62    }
63
64    /// `y` at `x = [x₁, x₂, x₃]`; `NaN` for a slice of another length.
65    pub fn evaluate(&self, x: &[f64]) -> f64 {
66        let &[x1, x2, x3] = x else {
67            return f64::NAN;
68        };
69        let s2 = x2.sin();
70        x1.sin() + self.a * s2 * s2 + self.b * x3.powi(4) * x1.sin()
71    }
72
73    /// Its mean, `a/2`.
74    pub fn mean(&self) -> f64 {
75        0.5 * self.a
76    }
77
78    /// Its variance, `V`.
79    pub fn variance(&self) -> f64 {
80        let (a, b) = (self.a, self.b);
81        0.5 + a * a / 8.0 + b * PI.powi(4) / 5.0 + b * b * PI.powi(8) / 18.0
82    }
83
84    /// `V₁₃`, the part of the variance `x₁` and `x₃` cause together.
85    fn v13(&self) -> f64 {
86        8.0 * self.b * self.b * PI.powi(8) / 225.0
87    }
88
89    /// The first-order indices, `S₁`, `S₂`, `S₃`.
90    pub fn first_order(&self) -> [f64; 3] {
91        let v = self.variance();
92        let v1 = 0.5 * (1.0 + self.b * PI.powi(4) / 5.0).powi(2);
93        let v2 = self.a * self.a / 8.0;
94        [v1 / v, v2 / v, 0.0]
95    }
96
97    /// The total indices, `S_T₁`, `S_T₂`, `S_T₃`.
98    pub fn total(&self) -> [f64; 3] {
99        let [s1, s2, _] = self.first_order();
100        let s13 = self.v13() / self.variance();
101        [s1 + s13, s2, s13]
102    }
103}
104
105/// Sobol's g function, `y = Πᵢ gᵢ(xᵢ)` with `gᵢ = (|4xᵢ − 2| + aᵢ)/(1 + aᵢ)`, each `xᵢ` uniform
106/// on `[0, 1]`.
107///
108/// Each `gᵢ` has mean 1 and variance `Vᵢ = 1/(3 (1 + aᵢ)²)`, so (A. Saltelli and others, 2010,
109/// the module [`super::sobol`]'s reference, Appendix A, p. 268, eqs. (31) and (32))
110///
111/// - `V = Πᵢ (1 + Vᵢ) − 1`, computed as `exp(Σᵢ ln(1 + Vᵢ)) − 1` with `ln_1p` and `exp_m1`, so
112///   that a large `aᵢ` (a tiny `Vᵢ`) isn't lost in the `1 +`,
113/// - `Sᵢ = Vᵢ/V`, and
114/// - `S_Tᵢ = Vᵢ Πⱼ≠ᵢ (1 + Vⱼ) / V`.
115#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
116#[serde(try_from = "SobolGData")]
117pub struct SobolG {
118    a: Vec<f64>,
119}
120
121/// The serialized form of a [`SobolG`].
122#[derive(Deserialize)]
123#[serde(deny_unknown_fields)]
124struct SobolGData {
125    a: Vec<f64>,
126}
127
128impl TryFrom<SobolGData> for SobolG {
129    type Error = AnalysisError;
130
131    fn try_from(data: SobolGData) -> Result<Self, AnalysisError> {
132        Self::new(data.a)
133    }
134}
135
136impl SobolG {
137    /// The largest `aᵢ`: `3 (1 + aᵢ)²` overflows past about 7.7e153, and at 1e100 a factor's
138    /// share is already 1e-200 of the next one's.
139    pub const MAX_A: f64 = 1e100;
140
141    /// The function with one factor for each of `a`.
142    ///
143    /// # Errors
144    ///
145    /// - [`AnalysisError::TooFew`] for an empty `a`.
146    /// - [`AnalysisError::Domain`] for an `aᵢ` that isn't between 0 and [`SobolG::MAX_A`].
147    pub fn new(a: Vec<f64>) -> Result<Self, AnalysisError> {
148        if a.is_empty() {
149            return Err(AnalysisError::TooFew {
150                what: "g function's factors",
151                count: 0,
152                minimum: 1,
153            });
154        }
155        if let Some(&value) = a.iter().find(|&&ai| !(0.0..=Self::MAX_A).contains(&ai)) {
156            return Err(AnalysisError::Domain {
157                what: "g function's a (from 0 to 1e100)",
158                value,
159            });
160        }
161        Ok(Self { a })
162    }
163
164    /// The `aᵢ`.
165    pub fn a(&self) -> &[f64] {
166        &self.a
167    }
168
169    /// Its factors, `x1` to `xk`, each on `[0, 1]`.
170    ///
171    /// # Errors
172    ///
173    /// None in practice; [`Factor::new`]'s, as the type requires.
174    pub fn factors(&self) -> Result<Vec<Factor>, AnalysisError> {
175        (1..=self.a.len())
176            .map(|i| Factor::new(format!("x{i}"), 0.0, 1.0))
177            .collect()
178    }
179
180    /// `y` at `x`; `NaN` for a slice whose length isn't the number of factors.
181    pub fn evaluate(&self, x: &[f64]) -> f64 {
182        if x.len() != self.a.len() {
183            return f64::NAN;
184        }
185        x.iter()
186            .zip(&self.a)
187            .map(|(&xi, &ai)| ((4.0 * xi - 2.0).abs() + ai) / (1.0 + ai))
188            .product()
189    }
190
191    /// Its mean, 1.
192    pub fn mean(&self) -> f64 {
193        1.0
194    }
195
196    /// Each factor's part of the variance alone, `Vᵢ = 1/(3 (1 + aᵢ)²)`.
197    fn parts(&self) -> Vec<f64> {
198        self.a
199            .iter()
200            .map(|ai| 1.0 / (3.0 * (1.0 + ai) * (1.0 + ai)))
201            .collect()
202    }
203
204    /// Its variance, `V = Πᵢ (1 + Vᵢ) − 1`.
205    pub fn variance(&self) -> f64 {
206        self.log_product().exp_m1()
207    }
208
209    /// `ln Πᵢ (1 + Vᵢ)`.
210    fn log_product(&self) -> f64 {
211        self.parts().iter().map(|v| v.ln_1p()).sum()
212    }
213
214    /// The first-order indices, `Sᵢ = Vᵢ/V`.
215    pub fn first_order(&self) -> Vec<f64> {
216        let v = self.variance();
217        self.parts().iter().map(|vi| vi / v).collect()
218    }
219
220    /// The total indices, `S_Tᵢ = Vᵢ Πⱼ≠ᵢ (1 + Vⱼ) / V`.
221    pub fn total(&self) -> Vec<f64> {
222        let log_all = self.log_product();
223        let v = log_all.exp_m1();
224        self.parts()
225            .iter()
226            .map(|vi| vi * (log_all - vi.ln_1p()).exp() / v)
227            .collect()
228    }
229}
230
231#[cfg(test)]
232mod tests {
233    use super::*;
234
235    #[test]
236    fn ishigamis_indices_are_the_published_ones() {
237        // Sobol' and Levitan's b = 0.05 case (p. 59, Example 6.3): S₁ = 0.219, S₂ = 0.687,
238        // S₁₃ = 0.0946, to their printing.
239        let f = Ishigami { a: 7.0, b: 0.05 };
240        let [s1, s2, s3] = f.first_order();
241        let [t1, _, t3] = f.total();
242        assert!((s1 - 0.219).abs() < 5e-4, "{s1}");
243        assert!((s2 - 0.687).abs() < 5e-4, "{s2}");
244        assert_eq!(s3, 0.0);
245        assert!((t3 - 0.0946).abs() < 5e-5, "{t3}");
246        assert!((t1 - s1 - t3).abs() < 1e-15);
247        // a = 7, b = 0.1: S₁ = 0.3139, S₂ = 0.4424 (Marrel and others, 2009, quoting them).
248        let [s1, s2, _] = Ishigami::STANDARD.first_order();
249        assert!((s1 - 0.3139).abs() < 5e-5 && (s2 - 0.4424).abs() < 5e-5);
250        // The parts add up to the whole.
251        let f = Ishigami::STANDARD;
252        let [s1, s2, _] = f.first_order();
253        let [_, _, t3] = f.total();
254        assert!((s1 + s2 + t3 - 1.0).abs() < 1e-15);
255    }
256
257    #[test]
258    fn ishigamis_mean_and_variance_are_its_integrals() {
259        // Midpoint rule on a 120³ grid: smooth and periodic in x₁ and x₂, so it converges fast;
260        // x₃⁴ is not periodic, so allow its O(h²) error.
261        let f = Ishigami::STANDARD;
262        let n = 120;
263        let h = 2.0 * PI / f64::from(n);
264        let (mut sum, mut squares) = (0.0, 0.0);
265        for i in 0..n {
266            for j in 0..n {
267                for k in 0..n {
268                    let x = |m: i32| -PI + (f64::from(m) + 0.5) * h;
269                    let y = f.evaluate(&[x(i), x(j), x(k)]);
270                    sum += y;
271                    squares += y * y;
272                }
273            }
274        }
275        let count = f64::from(n).powi(3);
276        let mean = sum / count;
277        let variance = squares / count - mean * mean;
278        assert!((mean - f.mean()).abs() < 1e-9, "{mean}");
279        assert!(
280            (variance - f.variance()).abs() / f.variance() < 1e-3,
281            "{variance}"
282        );
283    }
284
285    #[test]
286    fn the_g_functions_indices_add_up_and_order_by_a() {
287        let g = SobolG::new(vec![0.0, 1.0, 4.5, 9.0, 99.0, 99.0, 99.0, 99.0]).unwrap();
288        let s = g.first_order();
289        let t = g.total();
290        assert!(s.iter().sum::<f64>() < 1.0);
291        assert!(t.iter().sum::<f64>() > 1.0);
292        for i in 0..8 {
293            assert!(s[i] <= t[i]);
294        }
295        for i in 0..3 {
296            assert!(t[i] > t[i + 1]);
297        }
298        // A huge a: its tiny part isn't lost in 1 + Vᵢ.
299        let tiny = SobolG::new(vec![1e9]).unwrap();
300        let part = 1.0 / (3.0 * (1.0 + 1e9) * (1.0 + 1e9));
301        assert!((tiny.variance() / part - 1.0).abs() < 1e-12);
302        assert!((tiny.first_order()[0] - 1.0).abs() < 1e-12);
303        // One factor: all its variance is its own, 1/3 at a = 0.
304        let one = SobolG::new(vec![0.0]).unwrap();
305        assert!((one.variance() - 1.0 / 3.0).abs() < 1e-15);
306        // V = exp_m1(ln_1p(1/3)) may differ from 1/3 in its last bits on another platform's
307        // ln_1p and exp_m1.
308        assert!((one.first_order()[0] - 1.0).abs() < 1e-14);
309        assert!((one.total()[0] - 1.0).abs() < 1e-14);
310    }
311
312    #[test]
313    fn the_g_functions_mean_and_variance_are_its_integrals() {
314        // Two factors on a 4000² midpoint grid: |4x − 2| is linear on each half, so the rule is
315        // exact on each cell but the kink's, which falls on a cell edge.
316        let g = SobolG::new(vec![0.5, 2.0]).unwrap();
317        let n = 4000;
318        let h = 1.0 / f64::from(n);
319        let (mut sum, mut squares) = (0.0, 0.0);
320        for i in 0..n {
321            for j in 0..n {
322                let y = g.evaluate(&[(f64::from(i) + 0.5) * h, (f64::from(j) + 0.5) * h]);
323                sum += y;
324                squares += y * y;
325            }
326        }
327        let count = f64::from(n * n);
328        let mean = sum / count;
329        let variance = squares / count - mean * mean;
330        assert!((mean - 1.0).abs() < 1e-12, "{mean}");
331        assert!((variance - g.variance()).abs() < 1e-7, "{variance}");
332    }
333
334    #[test]
335    fn a_point_of_the_wrong_length_is_nan_and_bad_parameters_are_refused() {
336        assert!(Ishigami::STANDARD.evaluate(&[0.0, 0.0]).is_nan());
337        let g = SobolG::new(vec![1.0, 2.0]).unwrap();
338        assert!(g.evaluate(&[0.5]).is_nan());
339        assert_eq!(g.evaluate(&[0.5, 0.5]), (1.0 / 2.0) * (2.0 / 3.0));
340        assert!(matches!(
341            SobolG::new(vec![]),
342            Err(AnalysisError::TooFew {
343                what: "g function's factors",
344                count: 0,
345                minimum: 1
346            })
347        ));
348        match SobolG::new(vec![1.0, -0.5]) {
349            Err(AnalysisError::Domain { value, .. }) => assert_eq!(value, -0.5),
350            other => panic!("{other:?}"),
351        }
352        let refused = serde_json::from_str::<SobolG>(r#"{"a":[1.0,-1.0]}"#)
353            .unwrap_err()
354            .to_string();
355        assert!(refused.contains("g function's a"), "{refused}");
356        assert!(SobolG::new(vec![SobolG::MAX_A]).is_ok());
357        match SobolG::new(vec![1e101]) {
358            Err(AnalysisError::Domain { what, value }) => {
359                assert_eq!((what, value), ("g function's a (from 0 to 1e100)", 1e101));
360            }
361            other => panic!("{other:?}"),
362        }
363        let big = SobolG::new(vec![SobolG::MAX_A, SobolG::MAX_A]).unwrap();
364        assert!(big.first_order().iter().all(|s| (s - 0.5).abs() < 1e-12));
365        assert_eq!(
366            serde_json::from_str::<SobolG>(r#"{"a":[1.0,2.0]}"#).unwrap(),
367            g
368        );
369    }
370}
hpr_analysis/sensitivity/benchmark.rs

hpr_analysis/sensitivity/
benchmark.rs
