pairs.rs

1use crate::moire::{layer, Layer, Spec};2use mrlycore::errors::{value_error, Result};3use mrlynum::factor::lcm;45fn odd_blocks(span: u64, block: u64) -> u64 {6    block * (span / (2 * block)) + (span % (2 * block)).saturating_sub(block)7}89fn overlap(m: u64, n: u64) -> (u64, u64) {10    let grid = lcm(m as usize, n as usize) as u64;11    let (a, b) = (grid / m, grid / n);12    let both = (1..n)13        .step_by(2)14        .map(|k| odd_blocks((k + 1) * b, a) - odd_blocks(k * b, a))15        .sum();16    (both, grid)17}1819/// Returns the exact Pearson correlation of the flat carpet layers at two scales, area-weighted on their lcm grid.20///21/// A layer at scale n is lit where the row and the column of the n by n grid are not both odd.22/// The correlation is exactly zero when the two odd scales are coprime, and zero by convention below scale two, where a layer is constant.23///24/// ```25/// assert_eq!(mrlylab::moire::pairs::correlation(3, 5), 0.0);26/// assert!(mrlylab::moire::pairs::correlation(3, 9) > 0.0);27/// ```28pub fn correlation(m: usize, n: usize) -> f64 {29    let (m, n) = (m as u64, n as u64);30    let (hm, hn) = (m / 2, n / 2);31    if hm == 0 || hn == 0 {32        return 0.0;33    }34    let (both, grid) = overlap(m, n);35    let cross = (both * m * n) as i128;36    let solo = (hm * hn * grid) as i128;37    let gap = cross - solo;38    if gap == 0 {39        return 0.0;40    }41    let scale = (grid * m * n) as f64;42    let covariance = (gap as f64 / scale) * ((cross + solo) as f64 / scale);43    let variance = |half: u64, side: u64| {44        let share = (half * half) as f64 / (side * side) as f64;45        share * (1.0 - share)46    };47    covariance / (variance(hm, m) * variance(hn, n)).sqrt()48}4950/// The witness row of an odd scale: its correlation with every earlier odd scale from three, and the verdict the row gives.51#[derive(Clone, Debug, PartialEq)]52pub struct Witness {53    /// The scale on trial.54    pub scale: usize,55    /// The earlier odd scales, three up to the scale less two.56    pub scales: Vec<usize>,57    /// The exact correlation with each earlier scale.58    pub row: Vec<f64>,59    /// The largest correlation in the row, zero for an empty row.60    pub max: f64,61    /// The earlier scale carrying the largest correlation, zero when the row is clear.62    pub at: usize,63    /// Whether the row is exactly clear, which is the scale being prime.64    pub prime: bool,65}6667/// Puts an odd scale of three or more on trial against every earlier odd scale, or an error for another scale.68///69/// ```70/// let trial = mrlylab::moire::pairs::witness(9).unwrap();71/// assert_eq!((trial.scales, trial.at, trial.prime), (vec![3, 5, 7], 3, false));72/// ```73pub fn witness(scale: usize) -> Result<Witness> {74    if scale < 3 || scale.is_multiple_of(2) {75        return value_error("the stack has odd scales from three.");76    }77    let scales: Vec<usize> = (3..scale).step_by(2).collect();78    let row: Vec<f64> = scales.iter().map(|&m| correlation(m, scale)).collect();79    let (mut max, mut at) = (0.0, 0);80    for (&m, &r) in scales.iter().zip(&row) {81        if r > max {82            (max, at) = (r, m);83        }84    }85    let prime = row.iter().all(|&r| r == 0.0);86    Ok(Witness {87        scale,88        scales,89        row,90        max,91        at,92        prime,93    })94}9596/// Returns the Pearson correlation of two rendered carpet layers on their lcm grid, sampled rather than integrated.97pub fn sampled(m: usize, n: usize) -> f64 {98    let size = lcm(m, n);99    let mask = |number| {100        let params = Layer {101            size,102            ..Layer::new(Spec::new(7, 2, 2), number)103        };104        layer(&params).unwrap()105    };106    let (a, b) = (mask(m), mask(n));107    let mean = |v: &[bool]| v.iter().filter(|&&x| x).count() as f64 / v.len() as f64;108    let (ea, eb) = (mean(&a), mean(&b));109    let eab = a.iter().zip(&b).filter(|(&x, &y)| x && y).count() as f64 / a.len() as f64;110    (eab - ea * eb) / (ea * (1.0 - ea) * eb * (1.0 - eb)).sqrt()111}112113#[cfg(test)]114mod tests {115    use super::*;116117    fn brute(m: u64, n: u64) -> (u64, u64) {118        let grid = lcm(m as usize, n as usize) as u64;119        let (a, b) = (grid / m, grid / n);120        let both = (0..grid)121            .filter(|&j| !(j / a).is_multiple_of(2) && !(j / b).is_multiple_of(2))122            .count() as u64;123        (both, grid)124    }125126    #[test]127    fn the_closed_count_matches_the_lcm_grid() {128        for m in 2..30u64 {129            for n in 2..30u64 {130                assert_eq!(overlap(m, n), brute(m, n), "{m} {n}");131            }132        }133    }134}