pairs.rs

1use crate::core::error::{value_error, Result};2use crate::math::moire::{layer, Layer, Spec};3use crate::num::factor::lcm;4use serde::{Deserialize, Serialize};56fn odd_blocks(span: u64, block: u64) -> u64 {7    block * (span / (2 * block)) + (span % (2 * block)).saturating_sub(block)8}910fn overlap(m: u64, n: u64) -> (u64, u64) {11    let grid = lcm(m as usize, n as usize) as u64;12    let (a, b) = (grid / m, grid / n);13    let both = (1..n)14        .step_by(2)15        .map(|k| odd_blocks((k + 1) * b, a) - odd_blocks(k * b, a))16        .sum();17    (both, grid)18}1920/// Returns the exact Pearson correlation of the flat carpet layers at two scales, area-weighted on their lcm grid.21///22/// A layer at scale n is lit where the row and the column of the n by n grid are not both odd.23/// The correlation is exactly zero when the two odd scales are coprime, and zero by convention below scale two, where a layer is constant.24///25/// ```26/// assert_eq!(mrlyrs::math::moire::pairs::correlation(3, 5), 0.0);27/// assert!(mrlyrs::math::moire::pairs::correlation(3, 9) > 0.0);28/// ```29pub fn correlation(m: usize, n: usize) -> f64 {30    let (m, n) = (m as u64, n as u64);31    let (hm, hn) = (m / 2, n / 2);32    if hm == 0 || hn == 0 {33        return 0.0;34    }35    let (both, grid) = overlap(m, n);36    let cross = (both * m * n) as i128;37    let solo = (hm * hn * grid) as i128;38    let gap = cross - solo;39    if gap == 0 {40        return 0.0;41    }42    let scale = (grid * m * n) as f64;43    let covariance = (gap as f64 / scale) * ((cross + solo) as f64 / scale);44    let variance = |half: u64, side: u64| {45        let share = (half * half) as f64 / (side * side) as f64;46        share * (1.0 - share)47    };48    covariance / (variance(hm, m) * variance(hn, n)).sqrt()49}5051/// The witness row of an odd scale: its correlation with every earlier odd scale from three, and the verdict the row gives.52#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]53pub struct Witness {54    /// The scale on trial.55    pub scale: usize,56    /// The earlier odd scales, three up to the scale less two.57    pub scales: Vec<usize>,58    /// The exact correlation with each earlier scale.59    pub row: Vec<f64>,60    /// The largest correlation in the row, zero for an empty row.61    pub max: f64,62    /// The earlier scale carrying the largest correlation, zero when the row is clear.63    pub at: usize,64    /// Whether the row is exactly clear, which is the scale being prime.65    pub prime: bool,66}6768/// Puts an odd scale of three or more on trial against every earlier odd scale.69///70/// ```71/// let trial = mrlyrs::math::moire::pairs::witness(9).unwrap();72/// assert_eq!((trial.scales, trial.at, trial.prime), (vec![3, 5, 7], 3, false));73/// ```74///75/// # Errors76///77/// Errors for a scale that is not odd and three or more.78pub fn witness(scale: usize) -> Result<Witness> {79    if scale < 3 || scale.is_multiple_of(2) {80        return value_error("the stack has odd scales from three.");81    }82    let scales: Vec<usize> = (3..scale).step_by(2).collect();83    let row: Vec<f64> = scales.iter().map(|&m| correlation(m, scale)).collect();84    let (mut max, mut at) = (0.0, 0);85    for (&m, &r) in scales.iter().zip(&row) {86        if r > max {87            (max, at) = (r, m);88        }89    }90    let prime = row.iter().all(|&r| r == 0.0);91    Ok(Witness {92        scale,93        scales,94        row,95        max,96        at,97        prime,98    })99}100101/// Returns the Pearson correlation of two rendered carpet layers on their lcm grid, sampled rather than integrated.102///103/// # Errors104///105/// Errors when either scale is zero.106pub fn sampled(m: usize, n: usize) -> Result<f64> {107    if m == 0 || n == 0 {108        return value_error("a sampled pair needs two scales of at least one.");109    }110    let size = lcm(m, n);111    let mask = |number| {112        let params = Layer {113            size,114            ..Layer::new(Spec::new(7, 2, 2), number)115        };116        layer(&params)117    };118    let (a, b) = (mask(m)?, mask(n)?);119    let mean = |v: &[bool]| v.iter().filter(|&&x| x).count() as f64 / v.len() as f64;120    let (ea, eb) = (mean(&a), mean(&b));121    let eab = a.iter().zip(&b).filter(|(&x, &y)| x && y).count() as f64 / a.len() as f64;122    Ok((eab - ea * eb) / (ea * (1.0 - ea) * eb * (1.0 - eb)).sqrt())123}124125#[cfg(test)]126mod tests {127    use super::*;128129    fn brute(m: u64, n: u64) -> (u64, u64) {130        let grid = lcm(m as usize, n as usize) as u64;131        let (a, b) = (grid / m, grid / n);132        let both = (0..grid)133            .filter(|&j| !(j / a).is_multiple_of(2) && !(j / b).is_multiple_of(2))134            .count() as u64;135        (both, grid)136    }137138    #[test]139    fn the_closed_count_matches_the_lcm_grid() {140        for m in 2..30u64 {141            for n in 2..30u64 {142                assert_eq!(overlap(m, n), brute(m, n), "{m} {n}");143            }144        }145    }146147    #[test]148    fn refuses_a_zero_scale() {149        assert!(sampled(0, 3).is_err());150        assert!(sampled(3, 0).is_err());151    }152}