frame.rs

1use crate::lattice::{cell, column, row, Rule};2use crate::sums::{mean, odds};34#[derive(Clone, Copy)]5pub enum Line {6    X,7    Z,8    W,9    D,10    K,11}1213impl Line {14    pub fn name(self) -> &'static str {15        match self {16            Line::X => "X",17            Line::Z => "Z",18            Line::W => "X+Z",19            Line::D => "X-Z",20            Line::K => "2X+Z",21        }22    }2324    fn at(self, x: f64, z: f64) -> f64 {25        match self {26            Line::X => x,27            Line::Z => z,28            Line::W => x + z,29            Line::D => x - z,30            Line::K => 2.0 * x + z,31        }32    }33}3435pub struct Frame {36    pub side: usize,37    pub layers: usize,38    pub values: Vec<f64>,39    pub whole: Vec<bool>,40    pub points: Vec<f64>,41}4243pub fn points(side: usize) -> Vec<f64> {44    (0..side)45        .map(|step| (step as f64 + 0.5) / side as f64)46        .collect()47}4849pub fn stack(limit: usize, rule: &Rule, side: usize) -> Frame {50    let points = points(side);51    let mut values = vec![0.0f64; side * side];52    let mut whole = vec![true; side * side];53    let mut layers = 0;54    let mut line = vec![0.0f64; side];55    let mut reach = vec![true; side];56    for number in odds(limit) {57        let n = number as i64;58        let columns: Vec<i64> = points.iter().map(|p| column(*p, n)).collect();59        let rows: Vec<i64> = points.iter().map(|p| row(*p, n)).collect();60        let mut start = 0;61        while start < side {62            let x = columns[start];63            let mut end = start;64            while end < side && columns[end] == x {65                end += 1;66            }67            for (col, z) in rows.iter().enumerate() {68                let found = cell(rule, n, x, *z);69                line[col] = f64::from(found.unwrap_or(false));70                reach[col] = found.is_some();71            }72            let clipped = reach.iter().any(|ok| !ok);73            for slot in start..end {74                let base = slot * side;75                for (value, add) in values[base..base + side].iter_mut().zip(&line) {76                    *value += add;77                }78                if clipped {79                    for (keep, ok) in whole[base..base + side].iter_mut().zip(&reach) {80                        *keep &= ok;81                    }82                }83            }84            start = end;85        }86        layers += 1;87    }88    for value in values.iter_mut() {89        *value /= layers as f64;90    }91    Frame {92        side,93        layers,94        values,95        whole,96        points,97    }98}99100impl Frame {101    pub fn window(&self, line: Line, keep: impl Fn(f64) -> bool) -> f64 {102        let mut kept = Vec::new();103        for row in 0..self.side {104            for col in 0..self.side {105                let at = row * self.side + col;106                if self.whole[at] && keep(line.at(self.points[row], self.points[col])) {107                    kept.push(self.values[at]);108                }109            }110        }111        mean(&kept)112    }113114    pub fn hexmean(&self) -> f64 {115        let kept: Vec<f64> = self116            .values117            .iter()118            .zip(&self.whole)119            .filter(|(_, keep)| **keep)120            .map(|(value, _)| *value)121            .collect();122        mean(&kept)123    }124125    pub fn band(&self, line: Line, position: f64, width: f64) -> [f64; 3] {126        let left = self.window(line, |at| at >= position - width && at < position);127        let right = self.window(line, |at| at > position && at <= position + width);128        let centre = self.window(line, |at| (at - position).abs() <= width);129        [left, right, centre]130    }131}132133pub fn plateau(limit: usize, rule: &Rule, side: usize) -> (f64, f64) {134    let frame = stack(limit, rule, side);135    let background = frame.hexmean();136    let mut sums = [vec![0.0f64; side], vec![0.0f64; side]];137    let mut whole = [vec![true; side], vec![true; side]];138    for number in odds(limit) {139        let n = number as i64;140        for (slot, x) in [n - 1, n].into_iter().enumerate() {141            for (at, point) in frame.points.iter().enumerate() {142                match cell(rule, n, x, row(*point, n)) {143                    Some(fill) => sums[slot][at] += f64::from(fill),144                    None => whole[slot][at] = false,145                }146            }147        }148    }149    let strip = |slot: usize| {150        let kept: Vec<f64> = sums[slot]151            .iter()152            .zip(&whole[slot])153            .filter(|(_, keep)| **keep)154            .map(|(value, _)| value / frame.layers as f64)155            .collect();156        mean(&kept) - background157    };158    (strip(0), strip(1))159}160161pub fn quarter(rule: &Rule) {162    println!("carpet ideal frame at N = 55, 3601 samples per axis, one-sided bands of width 0.018");163    let frame = stack(55, rule, 3601);164    let background = frame.hexmean();165    println!("  hexagon mean {background:.6}");166    let spec = [167        (Line::X, 0.25),168        (Line::X, 0.75),169        (Line::X, 1.0 / 3.0),170        (Line::X, 0.2),171        (Line::X, 0.5),172        (Line::Z, 0.25),173        (Line::D, 0.25),174        (Line::W, 1.25),175    ];176    for (line, position) in spec {177        let [left, right, _] = frame.band(line, position, 0.018);178        println!(179            "  {} = {position:.4}: left {:+.6}  right {:+.6}",180            line.name(),181            left - background,182            right - background183        );184    }185    println!("quarter line X = 1/4, one-sided strips straight at the line, 3601 samples");186    for limit in [151usize, 301, 601, 1201] {187        let (left, right) = plateau(limit, rule, 3601);188        println!(189            "  N = {limit:4}: left {left:+.5}  right {right:+.5}  mean magnitude {:.5}  (limit 1/8)",190            (left.abs() + right.abs()) / 2.0191        );192    }193}194195pub fn void_star(rule: &Rule) {196    println!(197        "void ideal frame at N = 55, 3601 samples per axis, centred bands of half-width 0.004"198    );199    let frame = stack(55, rule, 3601);200    let background = frame.hexmean();201    println!("  hexagon mean {background:.6}  (limit 1/4)");202    for (line, position) in [203        (Line::X, 0.5),204        (Line::Z, 0.5),205        (Line::W, 1.0),206        (Line::D, 0.0),207        (Line::K, 1.5),208    ] {209        let [_, _, centre] = frame.band(line, position, 0.004);210        println!(211            "  {} = {position:.2}: ink {centre:.6}  ratio to mean {:.3}",212            line.name(),213            centre / background214        );215    }216}217218pub fn fading(rule: &Rule) {219    println!("carpet main diagonal X = Z, band half-width 0.01, 2801 samples per axis");220    let mut scaled = Vec::new();221    for count in [28usize, 100, 200, 400] {222        let frame = stack(2 * count - 1, rule, 2801);223        let excess = frame.window(Line::D, |at| at.abs() <= 0.01) - frame.hexmean();224        println!(225            "  L = {count:3}: excess {excess:+.6}  excess * L {:+.4}",226            excess * count as f64227        );228        scaled.push(excess * count as f64);229    }230    println!(231        "  slope of excess * L against ln L: {:+.4} (100 to 200)  {:+.4} (200 to 400)",232        (scaled[2] - scaled[1]) / 2f64.ln(),233        (scaled[3] - scaled[2]) / 2f64.ln()234    );235}