spin.rs

1use mrlylab::moire::{stack, Combine, Field, Lattice, Spec};2use mrlymath::bang::corners_to_code;3use mrlynum::spin::{profile, reach};45const RINGS: [(usize, usize, f64); 3] = [(3, 5, 0.38), (5, 7, -0.33), (9, 13, 0.38)];67const SIZE: usize = 135;89const STEPS: usize = 400;1011// RING PROFILES1213pub fn the_coprime_law_survives_the_spin() -> Result<(), String> {14    let far = reach(SIZE);15    let disc: Vec<f64> = (0..STEPS)16        .map(|k| {17            let radius = k as f64 * far / (STEPS - 1) as f64;18            if radius <= SIZE as f64 / 2.0 {19                radius20            } else {21                0.022            }23        })24        .collect();25    let ones = vec![1.0; SIZE * SIZE];26    for (m, n, ring) in RINGS {27        let (a, b) = (layer(m)?, layer(n)?);28        let flat = pearson(&a.as_f64(), &b.as_f64(), &ones);29        if flat.abs() >= 0.03 {30            return Err(format!("the flat layers {m} and {n} correlate at {flat}"));31        }32        let spun = pearson(&profiled(&a), &profiled(&b), &disc);33        if (spun - ring).abs() >= 0.02 {34            return Err(format!("the ring profiles {m} and {n} correlate at {spun}"));35        }36    }37    Ok(())38}3940fn layer(scale: usize) -> Result<Field, String> {41    let spec = Spec::new(corners_to_code(&[vec![0, 0]], 2, 2), 2, 2);42    stack(spec, &[scale], Combine::Sum, 1, Lattice::Square, SIZE, &[])43        .map_err(|_| format!("the layer at {scale} does not build"))44}4546fn profiled(field: &Field) -> Vec<f64> {47    profile(&field.data, SIZE, STEPS)48        .iter()49        .map(|&value| value as f64)50        .collect()51}5253fn pearson(a: &[f64], b: &[f64], weight: &[f64]) -> f64 {54    let total: f64 = weight.iter().sum();55    let mean = |x: &[f64]| x.iter().zip(weight).map(|(v, w)| v * w).sum::<f64>() / total;56    let (ma, mb) = (mean(a), mean(b));57    let covariance: f64 = a58        .iter()59        .zip(b)60        .zip(weight)61        .map(|((x, y), w)| w * (x - ma) * (y - mb))62        .sum();63    let spread = |x: &[f64], m: f64| {64        x.iter()65            .zip(weight)66            .map(|(v, w)| w * (v - m) * (v - m))67            .sum::<f64>()68    };69    covariance / (spread(a, ma) * spread(b, mb)).sqrt()70}