ggplot_rs/scale/
transform.rs

1use crate::data::Value;
2
3/// Scale transformation types.
4#[derive(Clone, Debug)]
5pub enum ScaleTransform {
6    Identity,
7    Log10,
8    Log2,
9    Ln,
10    Sqrt,
11    Reverse,
12    /// Logit: `ln(p / (1 - p))`, for proportions in (0, 1).
13    Logit,
14    /// Probit: the inverse normal CDF, for proportions in (0, 1).
15    Probit,
16    /// Sign-preserving pseudo-log (`asinh(x / 2)`) — handles zero and negatives.
17    PseudoLog,
18    /// Reciprocal `1 / x`.
19    Reciprocal,
20    /// Exponential (`exp`); axis labels are therefore spaced logarithmically.
21    Exp,
22    /// Box–Cox with the given lambda: `(x^λ − 1) / λ`, or `ln(x)` at λ = 0 (x > 0).
23    BoxCox(f64),
24}
25
26impl ScaleTransform {
27    /// Apply the forward transformation.
28    pub fn apply(&self, value: f64) -> f64 {
29        match self {
30            ScaleTransform::Identity => value,
31            ScaleTransform::Log10 => {
32                if value > 0.0 {
33                    value.log10()
34                } else {
35                    f64::NEG_INFINITY
36                }
37            }
38            ScaleTransform::Log2 => {
39                if value > 0.0 {
40                    value.log2()
41                } else {
42                    f64::NEG_INFINITY
43                }
44            }
45            ScaleTransform::Ln => {
46                if value > 0.0 {
47                    value.ln()
48                } else {
49                    f64::NEG_INFINITY
50                }
51            }
52            ScaleTransform::Sqrt => {
53                if value >= 0.0 {
54                    value.sqrt()
55                } else {
56                    f64::NAN
57                }
58            }
59            ScaleTransform::Reverse => -value,
60            ScaleTransform::Logit => {
61                if value <= 0.0 {
62                    f64::NEG_INFINITY
63                } else if value >= 1.0 {
64                    f64::INFINITY
65                } else {
66                    (value / (1.0 - value)).ln()
67                }
68            }
69            ScaleTransform::Probit => qnorm(value),
70            ScaleTransform::PseudoLog => (value / 2.0).asinh(),
71            ScaleTransform::Reciprocal => {
72                if value != 0.0 {
73                    1.0 / value
74                } else {
75                    f64::NAN
76                }
77            }
78            ScaleTransform::Exp => value.exp(),
79            ScaleTransform::BoxCox(lambda) => {
80                if value <= 0.0 {
81                    f64::NAN
82                } else if lambda.abs() < 1e-9 {
83                    value.ln()
84                } else {
85                    (value.powf(*lambda) - 1.0) / lambda
86                }
87            }
88        }
89    }
90
91    /// Apply the inverse transformation.
92    pub fn inverse(&self, value: f64) -> f64 {
93        match self {
94            ScaleTransform::Identity => value,
95            ScaleTransform::Log10 => 10f64.powf(value),
96            ScaleTransform::Log2 => 2f64.powf(value),
97            ScaleTransform::Ln => value.exp(),
98            ScaleTransform::Sqrt => value * value,
99            ScaleTransform::Reverse => -value,
100            ScaleTransform::Logit => 1.0 / (1.0 + (-value).exp()),
101            ScaleTransform::Probit => pnorm(value),
102            ScaleTransform::PseudoLog => 2.0 * value.sinh(),
103            ScaleTransform::Reciprocal => {
104                if value != 0.0 {
105                    1.0 / value
106                } else {
107                    f64::NAN
108                }
109            }
110            ScaleTransform::Exp => value.ln(),
111            ScaleTransform::BoxCox(lambda) => {
112                if lambda.abs() < 1e-9 {
113                    value.exp()
114                } else {
115                    (value * lambda + 1.0).powf(1.0 / lambda)
116                }
117            }
118        }
119    }
120
121    /// Transform a Value.
122    pub fn transform_value(&self, value: &Value) -> Value {
123        match value.as_f64() {
124            Some(f) => {
125                let t = self.apply(f);
126                // ±Inf in → ±Inf out is kept: it means "panel edge" (geom_rect).
127                if t.is_finite() || (f.is_infinite() && t.is_infinite()) {
128                    Value::Float(t)
129                } else {
130                    Value::Na
131                }
132            }
133            None => value.clone(),
134        }
135    }
136
137    pub fn is_identity(&self) -> bool {
138        matches!(self, ScaleTransform::Identity)
139    }
140
141    /// "Nice" tick positions in *data* (pre-transform) space — used for
142    /// log/sqrt axes, where ticks should read as round numbers in the
143    /// original units rather than evenly spaced in transformed space (which
144    /// gives half-decades like `316.23` under log10, or squares like `625`
145    /// under sqrt). Mirrors ggplot2/scales: `trans$breaks(trans$inverse(limits))`;
146    /// the caller forward-transforms the result back onto the scale
147    /// (`ScaleContinuous::breaks`).
148    ///
149    /// `data_min`/`data_max` must be finite with `data_max > data_min`.
150    pub(crate) fn breaks_in_data_space(&self, data_min: f64, data_max: f64, n: usize) -> Vec<f64> {
151        match self {
152            ScaleTransform::Log10 => super::util::log_breaks(data_min, data_max, n, 10.0),
153            ScaleTransform::Log2 => super::util::log_breaks(data_min, data_max, n, 2.0),
154            ScaleTransform::Ln => {
155                super::util::log_breaks(data_min, data_max, n, std::f64::consts::E)
156            }
157            // sqrt_trans breaks in `scales` are extended_breaks() over the
158            // untransformed data — not evenly spaced in sqrt space.
159            ScaleTransform::Sqrt => super::util::extended_breaks(data_min, data_max, n),
160            _ => super::util::extended_breaks(data_min, data_max, n),
161        }
162    }
163
164    /// Whether [`breaks_in_data_space`](Self::breaks_in_data_space) applies —
165    /// i.e. this transform's ticks should be chosen in data space rather than
166    /// transformed space. Only log/sqrt change behaviour; other transforms
167    /// keep their previous (transformed-space) break placement.
168    pub(crate) fn uses_data_space_breaks(&self) -> bool {
169        matches!(
170            self,
171            ScaleTransform::Log10
172                | ScaleTransform::Log2
173                | ScaleTransform::Ln
174                | ScaleTransform::Sqrt
175        )
176    }
177}
178
179/// Standard normal CDF, via the Abramowitz & Stegun 7.1.26 erf approximation.
180fn pnorm(x: f64) -> f64 {
181    0.5 * (1.0 + erf(x / std::f64::consts::SQRT_2))
182}
183
184fn erf(x: f64) -> f64 {
185    let sign = if x < 0.0 { -1.0 } else { 1.0 };
186    let x = x.abs();
187    let t = 1.0 / (1.0 + 0.327_591_1 * x);
188    let y = 1.0
189        - (((((1.061_405_429 * t - 1.453_152_027) * t) + 1.421_413_741) * t - 0.284_496_736) * t
190            + 0.254_829_592)
191            * t
192            * (-x * x).exp();
193    sign * y
194}
195
196/// Inverse normal CDF (probit), Abramowitz & Stegun 26.2.23 rational approximation.
197fn qnorm(p: f64) -> f64 {
198    if p <= 0.0 {
199        return f64::NEG_INFINITY;
200    }
201    if p >= 1.0 {
202        return f64::INFINITY;
203    }
204    if p < 0.5 {
205        -rational_approx((-2.0 * p.ln()).sqrt())
206    } else if p > 0.5 {
207        rational_approx((-2.0 * (1.0 - p).ln()).sqrt())
208    } else {
209        0.0
210    }
211}
212
213fn rational_approx(t: f64) -> f64 {
214    let c0 = 2.515_517;
215    let c1 = 0.802_853;
216    let c2 = 0.010_328;
217    let d1 = 1.432_788;
218    let d2 = 0.189_269;
219    let d3 = 0.001_308;
220    t - (c0 + c1 * t + c2 * t * t) / (1.0 + d1 * t + d2 * t * t + d3 * t * t * t)
221}
222
223#[cfg(test)]
224mod tests {
225    use super::ScaleTransform::*;
226
227    fn roundtrip(t: super::ScaleTransform, v: f64, tol: f64) {
228        let back = t.inverse(t.apply(v));
229        assert!(
230            (back - v).abs() < tol,
231            "{t:?}: {v} -> {} -> {back}",
232            t.apply(v)
233        );
234    }
235
236    #[test]
237    fn transforms_roundtrip() {
238        roundtrip(Logit, 0.3, 1e-9);
239        roundtrip(Probit, 0.3, 1e-2); // approximation
240        roundtrip(PseudoLog, -4.0, 1e-9);
241        roundtrip(PseudoLog, 0.0, 1e-9);
242        roundtrip(Reciprocal, 2.5, 1e-9);
243        roundtrip(Exp, 1.7, 1e-9);
244        roundtrip(BoxCox(0.5), 4.0, 1e-9);
245        roundtrip(BoxCox(0.0), 4.0, 1e-9); // lambda 0 == ln
246    }
247
248    #[test]
249    fn transforms_domain_edges() {
250        assert_eq!(Logit.apply(0.0), f64::NEG_INFINITY);
251        assert_eq!(Logit.apply(1.0), f64::INFINITY);
252        assert!(Reciprocal.apply(0.0).is_nan());
253        assert!(BoxCox(0.5).apply(-1.0).is_nan());
254        assert_eq!(Probit.apply(0.5), 0.0);
255        // Box-Cox at lambda 0 is ln.
256        assert!((BoxCox(0.0).apply(std::f64::consts::E) - 1.0).abs() < 1e-9);
257    }
258}
ggplot_rs/scale/transform.rs

ggplot_rs/scale/
transform.rs
