unfolding.rs

1use crate::numerics::mean;2use faer::prelude::SolveLstsq;3use faer::Mat;45pub struct Spacings {6    pub values: Vec<f64>,7    pub negatives: usize,8}910fn spacings(unfolded: &[f64]) -> Spacings {11    let raw: Vec<f64> = unfolded.windows(2).map(|pair| pair[1] - pair[0]).collect();12    let negatives = raw.iter().filter(|s| **s < 0.0).count();13    let clamped: Vec<f64> = raw.iter().map(|s| s.max(0.0)).collect();14    let scale = mean(&clamped);15    Spacings {16        values: clamped.iter().map(|s| s / scale).collect(),17        negatives,18    }19}2021fn chebyshev(x: f64, degree: usize) -> Vec<f64> {22    let mut out = vec![1.0, x];23    for k in 2..=degree {24        out.push(2.0 * x * out[k - 1] - out[k - 2]);25    }26    out.truncate(degree + 1);27    out28}2930pub fn polynomial(values: &[f64], degree: usize) -> Spacings {31    let n = values.len();32    let (low, high) = (values[0], values[n - 1]);33    let scaled: Vec<f64> = values34        .iter()35        .map(|v| 2.0 * (v - low) / (high - low) - 1.0)36        .collect();37    let mut design = Mat::<f64>::zeros(n, degree + 1);38    let mut target = Mat::<f64>::zeros(n, 1);39    for (i, x) in scaled.iter().enumerate() {40        for (k, basis) in chebyshev(*x, degree).iter().enumerate() {41            design[(i, k)] = *basis;42        }43        target[(i, 0)] = i as f64 + 0.5;44    }45    let coefficients = design.qr().solve_lstsq(&target);46    let unfolded: Vec<f64> = scaled47        .iter()48        .map(|x| {49            chebyshev(*x, degree)50                .iter()51                .enumerate()52                .map(|(k, basis)| coefficients[(k, 0)] * basis)53                .sum()54        })55        .collect();56    spacings(&unfolded)57}5859pub fn window(values: &[f64], half: usize) -> Spacings {60    let n = values.len();61    let mut unfolded = vec![0.0; n];62    for i in 1..n {63        let lo = i.saturating_sub(half);64        let hi = (i + half).min(n - 1);65        let width = values[hi] - values[lo];66        let gap = values[i] - values[i - 1];67        let step = if gap > 0.0 {68            gap * (hi - lo) as f64 / width69        } else {70            0.071        };72        unfolded[i] = unfolded[i - 1] + step;73    }74    spacings(&unfolded)75}