ggplot_rs/scale/
transform.rs

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

ggplot_rs/scale/
transform.rs
