ggplot_rs/
helpers.rs

1//! ggplot2's helper functions (#414): discretising a continuous variable
2//! (`cut_interval`, `cut_number`, `cut_width`), `resolution`, and the
3//! `stat_summary` summary functions (`mean_se`, `mean_cl_normal`,
4//! `mean_cl_boot`, `mean_sdl`, `median_hilow`) as plain functions.
5
6use crate::data::Value;
7use crate::stat::summary::SummaryData;
8
9/// A discretised variable: one interval label per input value (`Value::Na`
10/// for missing / non-finite input) and the interval labels in order — use
11/// `values` as a `group` / `fill` column and `levels` for the scale order.
12#[derive(Clone, Debug, PartialEq)]
13pub struct Cut {
14    pub values: Vec<Value>,
15    pub levels: Vec<String>,
16}
17
18impl From<Cut> for Vec<Value> {
19    fn from(c: Cut) -> Self {
20        c.values
21    }
22}
23
24/// C's `%.{digits}g` (what R's `formatC(x, digits)` prints).
25fn fmt_g(v: f64, digits: usize) -> String {
26    if v == 0.0 {
27        return "0".into();
28    }
29    let digits = digits.max(1);
30    let sci = format!("{:.*e}", digits - 1, v);
31    let (mant, exp) = sci.split_once('e').unwrap();
32    let exp: i32 = exp.parse().unwrap();
33    let trim = |s: &str| -> String {
34        if s.contains('.') {
35            s.trim_end_matches('0').trim_end_matches('.').to_string()
36        } else {
37            s.to_string()
38        }
39    };
40    if exp < -4 || exp >= digits as i32 {
41        let sign = if exp < 0 { '-' } else { '+' };
42        format!("{}e{sign}{:02}", trim(mant), exp.abs())
43    } else {
44        let decimals = (digits as i32 - 1 - exp).max(0) as usize;
45        trim(&format!("{v:.decimals$}"))
46    }
47}
48
49/// R's `cut(x, breaks, include.lowest = TRUE, right)` with its default
50/// `dig.lab = 3` labels.
51fn cut(x: &[f64], breaks: &[f64], right: bool) -> Cut {
52    let nb = breaks.len();
53    if nb < 2 {
54        return Cut {
55            values: vec![Value::Na; x.len()],
56            levels: Vec::new(),
57        };
58    }
59    // R: widen dig.lab until neighbouring break labels differ.
60    let mut labels: Vec<String> = Vec::new();
61    for dig in 3..=12 {
62        labels = breaks.iter().map(|b| fmt_g(*b, dig)).collect();
63        if labels.windows(2).all(|w| w[0] != w[1]) {
64            break;
65        }
66    }
67    let levels: Vec<String> = (0..nb - 1)
68        .map(|i| {
69            let (lo, hi) = (&labels[i], &labels[i + 1]);
70            match (right, i == 0, i == nb - 2) {
71                (true, true, _) => format!("[{lo},{hi}]"),
72                (true, false, _) => format!("({lo},{hi}]"),
73                (false, _, true) => format!("[{lo},{hi}]"),
74                (false, _, false) => format!("[{lo},{hi})"),
75            }
76        })
77        .collect();
78    let values = x
79        .iter()
80        .map(|&v| {
81            if !v.is_finite() {
82                return Value::Na;
83            }
84            let k = (0..nb - 1).find(|&i| {
85                let (lo, hi) = (breaks[i], breaks[i + 1]);
86                if right {
87                    (v > lo || (i == 0 && v >= lo)) && v <= hi
88                } else {
89                    v >= lo && (v < hi || (i == nb - 2 && v <= hi))
90                }
91            });
92            k.map_or(Value::Na, |k| Value::Str(levels[k].clone()))
93        })
94        .collect();
95    Cut { values, levels }
96}
97
98fn finite_range(x: &[f64]) -> Option<(f64, f64)> {
99    let f = x.iter().copied().filter(|v| v.is_finite());
100    let lo = f.clone().fold(f64::INFINITY, f64::min);
101    let hi = f.fold(f64::NEG_INFINITY, f64::max);
102    lo.is_finite().then_some((lo, hi))
103}
104
105/// `n` intervals of equal width over the range of `x` — ggplot2's
106/// `cut_interval(x, n)`.
107pub fn cut_interval(x: &[f64], n: usize) -> Cut {
108    let Some((lo, hi)) = finite_range(x) else {
109        return cut(x, &[], true);
110    };
111    let n = n.max(1);
112    let breaks: Vec<f64> = (0..=n)
113        .map(|i| lo + (hi - lo) * i as f64 / n as f64)
114        .collect();
115    cut(x, &breaks, true)
116}
117
118/// Intervals of width `length` on a grid of multiples of `length` covering
119/// `x` — ggplot2's `cut_interval(x, length = …)`.
120pub fn cut_interval_length(x: &[f64], length: f64) -> Cut {
121    let Some((lo, hi)) = finite_range(x) else {
122        return cut(x, &[], true);
123    };
124    if !(length.is_finite() && length > 0.0) {
125        return cut(x, &[], true);
126    }
127    // scales::fullseq
128    let breaks: Vec<f64> = if hi - lo < 1e-12 {
129        vec![lo - length / 2.0, lo + length / 2.0]
130    } else {
131        let (a, b) = ((lo / length).floor(), (hi / length).ceil());
132        (0..=(b - a) as i64)
133            .map(|k| (a + k as f64) * length)
134            .collect()
135    };
136    cut(x, &breaks, true)
137}
138
139/// `n` intervals with (about) equal numbers of observations — ggplot2's
140/// `cut_number(x, n)` (R type-7 quantile breaks).
141pub fn cut_number(x: &[f64], n: usize) -> Cut {
142    let mut s: Vec<f64> = x.iter().copied().filter(|v| v.is_finite()).collect();
143    if s.is_empty() {
144        return cut(x, &[], true);
145    }
146    s.sort_by(f64::total_cmp);
147    let n = n.max(1);
148    let q = |p: f64| {
149        let h = (s.len() - 1) as f64 * p;
150        let (a, b) = (h.floor() as usize, h.ceil() as usize);
151        s[a] + (h - a as f64) * (s[b] - s[a])
152    };
153    let breaks: Vec<f64> = (0..=n).map(|i| q(i as f64 / n as f64)).collect();
154    cut(x, &breaks, true)
155}
156
157/// Options for [`cut_width`]: at most one of `center` / `boundary`, and
158/// which side of each interval is closed.
159#[derive(Clone, Copy, Debug, Default, PartialEq)]
160pub struct CutWidth {
161    /// Centre of one bin (ggplot2's `center`).
162    pub center: Option<f64>,
163    /// An edge of one bin (ggplot2's `boundary`; default `width / 2`).
164    pub boundary: Option<f64>,
165    /// Close intervals on the left (`closed = "left"`); default right.
166    pub closed_left: bool,
167}
168
169/// Intervals of width `width` — ggplot2's `cut_width(x, width, center,
170/// boundary, closed)`.
171pub fn cut_width(x: &[f64], width: f64, opts: CutWidth) -> Cut {
172    let Some((lo, _)) = finite_range(x) else {
173        return cut(x, &[], true);
174    };
175    if !(width.is_finite() && width > 0.0) {
176        return cut(x, &[], true);
177    }
178    let boundary = match (opts.boundary, opts.center) {
179        (Some(b), _) => b,
180        (None, Some(c)) => c - width / 2.0,
181        (None, None) => width / 2.0,
182    };
183    let min_x = boundary + ((lo - boundary) / width).floor() * width;
184    let max_x = x
185        .iter()
186        .copied()
187        .filter(|v| v.is_finite())
188        .fold(f64::NEG_INFINITY, f64::max)
189        + (1.0 - 1e-8) * width;
190    let mut breaks = Vec::new();
191    let mut b = min_x;
192    while b <= max_x + 1e-12 {
193        breaks.push(b);
194        b = min_x + breaks.len() as f64 * width;
195    }
196    cut(x, &breaks, !opts.closed_left)
197}
198
199/// The smallest non-zero gap between the (unique) values of `x` — ggplot2's
200/// `resolution(x, zero)`; `zero = true` also counts the distance to 0.
201/// Constant data → 1.
202pub fn resolution(x: &[f64], zero: bool) -> f64 {
203    let mut v: Vec<f64> = x.iter().copied().filter(|v| v.is_finite()).collect();
204    let Some((lo, hi)) = finite_range(&v) else {
205        return 1.0;
206    };
207    if hi - lo < 1e-12 * hi.abs().max(1.0) {
208        return 1.0;
209    }
210    if zero {
211        v.push(0.0);
212    }
213    v.sort_by(f64::total_cmp);
214    v.dedup();
215    let tol = f64::EPSILON.sqrt();
216    v.windows(2)
217        .map(|w| w[1] - w[0])
218        .filter(|d| *d > tol)
219        .fold(f64::INFINITY, f64::min)
220        .min(f64::MAX)
221}
222
223/// `(y, ymin, ymax)` = mean ± `mult` standard errors (`mean_se(x, mult)`).
224pub fn mean_se(x: &[f64], mult: f64) -> (f64, f64, f64) {
225    let (m, lo, _) = SummaryData::MeanSe.apply3(x);
226    let se = (m - lo) * mult;
227    (m, m - se, m + se)
228}
229
230/// Mean with a Student-t confidence interval (`mean_cl_normal(x, conf.int)`).
231pub fn mean_cl_normal(x: &[f64], conf_int: f64) -> (f64, f64, f64) {
232    SummaryData::MeanClNormal { level: conf_int }.apply3(x)
233}
234
235/// Mean with a percentile-bootstrap interval from `b` resamples
236/// (`mean_cl_boot(x, conf.int, B)`; seeded, so reproducible — but not R's
237/// random stream).
238pub fn mean_cl_boot(x: &[f64], conf_int: f64, b: usize) -> (f64, f64, f64) {
239    SummaryData::MeanClBoot { level: conf_int, b }.apply3(x)
240}
241
242/// Mean ± `mult` standard deviations (`mean_sdl(x, mult)`, default 2).
243pub fn mean_sdl(x: &[f64], mult: f64) -> (f64, f64, f64) {
244    SummaryData::MeanSdl { mult }.apply3(x)
245}
246
247/// Median with the outer `(1 ± conf_int) / 2` sample quantiles
248/// (`median_hilow(x, conf.int)`).
249pub fn median_hilow(x: &[f64], conf_int: f64) -> (f64, f64, f64) {
250    SummaryData::MedianHilow { level: conf_int }.apply3(x)
251}
252
253#[cfg(test)]
254mod tests {
255    use super::*;
256
257    const X: [f64; 10] = [0.5, 1.2, 3.3, 4.7, 5.5, 7.25, 8.1, 9.9, 10.0, f64::NAN];
258
259    fn labels(c: &Cut) -> Vec<String> {
260        c.values
261            .iter()
262            .map(|v| v.as_str().map_or("NA".to_string(), str::to_string))
263            .collect()
264    }
265
266    fn r(s: &str) -> Vec<String> {
267        s.split(", ")
268            .map(|t| t.trim_matches('"').to_string())
269            .collect()
270    }
271
272    // Expected values: ggplot2 4.0.3 / R 4.6.0, printed with as.character().
273    #[test]
274    fn cut_labels_match_ggplot2() {
275        assert_eq!(
276            labels(&cut_interval(&X, 3)),
277            r(
278                r#""[0.5,3.67]", "[0.5,3.67]", "[0.5,3.67]", "(3.67,6.83]", "(3.67,6.83]", "(6.83,10]", "(6.83,10]", "(6.83,10]", "(6.83,10]", "NA""#
279            )
280        );
281        assert_eq!(
282            labels(&cut_interval_length(&X, 2.5)),
283            r(
284                r#""[0,2.5]", "[0,2.5]", "(2.5,5]", "(2.5,5]", "(5,7.5]", "(5,7.5]", "(7.5,10]", "(7.5,10]", "(7.5,10]", "NA""#
285            )
286        );
287        assert_eq!(
288            labels(&cut_number(&X, 3)),
289            r(
290                r#""[0.5,4.23]", "[0.5,4.23]", "[0.5,4.23]", "(4.23,7.53]", "(4.23,7.53]", "(4.23,7.53]", "(7.53,10]", "(7.53,10]", "(7.53,10]", "NA""#
291            )
292        );
293        assert_eq!(
294            labels(&cut_width(&X, 2.0, CutWidth::default())),
295            r(
296                r#""[-1,1]", "(1,3]", "(3,5]", "(3,5]", "(5,7]", "(7,9]", "(7,9]", "(9,11]", "(9,11]", "NA""#
297            )
298        );
299        let b0 = CutWidth {
300            boundary: Some(0.0),
301            ..Default::default()
302        };
303        assert_eq!(
304            labels(&cut_width(&X, 2.0, b0)),
305            r(
306                r#""[0,2]", "[0,2]", "(2,4]", "(4,6]", "(4,6]", "(6,8]", "(8,10]", "(8,10]", "(8,10]", "NA""#
307            )
308        );
309        let left = CutWidth {
310            center: Some(0.0),
311            closed_left: true,
312            ..Default::default()
313        };
314        assert_eq!(
315            labels(&cut_width(&X, 2.0, left)),
316            r(
317                r#""[-1,1)", "[1,3)", "[3,5)", "[3,5)", "[5,7)", "[7,9)", "[7,9)", "[9,11]", "[9,11]", "NA""#
318            )
319        );
320        assert_eq!(
321            labels(&cut_interval(&[1000.0, 2500.0, 12345.0], 2)),
322            r(r#""[1e+03,6.67e+03]", "[1e+03,6.67e+03]", "(6.67e+03,1.23e+04]""#)
323        );
324        assert_eq!(
325            cut_interval(&X, 3).levels,
326            ["[0.5,3.67]", "(3.67,6.83]", "(6.83,10]"]
327        );
328    }
329
330    #[test]
331    fn resolution_matches_ggplot2() {
332        assert_eq!(resolution(&[1.0, 1.5, 3.0, 4.0], true), 0.5);
333        assert_eq!(resolution(&[2.5, 3.5, 4.75], true), 1.0);
334        assert_eq!(resolution(&[1.0, 2.0, 3.0, 4.0, 5.0], true), 1.0);
335        assert_eq!(resolution(&[10.0, 20.0, 30.0], false), 10.0);
336    }
337
338    #[test]
339    fn summary_functions_match_ggplot2() {
340        let y = [2.1, 3.4, 1.9, 5.6, 4.4, 3.3, 2.8];
341        let close = |got: (f64, f64, f64), want: (f64, f64, f64)| {
342            for (a, b) in [(got.0, want.0), (got.1, want.1), (got.2, want.2)] {
343                assert!((a - b).abs() < 1e-5, "{got:?} vs {want:?}");
344            }
345        };
346        close(mean_se(&y, 1.0), (3.357143, 2.865858, 3.848427));
347        close(mean_se(&y, 2.0), (3.357143, 2.374574, 4.339712));
348        close(mean_cl_normal(&y, 0.95), (3.357143, 2.155013, 4.559273));
349        close(mean_cl_normal(&y, 0.9), (3.357143, 2.402488, 4.311797));
350        close(mean_sdl(&y, 2.0), (3.357143, 0.7575092, 5.956777));
351        close(mean_sdl(&y, 1.0), (3.357143, 2.057326, 4.65696));
352        close(median_hilow(&y, 0.95), (3.3, 1.93, 5.42));
353        close(median_hilow(&y, 0.5), (3.3, 2.45, 3.9));
354        let (m, lo, hi) = mean_cl_boot(&y, 0.95, 1000);
355        assert!((m - 3.357143).abs() < 1e-5 && lo < m && hi > m);
356    }
357}
ggplot_rs/helpers.rs

ggplot_rs/
helpers.rs
