hpr_analysis/
sensitivity.rs

1//! Sensitivity analysis: which of a model's uncertain inputs move its output most.
2//!
3//! **Guide:** [Sensitivity analysis][guide] runs both methods on test functions whose answers are
4//! known, then on a rocket's apogee, and says how far to trust them.
5//!
6//! [guide]: https://nrdptel.github.io/hpr-sim/sensitivity.html
7//!
8//! Two methods, both over *factors*: inputs each spread evenly between a low and a high value
9//! ([`Factor`]), independent of each other.
10//!
11//! - [`morris`]: Morris's screening. A few cheap runs, stepping one factor at a time along random
12//!   paths through a grid, rank the factors by how much a step in each moves the output. It
13//!   finds the factors that don't matter, at about ten runs per factor.
14//! - [`sobol`]: Sobol' indices. Many runs split the output's variance among the factors: the
15//!   share each causes alone (its *first-order* index) and the share it has any part in, with
16//!   the others (its *total* index). It costs thousands of runs per factor.
17//! - [`benchmark`]: two test functions whose indices are known in closed form, Ishigami and
18//!   Homma's and Sobol's g, which the tests hold both methods to.
19//!
20//! Each method lays out its points first (a *design*), takes the model's output at each, in the
21//! design's order, and then analyses them. The model can be anything: a flight, flown with the
22//! factors set into its inputs, or a function. So the runs can be made however the caller likes,
23//! on many threads or on many machines; each method also has a shortcut that takes a closure.
24//!
25//! # Reproducibility
26//!
27//! A design is drawn from a seed. Each Morris path and each Sobol' sample row has its own random
28//! stream ([`SeededRng::for_stream`](hpr_core::random::SeededRng::for_stream)), keyed by the
29//! seed and its index, so path or row `k` is the same however many are drawn. Every sum runs in
30//! the design's order, so on one platform an analysis is bit for bit the same every time.
31//!
32//! # Left out
33//!
34//! Factors are uniform and independent: no other distributions, no correlations. A normal input
35//! can be given as a range about its mean (say ±2 standard deviations), which spreads it more
36//! evenly than it is. Sobol' samples are pseudo-random, not quasi-random. Second-order Sobol'
37//! indices, and Campolongo's choice of the most spread-out Morris paths among many, are not
38//! computed.
39
40pub mod benchmark;
41pub mod morris;
42pub mod sobol;
43
44use std::collections::BTreeSet;
45
46use serde::{Deserialize, Serialize};
47
48use crate::error::AnalysisError;
49
50/// The most points a design lays out: 1,048,576, far more than flights can be flown for. With
51/// [`MAX_DESIGN_VALUES`] it keeps a design's points under about 200 MiB (a point costs its
52/// coordinates, 8 bytes each, and a list of them, about 40 bytes).
53pub const MAX_DESIGN_POINTS: usize = 1 << 20;
54
55/// The most coordinates (points times factors) a design holds: 16,777,216, 128 MiB.
56pub const MAX_DESIGN_VALUES: usize = 1 << 24;
57
58/// An uncertain input, spread evenly between `low` and `high`. It serializes as its three
59/// fields, and reads back through [`Factor::new`]'s checks.
60#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
61#[serde(try_from = "FactorData")]
62pub struct Factor {
63    name: String,
64    low: f64,
65    high: f64,
66}
67
68/// The serialized form of a [`Factor`].
69#[derive(Deserialize)]
70#[serde(deny_unknown_fields)]
71struct FactorData {
72    name: String,
73    low: f64,
74    high: f64,
75}
76
77impl TryFrom<FactorData> for Factor {
78    type Error = AnalysisError;
79
80    fn try_from(data: FactorData) -> Result<Self, AnalysisError> {
81        Self::new(data.name, data.low, data.high)
82    }
83}
84
85impl Factor {
86    /// A factor named `name`, spread evenly over `[low, high]`.
87    ///
88    /// # Errors
89    ///
90    /// [`AnalysisError::Domain`] for a bound that isn't finite, or a `high` not above `low`.
91    pub fn new(name: impl Into<String>, low: f64, high: f64) -> Result<Self, AnalysisError> {
92        if !low.is_finite() {
93            return Err(AnalysisError::Domain {
94                what: "factor's low value",
95                value: low,
96            });
97        }
98        if !(high.is_finite() && high > low) {
99            return Err(AnalysisError::Domain {
100                what: "factor's high value (finite, above the low one)",
101                value: high,
102            });
103        }
104        Ok(Self {
105            name: name.into(),
106            low,
107            high,
108        })
109    }
110
111    /// Its name.
112    pub fn name(&self) -> &str {
113        &self.name
114    }
115
116    /// The low end of its range.
117    pub fn low(&self) -> f64 {
118        self.low
119    }
120
121    /// The high end of its range.
122    pub fn high(&self) -> f64 {
123        self.high
124    }
125
126    /// The value a fraction `u` of the way from `low` to `high`: `low + u (high − low)`.
127    pub fn at(&self, u: f64) -> f64 {
128        self.low + u * (self.high - self.low)
129    }
130}
131
132/// Checks that there is at least one factor, and no two share a name.
133fn check_factors(factors: &[Factor]) -> Result<(), AnalysisError> {
134    if factors.is_empty() {
135        return Err(AnalysisError::TooFew {
136            what: "factors",
137            count: 0,
138            minimum: 1,
139        });
140    }
141    let mut names = BTreeSet::new();
142    for factor in factors {
143        if !names.insert(factor.name()) {
144            return Err(AnalysisError::DuplicateFactor(factor.name().to_owned()));
145        }
146    }
147    Ok(())
148}
149
150/// Checks that `outputs` has one finite value for each of a design's `points`.
151fn check_outputs(outputs: &[f64], points: usize) -> Result<(), AnalysisError> {
152    if outputs.len() != points {
153        return Err(AnalysisError::Length {
154            what: "outputs, against the design's points",
155            length: outputs.len(),
156            expected: points,
157        });
158    }
159    if let Some((index, &value)) = outputs.iter().enumerate().find(|(_, y)| !y.is_finite()) {
160        return Err(AnalysisError::Output { index, value });
161    }
162    Ok(())
163}
164
165/// Checks a design of `count` paths or rows, `per` points each, over `factors` factors, against
166/// [`MAX_DESIGN_POINTS`] and [`MAX_DESIGN_VALUES`].
167fn check_size(
168    what: &'static str,
169    count: usize,
170    per: usize,
171    factors: usize,
172) -> Result<(), AnalysisError> {
173    let points = count.saturating_mul(per);
174    if points > MAX_DESIGN_POINTS {
175        return Err(AnalysisError::Count {
176            what,
177            count: points,
178            limit: MAX_DESIGN_POINTS,
179        });
180    }
181    let values = points.saturating_mul(factors);
182    if values > MAX_DESIGN_VALUES {
183        return Err(AnalysisError::Count {
184            what: "design coordinates (points times factors)",
185            count: values,
186            limit: MAX_DESIGN_VALUES,
187        });
188    }
189    Ok(())
190}
191
192#[cfg(test)]
193mod tests {
194    use super::*;
195
196    #[test]
197    fn a_factor_maps_the_unit_interval_onto_its_range() {
198        let f = Factor::new("drag", 0.9, 1.1).unwrap();
199        assert_eq!(f.at(0.0), 0.9);
200        assert!((f.at(1.0) - 1.1).abs() < 1e-15);
201        assert!((f.at(0.5) - 1.0).abs() < 1e-15);
202    }
203
204    #[test]
205    fn a_factor_refuses_an_empty_or_reversed_range() {
206        let refused = |low: f64, high: f64| match Factor::new("x", low, high) {
207            Err(AnalysisError::Domain { what, value }) => (what, value),
208            other => panic!("{low}, {high}: {other:?}"),
209        };
210        assert_eq!(refused(f64::NAN, 1.0).0, "factor's low value");
211        let (what, value) = refused(1.0, 1.0);
212        assert_eq!(what, "factor's high value (finite, above the low one)");
213        assert_eq!(value, 1.0);
214        assert_eq!(refused(1.0, 0.0).1, 0.0);
215        assert!(refused(0.0, f64::INFINITY).1.is_infinite());
216    }
217
218    #[test]
219    fn a_factor_reads_back_through_its_checks() {
220        let f = Factor::new("wind", 0.0, 8.0).unwrap();
221        let json = serde_json::to_string(&f).unwrap();
222        assert_eq!(json, r#"{"name":"wind","low":0.0,"high":8.0}"#);
223        assert_eq!(serde_json::from_str::<Factor>(&json).unwrap(), f);
224        let refused = serde_json::from_str::<Factor>(r#"{"name":"w","low":2.0,"high":1.0}"#)
225            .unwrap_err()
226            .to_string();
227        assert!(refused.contains("factor's high value"), "{refused}");
228        let unknown = serde_json::from_str::<Factor>(r#"{"name":"w","low":0.0,"high":1.0,"sd":1}"#)
229            .unwrap_err()
230            .to_string();
231        assert!(unknown.contains("unknown field `sd`"), "{unknown}");
232    }
233
234    #[test]
235    fn factors_need_one_and_distinct_names() {
236        assert!(matches!(
237            check_factors(&[]),
238            Err(AnalysisError::TooFew {
239                what: "factors",
240                count: 0,
241                minimum: 1
242            })
243        ));
244        let x = Factor::new("x", 0.0, 1.0).unwrap();
245        let y = Factor::new("y", 0.0, 1.0).unwrap();
246        assert!(check_factors(&[x.clone(), y.clone()]).is_ok());
247        match check_factors(&[x.clone(), y, x]) {
248            Err(AnalysisError::DuplicateFactor(name)) => assert_eq!(name, "x"),
249            other => panic!("{other:?}"),
250        }
251    }
252
253    #[test]
254    fn a_design_is_held_to_its_size_limits() {
255        assert!(check_size("points", 1 << 18, 4, 16).is_ok());
256        match check_size("points", 1 << 19, 3, 1) {
257            Err(AnalysisError::Count { what, count, limit }) => {
258                assert_eq!((what, count, limit), ("points", 3 << 19, MAX_DESIGN_POINTS));
259            }
260            other => panic!("{other:?}"),
261        }
262        match check_size("points", usize::MAX, 2, 1) {
263            Err(AnalysisError::Count { count, .. }) => assert_eq!(count, usize::MAX),
264            other => panic!("{other:?}"),
265        }
266        match check_size("points", 1 << 18, 4, 17) {
267            Err(AnalysisError::Count { what, count, limit }) => {
268                assert_eq!(what, "design coordinates (points times factors)");
269                assert_eq!((count, limit), (17 << 20, MAX_DESIGN_VALUES));
270            }
271            other => panic!("{other:?}"),
272        }
273    }
274
275    #[test]
276    fn outputs_must_match_the_design_and_be_finite() {
277        assert!(check_outputs(&[1.0, 2.0], 2).is_ok());
278        match check_outputs(&[1.0], 2) {
279            Err(AnalysisError::Length {
280                length, expected, ..
281            }) => assert_eq!((length, expected), (1, 2)),
282            other => panic!("{other:?}"),
283        }
284        match check_outputs(&[1.0, f64::NAN, f64::INFINITY], 3) {
285            Err(AnalysisError::Output { index, value }) => {
286                assert_eq!(index, 1);
287                assert!(value.is_nan());
288            }
289            other => panic!("{other:?}"),
290        }
291    }
292}
hpr_analysis/sensitivity.rs

hpr_analysis/
sensitivity.rs
