flow_gate_core/transform/
logicle_impl.rs

1// Public API for downstream consumers (e.g. visualization libraries
2// that need axis tick generation). Items are not used by the core
3// library itself but are intentionally exposed.
4#![allow(dead_code)]
5
6use rayon::prelude::*;
7
8pub const LOGICLE_IMPL_VERSION: &str = "logicle_v1";
9
10const TAYLOR_LENGTH: usize = 16;
11const DEFAULT_FORWARD_LUT_BINS: usize = 16_384;
12const DEFAULT_INVERSE_LUT_BINS: usize = 8_192;
13const MAX_FORWARD_LUT_BINS: usize = 32_768;
14const LUT_ERROR_THRESHOLD: f64 = 1e-4;
15const MIN_SPAN: f64 = 1e-12;
16
17#[derive(Debug, Clone, Copy, PartialEq)]
18pub struct LogicleParams {
19    pub t: f64,
20    pub w: f64,
21    pub m: f64,
22    pub a: f64,
23}
24
25impl Default for LogicleParams {
26    fn default() -> Self {
27        Self {
28            t: 262_144.0,
29            w: 0.5,
30            m: 4.5,
31            a: 0.0,
32        }
33    }
34}
35
36impl LogicleParams {
37    pub fn validate(self) -> Result<Self, String> {
38        if !self.t.is_finite() || self.t <= 0.0 {
39            return Err("Logicle parameter T must be finite and > 0".to_string());
40        }
41        if !self.m.is_finite() || self.m <= 0.0 {
42            return Err("Logicle parameter M must be finite and > 0".to_string());
43        }
44        if !self.w.is_finite() || self.w < 0.0 {
45            return Err("Logicle parameter W must be finite and >= 0".to_string());
46        }
47        if !self.a.is_finite() {
48            return Err("Logicle parameter A must be finite".to_string());
49        }
50        if 2.0 * self.w > self.m {
51            return Err("Logicle parameter W is too large for M".to_string());
52        }
53        if -self.a > self.w || (self.a + self.w) > (self.m - self.w) {
54            return Err("Logicle parameters violate A/W/M constraints".to_string());
55        }
56        Ok(self)
57    }
58}
59
60#[derive(Debug, Clone, PartialEq, Eq)]
61pub enum AxisTickKind {
62    Major,
63    Minor,
64    Zero,
65}
66
67#[derive(Debug, Clone)]
68pub struct AxisTick {
69    pub signal_value: f64,
70    pub transformed_position: f64,
71    pub kind: AxisTickKind,
72}
73
74#[derive(Debug, Clone)]
75struct LogicleContext {
76    params: LogicleParams,
77    // Internal coefficients for the modified biexponential:
78    // B(y) = a*exp(b*y) - c*exp(-d*y) - f
79    a_coef: f64,
80    b_coef: f64,
81    c_coef: f64,
82    d_coef: f64,
83    f_coef: f64,
84    x0: f64,
85    x1: f64,
86    x_taylor: f64,
87    taylor: [f64; TAYLOR_LENGTH],
88}
89
90impl LogicleContext {
91    fn new(params: LogicleParams) -> Result<Self, String> {
92        let params = params.validate()?;
93
94        // Parks/Moore parameterization in Gating-ML / flowCore implementations.
95        let w = params.w / (params.m + params.a);
96        let x2 = params.a / (params.m + params.a);
97        let x1 = x2 + w;
98        let x0 = x2 + 2.0 * w;
99
100        let b_coef = (params.m + params.a) * std::f64::consts::LN_10;
101        let d_coef = solve_d(b_coef, w)?;
102        let c_over_a = (x0 * (b_coef + d_coef)).exp();
103        let m_over_a = (b_coef * x1).exp() - c_over_a / (d_coef * x1).exp();
104        let a_coef = params.t / (((b_coef).exp() - m_over_a) - c_over_a / (d_coef).exp());
105        let c_coef = c_over_a * a_coef;
106        let f_coef = -m_over_a * a_coef;
107
108        let x_taylor = x1 + w / 4.0;
109        let mut taylor = [0.0_f64; TAYLOR_LENGTH];
110        let mut pos_coef = a_coef * (b_coef * x1).exp();
111        let mut neg_coef = -c_coef / (d_coef * x1).exp();
112        for (i, slot) in taylor.iter_mut().enumerate() {
113            pos_coef *= b_coef / (i as f64 + 1.0);
114            neg_coef *= -d_coef / (i as f64 + 1.0);
115            *slot = pos_coef + neg_coef;
116        }
117        // By construction this should be exactly zero in the ideal math.
118        taylor[1] = 0.0;
119
120        Ok(Self {
121            params,
122            a_coef,
123            b_coef,
124            c_coef,
125            d_coef,
126            f_coef,
127            x0,
128            x1,
129            x_taylor,
130            taylor,
131        })
132    }
133
134    fn params(&self) -> LogicleParams {
135        self.params
136    }
137
138    fn series_biexponential(&self, scale: f64) -> f64 {
139        // Taylor expansion around x1 to stabilize near zero.
140        let x = scale - self.x1;
141        let mut sum = self.taylor[TAYLOR_LENGTH - 1] * x;
142        for i in (2..=(TAYLOR_LENGTH - 2)).rev() {
143            sum = (sum + self.taylor[i]) * x;
144        }
145        (sum * x + self.taylor[0]) * x
146    }
147
148    fn inverse_positive(&self, scale: f64) -> f64 {
149        if scale < self.x_taylor {
150            self.series_biexponential(scale)
151        } else {
152            (self.a_coef * (self.b_coef * scale).exp() + self.f_coef)
153                - self.c_coef / (self.d_coef * scale).exp()
154        }
155    }
156
157    fn eval_forward_function(&self, scale: f64, signal: f64) -> f64 {
158        self.inverse_positive(scale) - signal
159    }
160
161    fn forward_positive_halley(&self, signal: f64) -> Option<f64> {
162        let mut x = if signal < self.f_coef {
163            self.x1 + signal / self.taylor[0]
164        } else {
165            (signal / self.a_coef).ln() / self.b_coef
166        };
167
168        let mut tolerance = 3.0 * f64::EPSILON;
169        if x > 1.0 {
170            tolerance = 3.0 * x * f64::EPSILON;
171        }
172
173        for _ in 0..16 {
174            let ae2bx = self.a_coef * (self.b_coef * x).exp();
175            let ce2mdx = self.c_coef / (self.d_coef * x).exp();
176            let y = if x < self.x_taylor {
177                self.series_biexponential(x) - signal
178            } else {
179                (ae2bx + self.f_coef) - (ce2mdx + signal)
180            };
181            let dy = self.b_coef * ae2bx + self.d_coef * ce2mdx;
182            if !dy.is_finite() || dy.abs() < f64::EPSILON {
183                return None;
184            }
185            let ddy = self.b_coef * self.b_coef * ae2bx - self.d_coef * self.d_coef * ce2mdx;
186            let denom = dy * (1.0 - y * ddy / (2.0 * dy * dy));
187            if !denom.is_finite() || denom.abs() < f64::EPSILON {
188                return None;
189            }
190            let delta = y / denom;
191            if !delta.is_finite() {
192                return None;
193            }
194            x -= delta;
195            if delta.abs() < tolerance {
196                return Some(x);
197            }
198        }
199        None
200    }
201
202    fn forward_positive_bisection(&self, signal: f64) -> Option<f64> {
203        // Monotonic root solve fallback for robustness.
204        let mut lo = self.x1;
205        let mut hi = (self.x1 + 1.0).max(1.0);
206
207        let mut f_lo = self.eval_forward_function(lo, signal);
208        if !f_lo.is_finite() {
209            return None;
210        }
211        if f_lo >= 0.0 {
212            return Some(lo);
213        }
214
215        let mut f_hi = self.eval_forward_function(hi, signal);
216        let mut guard = 0;
217        while (!f_hi.is_finite() || f_hi < 0.0) && guard < 64 {
218            let span = (hi - lo).max(0.5);
219            hi += span * 1.8;
220            f_hi = self.eval_forward_function(hi, signal);
221            guard += 1;
222        }
223        if !f_hi.is_finite() || f_hi < 0.0 {
224            return None;
225        }
226
227        for _ in 0..80 {
228            let mid = 0.5 * (lo + hi);
229            let f_mid = self.eval_forward_function(mid, signal);
230            if !f_mid.is_finite() {
231                return None;
232            }
233            if f_mid.abs() < 1e-12 || (hi - lo).abs() < 1e-12 {
234                return Some(mid);
235            }
236            if f_mid < 0.0 {
237                lo = mid;
238                f_lo = f_mid;
239            } else {
240                hi = mid;
241                f_hi = f_mid;
242            }
243            if (f_hi - f_lo).abs() < 1e-20 {
244                break;
245            }
246        }
247        Some(0.5 * (lo + hi))
248    }
249
250    fn forward_positive(&self, signal: f64) -> f64 {
251        if signal == 0.0 {
252            return self.x1;
253        }
254        if let Some(value) = self.forward_positive_halley(signal) {
255            return value;
256        }
257        self.forward_positive_bisection(signal).unwrap_or(self.x1)
258    }
259
260    fn forward(&self, value: f64) -> f64 {
261        if !value.is_finite() {
262            return f64::NAN;
263        }
264        if value == 0.0 {
265            return self.x1;
266        }
267        if value < 0.0 {
268            let pos = self.forward_positive(-value);
269            return 2.0 * self.x1 - pos;
270        }
271        self.forward_positive(value)
272    }
273
274    fn inverse(&self, value: f64) -> f64 {
275        if !value.is_finite() {
276            return f64::NAN;
277        }
278        let negative = value < self.x1;
279        let reflected = if negative {
280            2.0 * self.x1 - value
281        } else {
282            value
283        };
284        let inverse = self.inverse_positive(reflected);
285        if negative {
286            -inverse
287        } else {
288            inverse
289        }
290    }
291
292    fn linear_half_signal(&self) -> f64 {
293        self.inverse(self.x0).abs()
294    }
295}
296
297fn solve_d(b: f64, w: f64) -> Result<f64, String> {
298    // Root solve for d from Moore/Parks biexponential constraints.
299    if w == 0.0 {
300        return Ok(b);
301    }
302    let tolerance = 2.0 * b * f64::EPSILON;
303    let mut d_lo = 0.0_f64;
304    let mut d_hi = b;
305    let mut d = 0.5 * (d_lo + d_hi);
306    let mut last_delta = d_hi - d_lo;
307    let f_b = -2.0 * b.ln() + w * b;
308    let mut f = 2.0 * d.ln() + w * d + f_b;
309    let mut last_f = f64::NAN;
310
311    for _ in 0..40 {
312        let df = 2.0 / d + w;
313        let newton_unsafe = ((d - d_hi) * df - f) * ((d - d_lo) * df - f) >= 0.0
314            || (1.9 * f).abs() > (last_delta * df).abs();
315
316        let delta = if newton_unsafe {
317            let delta = 0.5 * (d_hi - d_lo);
318            d = d_lo + delta;
319            delta
320        } else {
321            let delta = f / df;
322            d -= delta;
323            delta
324        };
325
326        if !delta.is_finite() {
327            return Err("Logicle d solve failed: non-finite update".to_string());
328        }
329        if delta.abs() < tolerance {
330            return Ok(d);
331        }
332        last_delta = delta;
333
334        f = 2.0 * d.ln() + w * d + f_b;
335        if f == 0.0 || f == last_f {
336            return Ok(d);
337        }
338        last_f = f;
339
340        if f < 0.0 {
341            d_lo = d;
342        } else {
343            d_hi = d;
344        }
345    }
346
347    Err("Logicle d solve failed to converge".to_string())
348}
349
350#[derive(Debug, Clone)]
351pub struct LogicleLut {
352    context: LogicleContext,
353    signal_min: f64,
354    signal_max: f64,
355    signal_step: f64,
356    forward_lut: Vec<f32>,
357    domain_min: f64,
358    domain_max: f64,
359    domain_step: f64,
360    inverse_lut: Vec<f32>,
361}
362
363impl LogicleLut {
364    pub fn build(
365        params: LogicleParams,
366        signal_min: f64,
367        signal_max: f64,
368        forward_bins: usize,
369        inverse_bins: usize,
370    ) -> Result<Self, String> {
371        let context = LogicleContext::new(params)?;
372        let lo = signal_min.min(signal_max);
373        let mut hi = signal_min.max(signal_max);
374        if (hi - lo).abs() < MIN_SPAN {
375            hi = lo + MIN_SPAN;
376        }
377
378        let forward_bins = forward_bins.max(2);
379        let signal_step = (hi - lo) / (forward_bins as f64 - 1.0);
380        let mut forward_lut = vec![0.0_f32; forward_bins];
381        for (i, slot) in forward_lut.iter_mut().enumerate() {
382            let signal = lo + signal_step * i as f64;
383            *slot = context.forward(signal) as f32;
384        }
385
386        let mut domain_min = forward_lut.first().copied().unwrap_or(0.0) as f64;
387        let mut domain_max = forward_lut.last().copied().unwrap_or(1.0) as f64;
388        if !domain_min.is_finite()
389            || !domain_max.is_finite()
390            || (domain_max - domain_min).abs() < MIN_SPAN
391        {
392            domain_min = context.forward(lo);
393            domain_max = context.forward(hi);
394            if (domain_max - domain_min).abs() < MIN_SPAN {
395                domain_max = domain_min + MIN_SPAN;
396            }
397        }
398
399        let inverse_bins = inverse_bins.max(2);
400        let domain_step = (domain_max - domain_min) / (inverse_bins as f64 - 1.0);
401        let mut inverse_lut = vec![0.0_f32; inverse_bins];
402        for (i, slot) in inverse_lut.iter_mut().enumerate() {
403            let domain = domain_min + domain_step * i as f64;
404            *slot = context.inverse(domain) as f32;
405        }
406
407        Ok(Self {
408            context,
409            signal_min: lo,
410            signal_max: hi,
411            signal_step,
412            forward_lut,
413            domain_min,
414            domain_max,
415            domain_step,
416            inverse_lut,
417        })
418    }
419
420    pub fn build_adaptive(
421        params: LogicleParams,
422        signal_min: f64,
423        signal_max: f64,
424    ) -> Result<Self, String> {
425        let mut lut = Self::build(
426            params,
427            signal_min,
428            signal_max,
429            DEFAULT_FORWARD_LUT_BINS,
430            DEFAULT_INVERSE_LUT_BINS,
431        )?;
432        if lut.max_forward_abs_error_sampled(256) > LUT_ERROR_THRESHOLD {
433            lut = Self::build(
434                params,
435                signal_min,
436                signal_max,
437                MAX_FORWARD_LUT_BINS,
438                DEFAULT_INVERSE_LUT_BINS,
439            )?;
440        }
441        Ok(lut)
442    }
443
444    pub fn params(&self) -> LogicleParams {
445        self.context.params()
446    }
447
448    pub fn forward(&self, signal: f64) -> f64 {
449        if !signal.is_finite() {
450            return f64::NAN;
451        }
452        if signal < self.signal_min || signal > self.signal_max {
453            return self.context.forward(signal);
454        }
455        let idx_f = (signal - self.signal_min) / self.signal_step;
456        interpolate_lut(idx_f, &self.forward_lut)
457    }
458
459    pub fn inverse(&self, domain: f64) -> f64 {
460        if !domain.is_finite() {
461            return f64::NAN;
462        }
463        if domain < self.domain_min || domain > self.domain_max {
464            return self.context.inverse(domain);
465        }
466        let idx_f = (domain - self.domain_min) / self.domain_step;
467        interpolate_lut(idx_f, &self.inverse_lut)
468    }
469
470    pub fn max_forward_abs_error_sampled(&self, samples: usize) -> f64 {
471        let sample_count = samples.max(16);
472        let span = (self.signal_max - self.signal_min).max(MIN_SPAN);
473        let mut max_err = 0.0_f64;
474        for i in 0..sample_count {
475            let signal = self.signal_min + span * (i as f64 / (sample_count - 1) as f64);
476            let direct = self.context.forward(signal);
477            let approx = self.forward(signal);
478            let err = (direct - approx).abs();
479            if err > max_err {
480                max_err = err;
481            }
482        }
483        max_err
484    }
485}
486
487fn interpolate_lut(index: f64, lut: &[f32]) -> f64 {
488    if lut.is_empty() || !index.is_finite() {
489        return f64::NAN;
490    }
491    if lut.len() == 1 {
492        return lut[0] as f64;
493    }
494    let low = index.floor().clamp(0.0, (lut.len() - 1) as f64) as usize;
495    let high = (low + 1).min(lut.len() - 1);
496    let frac = (index - low as f64).clamp(0.0, 1.0);
497    let a = lut[low] as f64;
498    let b = lut[high] as f64;
499    a + (b - a) * frac
500}
501
502pub fn logicle_forward(value: f64, params: LogicleParams) -> f64 {
503    match LogicleContext::new(params) {
504        Ok(ctx) => ctx.forward(value),
505        Err(_) => f64::NAN,
506    }
507}
508
509pub fn logicle_inverse(value: f64, params: LogicleParams) -> f64 {
510    match LogicleContext::new(params) {
511        Ok(ctx) => ctx.inverse(value),
512        Err(_) => f64::NAN,
513    }
514}
515
516pub fn apply_transform_logicle(values: &[f32], params: LogicleParams, out: &mut [f32]) {
517    if values.len() != out.len() {
518        return;
519    }
520    if values.is_empty() {
521        return;
522    }
523
524    let mut min_signal = f64::INFINITY;
525    let mut max_signal = f64::NEG_INFINITY;
526    for value in values.iter().copied() {
527        let v = value as f64;
528        if !v.is_finite() {
529            continue;
530        }
531        min_signal = min_signal.min(v);
532        max_signal = max_signal.max(v);
533    }
534    if !min_signal.is_finite() || !max_signal.is_finite() {
535        out.fill(f32::NAN);
536        return;
537    }
538
539    let lut = match LogicleLut::build_adaptive(params, min_signal, max_signal) {
540        Ok(lut) => lut,
541        Err(_) => {
542            out.fill(f32::NAN);
543            return;
544        }
545    };
546
547    out.par_iter_mut()
548        .zip(values.par_iter())
549        .for_each(|(dst, src)| {
550            let value = *src as f64;
551            *dst = if value.is_finite() {
552                lut.forward(value) as f32
553            } else {
554                f32::NAN
555            };
556        });
557}
558
559fn nice_linear_step(span: f64, target_ticks: usize) -> f64 {
560    let target = target_ticks.max(2) as f64;
561    let rough = (span / target).abs().max(MIN_SPAN);
562    let mag = 10.0_f64.powf(rough.log10().floor());
563    let residual = rough / mag;
564    let nice = if residual <= 1.0 {
565        1.0
566    } else if residual <= 2.0 {
567        2.0
568    } else if residual <= 5.0 {
569        5.0
570    } else {
571        10.0
572    };
573    nice * mag
574}
575
576fn add_log_decade_ticks(
577    min_signal: f64,
578    max_signal: f64,
579    sign: f64,
580    out: &mut Vec<(f64, AxisTickKind)>,
581) {
582    if min_signal <= 0.0 || max_signal <= 0.0 || max_signal < min_signal {
583        return;
584    }
585    let start_pow = min_signal.log10().floor() as i32;
586    let end_pow = max_signal.log10().ceil() as i32;
587    let range_lo = if sign >= 0.0 { min_signal } else { -max_signal };
588    let range_hi = if sign >= 0.0 { max_signal } else { -min_signal };
589    for p in start_pow..=end_pow {
590        let base = 10.0_f64.powi(p);
591        for (mult, kind) in [
592            (1.0_f64, AxisTickKind::Major),
593            (2.0_f64, AxisTickKind::Minor),
594            (5.0_f64, AxisTickKind::Minor),
595        ] {
596            let candidate = sign * base * mult;
597            if candidate >= range_lo && candidate <= range_hi {
598                out.push((candidate, kind));
599            }
600        }
601    }
602}
603
604pub fn logicle_axis_ticks(
605    params: LogicleParams,
606    domain_min: f64,
607    domain_max: f64,
608) -> Result<Vec<AxisTick>, String> {
609    let ctx = LogicleContext::new(params)?;
610    let d_lo = domain_min.min(domain_max);
611    let d_hi = domain_min.max(domain_max);
612    if !d_lo.is_finite() || !d_hi.is_finite() || (d_hi - d_lo).abs() < MIN_SPAN {
613        return Ok(Vec::new());
614    }
615
616    let signal_a = ctx.inverse(d_lo);
617    let signal_b = ctx.inverse(d_hi);
618    if !signal_a.is_finite() || !signal_b.is_finite() {
619        return Ok(Vec::new());
620    }
621    let s_lo = signal_a.min(signal_b);
622    let s_hi = signal_a.max(signal_b);
623
624    let mut candidates: Vec<(f64, AxisTickKind)> = Vec::new();
625
626    if s_lo <= 0.0 && s_hi >= 0.0 {
627        candidates.push((0.0, AxisTickKind::Zero));
628    }
629    if s_hi > 0.0 {
630        add_log_decade_ticks(s_lo.max(MIN_SPAN), s_hi, 1.0, &mut candidates);
631    }
632    if s_lo < 0.0 {
633        add_log_decade_ticks((-s_hi).max(MIN_SPAN), -s_lo, -1.0, &mut candidates);
634    }
635
636    // Ensure linear neighborhood around zero has readable ticks.
637    let linear_half = ctx.linear_half_signal().max(1.0);
638    let linear_extent = linear_half.min(s_hi.abs().max(s_lo.abs()));
639    let linear_step = nice_linear_step(2.0 * linear_extent, 8);
640    for i in -4..=4 {
641        let signal = i as f64 * linear_step;
642        if signal >= s_lo && signal <= s_hi {
643            let kind = if signal == 0.0 {
644                AxisTickKind::Zero
645            } else {
646                AxisTickKind::Minor
647            };
648            candidates.push((signal, kind));
649        }
650    }
651
652    let mut ticks = Vec::new();
653    for (signal, kind) in candidates {
654        let transformed = ctx.forward(signal);
655        if !transformed.is_finite() {
656            continue;
657        }
658        if transformed < d_lo - 1e-9 || transformed > d_hi + 1e-9 {
659            continue;
660        }
661        ticks.push(AxisTick {
662            signal_value: signal,
663            transformed_position: transformed,
664            kind,
665        });
666    }
667
668    ticks.sort_by(|a, b| {
669        let by_pos = a
670            .transformed_position
671            .partial_cmp(&b.transformed_position)
672            .unwrap_or(std::cmp::Ordering::Equal);
673        if by_pos != std::cmp::Ordering::Equal {
674            return by_pos;
675        }
676        axis_kind_priority(&a.kind).cmp(&axis_kind_priority(&b.kind))
677    });
678    ticks.dedup_by(|a, b| {
679        (a.transformed_position - b.transformed_position).abs() < 1e-7
680            || (a.signal_value - b.signal_value).abs() < 1e-7
681    });
682
683    Ok(ticks)
684}
685
686fn axis_kind_priority(kind: &AxisTickKind) -> i32 {
687    match kind {
688        AxisTickKind::Zero => 0,
689        AxisTickKind::Major => 1,
690        AxisTickKind::Minor => 2,
691    }
692}
693
694#[cfg(test)]
695mod tests {
696    use super::*;
697
698    #[test]
699    fn round_trip_forward_inverse_is_stable() {
700        let p = LogicleParams::default();
701        let values = [
702            -10_000.0, -2_000.0, -250.0, -5.0, 0.0, 5.0, 250.0, 5_000.0, 50_000.0, 250_000.0,
703        ];
704        for x in values {
705            let y = logicle_forward(x, p);
706            let x2 = logicle_inverse(y, p);
707            let abs_err = (x - x2).abs();
708            let rel = abs_err / x.abs().max(1.0);
709            assert!(rel < 1e-6, "round-trip error too high at {x}: got {x2}");
710        }
711    }
712
713    #[test]
714    fn forward_is_monotonic() {
715        let ctx = LogicleContext::new(LogicleParams::default()).expect("context");
716        let mut last = f64::NEG_INFINITY;
717        for i in -10_000..=10_000 {
718            let signal = i as f64 * 30.0;
719            let y = ctx.forward(signal);
720            assert!(y >= last, "non-monotonic at {signal}: {y} < {last}");
721            last = y;
722        }
723    }
724
725    #[test]
726    fn lut_error_is_small() {
727        let params = LogicleParams::default();
728        let lut = LogicleLut::build_adaptive(params, -25_000.0, 262_144.0).expect("lut");
729        assert!(lut.max_forward_abs_error_sampled(512) < 1e-3);
730    }
731}
flow_gate_core/transform/logicle_impl.rs

flow_gate_core/transform/
logicle_impl.rs
