Skip to main content

flow_gate_core/transform/
hyperlog.rs

1use crate::error::FlowGateError;
2use crate::traits::Transform;
3
4const TAYLOR_LENGTH: usize = 16;
5const LN_10: f64 = std::f64::consts::LN_10;
6
7#[derive(Debug, Clone, Copy, PartialEq)]
8pub struct HyperlogTransform {
9    pub t: f64,
10    pub w: f64,
11    pub m: f64,
12    pub a: f64,
13    w_n: f64,
14    x1: f64,
15    b_coef: f64,
16    a_coef: f64,
17    c_coef: f64,
18    f_coef: f64,
19    x_taylor: f64,
20    inverse_x0: f64,
21    taylor: [f64; TAYLOR_LENGTH],
22}
23
24impl HyperlogTransform {
25    pub fn new(t: f64, w: f64, m: f64, a: f64) -> Result<Self, FlowGateError> {
26        validate_hyperlog(t, w, m, a)?;
27
28        let w_n = w / (m + a);
29        let x2 = a / (m + a);
30        let x1 = x2 + w_n;
31        let x0 = x2 + 2.0 * w_n;
32        let b_coef = (m + a) * LN_10;
33
34        let e2bx0 = (b_coef * x0).exp();
35        let c_over_a = e2bx0 / w_n;
36        let f_over_a = (b_coef * x1).exp() + c_over_a * x1;
37
38        let a_coef = t / ((b_coef.exp() + c_over_a) - f_over_a);
39        let c_coef = c_over_a * a_coef;
40        let f_coef = f_over_a * a_coef;
41
42        let x_taylor = x1 + w_n / 4.0;
43        let mut taylor = [0.0_f64; TAYLOR_LENGTH];
44        let mut coef = a_coef * (b_coef * x1).exp();
45        for (i, entry) in taylor.iter_mut().enumerate() {
46            coef *= b_coef / ((i + 1) as f64);
47            *entry = coef;
48        }
49        taylor[0] += c_coef;
50
51        let inverse_x0 = if x0 < x_taylor {
52            taylor_series_at(&taylor, x1, x0)
53        } else {
54            (a_coef * (b_coef * x0).exp() + c_coef * x0) - f_coef
55        };
56
57        Ok(Self {
58            t,
59            w,
60            m,
61            a,
62            w_n,
63            x1,
64            b_coef,
65            a_coef,
66            c_coef,
67            f_coef,
68            x_taylor,
69            inverse_x0,
70            taylor,
71        })
72    }
73
74    #[inline]
75    fn taylor_series(&self, scale: f64) -> f64 {
76        taylor_series_at(&self.taylor, self.x1, scale)
77    }
78}
79
80impl Transform for HyperlogTransform {
81    fn apply(&self, value: f64) -> f64 {
82        if !value.is_finite() {
83            return f64::NAN;
84        }
85        if value == 0.0 {
86            return self.x1;
87        }
88
89        let negative = value < 0.0;
90        let signal = value.abs();
91
92        let mut x = if signal < self.inverse_x0 {
93            self.x1 + signal * self.w_n / self.inverse_x0
94        } else {
95            (signal / self.a_coef).ln() / self.b_coef
96        };
97        if !x.is_finite() {
98            return f64::NAN;
99        }
100
101        let tolerance = if x > 1.0 {
102            3.0 * x * f64::EPSILON
103        } else {
104            3.0 * f64::EPSILON
105        };
106
107        for _ in 0..40 {
108            let ae2bx = self.a_coef * (self.b_coef * x).exp();
109            let y = if x < self.x_taylor {
110                self.taylor_series(x) - signal
111            } else {
112                (ae2bx + self.c_coef * x) - (self.f_coef + signal)
113            };
114
115            let abe2bx = self.b_coef * ae2bx;
116            let dy = abe2bx + self.c_coef;
117            let ddy = self.b_coef * abe2bx;
118            let denom = dy * (1.0 - y * ddy / (2.0 * dy * dy));
119            if !denom.is_finite() || denom.abs() < f64::EPSILON {
120                return f64::NAN;
121            }
122            let delta = y / denom;
123            if !delta.is_finite() {
124                return f64::NAN;
125            }
126            x -= delta;
127            if delta.abs() < tolerance {
128                return if negative { 2.0 * self.x1 - x } else { x };
129            }
130        }
131        f64::NAN
132    }
133
134    fn invert(&self, scaled: f64) -> f64 {
135        if !scaled.is_finite() {
136            return f64::NAN;
137        }
138
139        let negative = scaled < self.x1;
140        let reflected = if negative {
141            2.0 * self.x1 - scaled
142        } else {
143            scaled
144        };
145
146        let inverse = if reflected < self.x_taylor {
147            self.taylor_series(reflected)
148        } else {
149            (self.a_coef * (self.b_coef * reflected).exp() + self.c_coef * reflected) - self.f_coef
150        };
151
152        if !inverse.is_finite() {
153            return f64::NAN;
154        }
155        if negative {
156            -inverse
157        } else {
158            inverse
159        }
160    }
161
162    fn transform_id(&self) -> &str {
163        "hyperlog"
164    }
165}
166
167#[inline]
168fn taylor_series_at(taylor: &[f64; TAYLOR_LENGTH], x1: f64, scale: f64) -> f64 {
169    let x = scale - x1;
170    let mut sum = taylor[TAYLOR_LENGTH - 1] * x;
171    for i in (0..(TAYLOR_LENGTH - 1)).rev() {
172        sum = (sum + taylor[i]) * x;
173    }
174    sum
175}
176
177fn validate_hyperlog(t: f64, w: f64, m: f64, a: f64) -> Result<(), FlowGateError> {
178    if !t.is_finite() || t <= 0.0 {
179        return Err(FlowGateError::InvalidTransformParam(
180            "Hyperlog parameter T must be finite and > 0".to_string(),
181        ));
182    }
183    if !m.is_finite() || m <= 0.0 {
184        return Err(FlowGateError::InvalidTransformParam(
185            "Hyperlog parameter M must be finite and > 0".to_string(),
186        ));
187    }
188    if !w.is_finite() || w <= 0.0 {
189        return Err(FlowGateError::InvalidTransformParam(
190            "Hyperlog parameter W must be finite and > 0".to_string(),
191        ));
192    }
193    if !a.is_finite() {
194        return Err(FlowGateError::InvalidTransformParam(
195            "Hyperlog parameter A must be finite".to_string(),
196        ));
197    }
198    if 2.0 * w > m {
199        return Err(FlowGateError::InvalidTransformParam(
200            "Hyperlog constraint violated: 2W <= M".to_string(),
201        ));
202    }
203    if -a > w || (a + w) > (m - w) {
204        return Err(FlowGateError::InvalidTransformParam(
205            "Hyperlog parameters violate A/W/M constraints".to_string(),
206        ));
207    }
208    Ok(())
209}