mass.rs

1use crate::design::BASE;2use mrlycore::Tensor;34pub struct Mass {5    pub radii: Vec<f64>,6    pub cumulative: Vec<u64>,7}89impl Mass {10    pub fn at(&self, radius: f64) -> u64 {11        match self12            .radii13            .binary_search_by(|probe| probe.partial_cmp(&radius).expect("finite radii"))14        {15            Ok(index) => self.cumulative[index],16            Err(0) => 0,17            Err(index) => self.cumulative[index - 1],18        }19    }2021    pub fn total(&self) -> u64 {22        *self.cumulative.last().unwrap_or(&0)23    }24}2526pub fn shells(grid: &Tensor, digit: (usize, usize)) -> Mass {27    let side = grid.shape[0];28    let bytes = grid.bytes();29    let anchor = (side as i64 * digit.0 as i64, side as i64 * digit.1 as i64);30    let mut keys: Vec<u32> = Vec::new();31    for row in 0..side {32        for col in 0..side {33            if bytes[row * side + col] == 0 {34                continue;35            }36            let dr = 2 * row as i64 + 1 - anchor.0;37            let dc = 2 * col as i64 + 1 - anchor.1;38            keys.push((dr * dr + dc * dc) as u32);39        }40    }41    keys.sort_unstable();42    let mut radii = Vec::new();43    let mut cumulative = Vec::new();44    let mut running = 0u64;45    for (index, key) in keys.iter().enumerate() {46        running += 1;47        if index + 1 == keys.len() || keys[index + 1] != *key {48            radii.push((*key as f64).sqrt() / 2.0);49            cumulative.push(running);50        }51    }52    Mass { radii, cumulative }53}5455/// The nearest filled cell to the fixed point of the digit, as four times its squared distance, an exact integer.56pub fn nearest_cell(grid: &Tensor, digit: (usize, usize)) -> u64 {57    let side = grid.shape[0];58    let bytes = grid.bytes();59    let anchor = (side as i64 * digit.0 as i64, side as i64 * digit.1 as i64);60    let reach = |at: i64, low: i64| -> i64 {61        let (a, b) = (2 * low - at, at - 2 * (low + 1));62        a.max(b).max(0)63    };64    let mut best = u64::MAX;65    for row in 0..side {66        for col in 0..side {67            if bytes[row * side + col] == 0 {68                continue;69            }70            let dr = reach(anchor.0, row as i64);71            let dc = reach(anchor.1, col as i64);72            best = best.min((dr * dr + dc * dc) as u64);73        }74    }75    best76}7778pub fn horizon(code: u128, digit: (usize, usize), table: &[(usize, usize)]) -> f64 {79    let point = (digit.0 as f64 / 2.0, digit.1 as f64 / 2.0);80    let mut best = f64::INFINITY;81    for bit in 0..9 {82        if code >> bit & 1 == 0 {83            continue;84        }85        let other = table[bit];86        if other == digit {87            continue;88        }89        let gap = |at: f64, index: usize| {90            let (low, high) = (index as f64 / 3.0, (index as f64 + 1.0) / 3.0);91            (low - at).max(at - high).max(0.0)92        };93        let (dr, dc) = (gap(point.0, other.0), gap(point.1, other.1));94        best = best.min((dr * dr + dc * dc).sqrt());95    }96    (3.0 * best).min(1.0)97}9899pub struct Ripple {100    pub slope: f64,101    pub curve: Vec<f64>,102    pub swing: f64,103    pub drift: f64,104    pub periods: usize,105}106107fn sample(mass: &Mass, low: f64, periods: usize, bins: usize) -> Vec<(f64, f64)> {108    let step = (BASE as f64).ln() / bins as f64;109    (0..periods * bins)110        .filter_map(|index| {111            let radius = low * (index as f64 * step).exp();112            let count = mass.at(radius);113            if count == 0 {114                None115            } else {116                Some((radius.ln(), (count as f64).ln()))117            }118        })119        .collect()120}121122fn fit(points: &[(f64, f64)]) -> (f64, f64) {123    let n = points.len() as f64;124    let sx: f64 = points.iter().map(|p| p.0).sum();125    let sy: f64 = points.iter().map(|p| p.1).sum();126    let sxx: f64 = points.iter().map(|p| p.0 * p.0).sum();127    let sxy: f64 = points.iter().map(|p| p.0 * p.1).sum();128    let slope = (n * sxy - sx * sy) / (n * sxx - sx * sx);129    (slope, (sy - slope * sx) / n)130}131132fn fold(points: &[(f64, f64)], dimension: f64, bins: usize) -> Vec<f64> {133    let mut sums = vec![0.0; bins];134    let mut hits = vec![0.0; bins];135    for (lx, ly) in points {136        let residual = ly - dimension * lx;137        let phase = (lx / (BASE as f64).ln()).rem_euclid(1.0);138        let bin = ((phase * bins as f64) as usize).min(bins - 1);139        sums[bin] += residual;140        hits[bin] += 1.0;141    }142    let curve: Vec<f64> = sums143        .iter()144        .zip(&hits)145        .map(|(sum, hit)| if *hit == 0.0 { f64::NAN } else { sum / hit })146        .collect();147    let live: Vec<f64> = curve.iter().copied().filter(|v| v.is_finite()).collect();148    let mean = live.iter().sum::<f64>() / live.len().max(1) as f64;149    curve.iter().map(|v| v - mean).collect()150}151152/// The whole powers of the base that fit between the window ends, counted by repeated multiplication so no logarithm can round a boundary the wrong way.153pub fn periods(low: f64, high: f64) -> usize {154    let mut count = 0usize;155    let mut edge = low * BASE as f64;156    while edge <= high * (1.0 + 1e-12) {157        count += 1;158        edge *= BASE as f64;159    }160    count161}162163pub fn ripple(mass: &Mass, low: f64, high: f64, dimension: f64, bins: usize) -> Ripple {164    let periods = periods(low, high).max(1);165    let points = sample(mass, low, periods, bins);166    let (slope, _) = fit(&points);167    let curve = fold(&points, dimension, bins);168    let swing = spread(&curve);169    let half = periods / 2;170    let drift = if half == 0 {171        f64::NAN172    } else {173        let early = fold(&sample(mass, low, half, bins), dimension, bins);174        let late = fold(175            &sample(176                mass,177                low * (BASE as f64).powi(half as i32),178                periods - half,179                bins,180            ),181            dimension,182            bins,183        );184        early185            .iter()186            .zip(&late)187            .map(|(a, b)| (a - b).abs())188            .fold(189                0.0f64,190                |best, gap| if gap.is_finite() { best.max(gap) } else { best },191            )192    };193    Ripple {194        slope,195        curve,196        swing,197        drift,198        periods,199    }200}201202pub fn spread(curve: &[f64]) -> f64 {203    let live: Vec<f64> = curve.iter().copied().filter(|v| v.is_finite()).collect();204    let low = live.iter().copied().fold(f64::INFINITY, f64::min);205    let high = live.iter().copied().fold(f64::NEG_INFINITY, f64::max);206    high - low207}208209pub fn distance(left: &[f64], right: &[f64]) -> f64 {210    left.iter()211        .zip(right)212        .map(|(a, b)| (a - b).abs())213        .fold(214            0.0f64,215            |best, gap| if gap.is_finite() { best.max(gap) } else { best },216        )217}218219pub fn scaling_error(mass: &Mass, low: f64, high: f64, fill: f64) -> f64 {220    let mut worst: f64 = 0.0;221    let mut radius = low;222    while radius * (BASE as f64) <= high {223        let near = mass.at(radius) as f64;224        let far = mass.at(radius * BASE as f64) as f64;225        if near > 0.0 {226            worst = worst.max((far / (fill * near) - 1.0).abs());227        }228        radius *= 1.05;229    }230    worst231}232233#[cfg(test)]234mod tests {235    use super::periods;236237    #[test]238    fn the_window_counts_whole_periods_at_every_level() {239        for level in 4..=12u32 {240            let side = 3f64.powi(level as i32);241            assert_eq!(242                periods(27.0, side),243                level as usize - 3,244                "corner level {level}"245            );246            assert_eq!(247                periods(27.0, side / 2.0),248                level as usize - 4,249                "centre level {level}"250            );251        }252        assert_eq!(periods(27.0, 26.0), 0);253        assert_eq!(periods(27.0, 81.0), 1);254        assert_eq!(periods(27.0, 80.9), 0);255    }256}