readings.rs

1use mrlyrs::core::error::Result;2use mrlyrs::core::tensor::Tensor;3use mrlyrs::life::{churn, entropy};4use mrlyrs::math::cell::graph::{core_graph, tunnel_graph};5use mrlyrs::math::graph::census::components;6use mrlyrs::math::two::{census, Cell2d};7use mrlyrs::num::fft::{peak_ring, radial_profile};8use std::f64::consts::PI;910/// The dihedral subgroup fixing a frame, by order and name.11#[derive(Clone, Copy, Debug, PartialEq, Eq)]12pub struct Symmetry {13    /// The number of square symmetries fixing the frame, 1 to 8.14    pub order: usize,15    /// The subgroup's name: d4, c4, d2, d2d, c2, m, md or 1.16    pub name: &'static str,17}1819/// The readings of one frame.20#[derive(Clone, Debug, PartialEq)]21pub struct Reading {22    /// The frame's side.23    pub side: usize,24    /// The number of filled cells.25    pub fill: usize,26    /// The 4-adjacent components of the filled cells on the plane.27    pub components: usize,28    /// The 4-adjacent components of the empty cells on the plane, the background included.29    pub holes: usize,30    /// The Euler characteristic of the filled cells.31    pub euler: i64,32    /// The least-squares slope of the dyadic box counts against the scale.33    pub box_slope: f64,34    /// The dihedral subgroup fixing the frame.35    pub symmetry: Symmetry,36    /// The ring of the strongest non-zero frequency read over the full Fourier square.37    pub ring: usize,38    /// The share of the non-zero power sitting on that ring.39    pub share: f64,40    /// The peak ring the ring cut reads, the corners dropped.41    pub ring_cut: usize,42    /// The frame's binary entropy in millibits.43    pub entropy: i64,44    /// The fraction of cells changed since the previous frame, zero without one.45    pub churn: f64,46    /// Every proper divisor of the side the block test cuts at.47    pub cuts: Vec<usize>,48    /// The smallest cut whose two factors pack into codes, as (d, outer, inner).49    pub factors: Option<(usize, u128, u128)>,50}5152/// Reads one frame, the churn taken against the previous frame when given.53pub fn read(frame: &Cell2d, previous: Option<&Cell2d>) -> Result<Reading> {54    let grid = frame.types();55    let side = grid.shape[0];56    let fill = grid.sum() as usize;57    let pieces = if fill == 0 {58        059    } else {60        components(&core_graph(frame)?)?61    };62    let holes = if fill == side * side {63        064    } else {65        components(&tunnel_graph(frame)?)?66    };67    let euler = census::euler(frame)?;68    let (ring, share, ring_cut) = spectrum_peaks(grid);69    let churn = match previous {70        Some(prev) => churn(&[prev.clone(), frame.clone()]),71        None => 0.0,72    };73    Ok(Reading {74        side,75        fill,76        components: pieces,77        holes,78        euler,79        box_slope: box_slope(grid),80        symmetry: symmetry(grid),81        ring,82        share,83        ring_cut,84        entropy: entropy(frame),85        churn,86        cuts: cuts(grid),87        factors: factors(grid),88    })89}9091/// Splits a square frame at a divisor d of its side: every non-zero d-block must be one tile, then the frame is outer (x) inner.92pub fn block_split(frame: &Tensor, d: usize) -> Option<(Tensor, Tensor)> {93    let side = frame.shape[0];94    if d == 0 || !side.is_multiple_of(d) {95        return None;96    }97    let n = side / d;98    let mut outer = Tensor::new(vec![d, d]);99    let mut inner: Option<Tensor> = None;100    for i in 0..d {101        for j in 0..d {102            let mut block = Tensor::new(vec![n, n]);103            let mut live = false;104            for p in 0..n {105                for q in 0..n {106                    if frame.get(&[i * n + p, j * n + q]).ok()? != 0 {107                        block.set(&[p, q], 1).ok()?;108                        live = true;109                    }110                }111            }112            if !live {113                continue;114            }115            outer.set(&[i, j], 1).ok()?;116            match &inner {117                None => inner = Some(block),118                Some(first) => {119                    if *first != block {120                        return None;121                    }122                }123            }124        }125    }126    inner.map(|block| (outer, block))127}128129/// Returns every proper divisor of the side at which the frame splits.130pub fn cuts(frame: &Tensor) -> Vec<usize> {131    let side = frame.shape[0];132    (2..side)133        .filter(|&d| side.is_multiple_of(d) && block_split(frame, d).is_some())134        .collect()135}136137/// Packs a 0/1 tile into its row-major code, None past 128 cells.138pub fn pack(tile: &Tensor) -> Option<u128> {139    if tile.size() > 128 {140        return None;141    }142    let mut code: u128 = 0;143    for i in 0..tile.size() {144        if tile.at(i) != 0 {145            code |= 1 << i;146        }147    }148    Some(code)149}150151/// Returns the smallest cut whose outer and inner factors both pack, as (d, outer code, inner code).152pub fn factors(frame: &Tensor) -> Option<(usize, u128, u128)> {153    for d in cuts(frame) {154        if let Some((outer, inner)) = block_split(frame, d) {155            if let (Some(a), Some(b)) = (pack(&outer), pack(&inner)) {156                return Some((d, a, b));157            }158        }159    }160    None161}162163/// Returns the first torus shift and cut at which the shifted frame factors, as (dr, dc, d, outer, inner).164pub fn shifted_factors(grid: &Tensor) -> Option<(usize, usize, usize, u128, u128)> {165    let n = grid.shape[0];166    if grid.sum() == 0 {167        return None;168    }169    let mut shifted = Tensor::new(vec![n, n]);170    for dr in 0..n {171        for dc in 0..n {172            for r in 0..n {173                for c in 0..n {174                    let v = grid.at(((r + dr) % n) * n + (c + dc) % n);175                    shifted.put(r * n + c, v);176                }177            }178            if let Some((d, a, b)) = factors(&shifted) {179                return Some((dr, dc, d, a, b));180            }181        }182    }183    None184}185186fn images(grid: &Tensor) -> Vec<Tensor> {187    let mut out = Vec::with_capacity(8);188    for k in 0..4 {189        let Ok(turned) = grid.rot90(k, (0, 1)) else {190            break;191        };192        if let Ok(flipped) = turned.flip(1) {193            out.push(flipped);194        }195        out.push(turned);196    }197    out198}199200/// Returns the dihedral subgroup fixing the frame.201pub fn symmetry(grid: &Tensor) -> Symmetry {202    let same = |t: Result<Tensor>| t.is_ok_and(|image| image.bytes().ok() == grid.bytes().ok());203    let r90 = same(grid.rot90(1, (0, 1)));204    let r180 = same(grid.rot90(2, (0, 1)));205    let fh = same(grid.flip(1));206    let fv = same(grid.flip(0));207    let t = same(grid.transpose(0, 1));208    let at = same(grid.rot90(2, (0, 1)).and_then(|t| t.transpose(0, 1)));209    let order = 1 + [r90, r180, r90, fh, fv, t, at]210        .iter()211        .filter(|&&b| b)212        .count();213    let name = match (order, r90, fh, t, r180) {214        (8, ..) => "d4",215        (4, true, ..) => "c4",216        (4, _, true, ..) => "d2",217        (4, ..) => "d2d",218        (2, _, _, _, true) => "c2",219        (2, ..) if fh || fv => "m",220        (2, ..) => "md",221        _ => "1",222    };223    Symmetry { order, name }224}225226/// Returns the least-squares slope of ln(box count) against ln(side over box) over dyadic boxes, zero on an empty frame.227pub fn box_slope(grid: &Tensor) -> f64 {228    let side = grid.shape[0];229    if grid.sum() == 0 || side < 2 {230        return 0.0;231    }232    let mut points = Vec::new();233    let mut b = 1;234    while b < side {235        let m = side.div_ceil(b);236        let mut boxes = vec![false; m * m];237        for r in 0..side {238            for c in 0..side {239                if grid.at(grid.index(&[r, c])) != 0 {240                    boxes[(r / b) * m + c / b] = true;241                }242            }243        }244        let count = boxes.iter().filter(|&&x| x).count();245        points.push(((side as f64 / b as f64).ln(), (count as f64).ln()));246        b *= 2;247    }248    let n = points.len() as f64;249    let (sx, sy): (f64, f64) = points250        .iter()251        .fold((0.0, 0.0), |(a, b), (x, y)| (a + x, b + y));252    let (sxx, sxy): (f64, f64) = points253        .iter()254        .fold((0.0, 0.0), |(a, b), (x, y)| (a + x * x, b + x * y));255    let denominator = n * sxx - sx * sx;256    if denominator.abs() < 1e-12 {257        return 0.0;258    }259    (n * sxy - sx * sy) / denominator260}261262/// Transforms a real n-square field exactly, rows then columns, returning the real and imaginary parts.263pub fn dft2(field: &[f64], n: usize) -> (Vec<f64>, Vec<f64>) {264    let cos: Vec<f64> = (0..n)265        .map(|k| (2.0 * PI * k as f64 / n as f64).cos())266        .collect();267    let sin: Vec<f64> = (0..n)268        .map(|k| (2.0 * PI * k as f64 / n as f64).sin())269        .collect();270    let mut re = vec![0.0; n * n];271    let mut im = vec![0.0; n * n];272    for r in 0..n {273        for kc in 0..n {274            let (mut a, mut b) = (0.0, 0.0);275            for c in 0..n {276                let v = field[r * n + c];277                if v != 0.0 {278                    let phase = (kc * c) % n;279                    a += v * cos[phase];280                    b -= v * sin[phase];281                }282            }283            re[r * n + kc] = a;284            im[r * n + kc] = b;285        }286    }287    let mut out_re = vec![0.0; n * n];288    let mut out_im = vec![0.0; n * n];289    for kc in 0..n {290        for kr in 0..n {291            let (mut a, mut b) = (0.0, 0.0);292            for r in 0..n {293                let phase = (kr * r) % n;294                let (x, y) = (re[r * n + kc], im[r * n + kc]);295                a += x * cos[phase] + y * sin[phase];296                b += y * cos[phase] - x * sin[phase];297            }298            out_re[kr * n + kc] = a;299            out_im[kr * n + kc] = b;300        }301    }302    (out_re, out_im)303}304305fn centred(values: &[f64], n: usize) -> Vec<f64> {306    let half = n / 2;307    let mut out = vec![0.0; n * n];308    for r in 0..n {309        for c in 0..n {310            out[((r + half) % n) * n + (c + half) % n] = values[r * n + c];311        }312    }313    out314}315316/// Returns the centred power square of a frame on its own torus.317pub fn power_square(grid: &Tensor) -> Vec<f64> {318    let n = grid.shape[0];319    let field: Vec<f64> = (0..grid.size()).map(|i| grid.at(i) as f64).collect();320    let (re, im) = dft2(&field, n);321    let power: Vec<f64> = re.iter().zip(&im).map(|(a, b)| a * a + b * b).collect();322    centred(&power, n)323}324325fn spectrum_peaks(grid: &Tensor) -> (usize, f64, usize) {326    let n = grid.shape[0];327    let power = power_square(grid);328    let half = n / 2;329    let ring_of = |r: usize, c: usize| {330        let dr = r as f64 - half as f64;331        let dc = c as f64 - half as f64;332        dr.hypot(dc).round() as usize333    };334    let floor = 1e-9 * power[half * n + half];335    let mut best = 0.0;336    let mut ring = 0;337    let mut total = 0.0;338    for r in 0..n {339        for c in 0..n {340            if r == half && c == half {341                continue;342            }343            let p = power[r * n + c];344            if p <= floor {345                continue;346            }347            total += p;348            if p > best {349                best = p;350                ring = ring_of(r, c);351            }352        }353    }354    if ring == 0 || total <= 0.0 {355        return (0, 0.0, 0);356    }357    let mut on_ring = 0.0;358    for r in 0..n {359        for c in 0..n {360            if !(r == half && c == half) && ring_of(r, c) == ring && power[r * n + c] > floor {361                on_ring += power[r * n + c];362            }363        }364    }365    let cut = peak_ring(&radial_profile(&power, n).expect("the power square is n by n"));366    (ring, on_ring / total, cut)367}368369/// Returns the first ring where the ring mean of a mask's signed transform on the canvas torus turns negative, zero when it never does.370pub fn first_negative_lobe(mask: &Tensor, canvas: usize) -> usize {371    let side = mask.shape[0];372    let centre = (side - 1) / 2;373    let mut field = vec![0.0; canvas * canvas];374    for r in 0..side {375        for c in 0..side {376            if mask.at(mask.index(&[r, c])) != 0 {377                let rr = (r + canvas - centre % canvas) % canvas;378                let cc = (c + canvas - centre % canvas) % canvas;379                field[rr * canvas + cc] += 1.0;380            }381        }382    }383    let (re, _) = dft2(&field, canvas);384    let profile =385        radial_profile(&centred(&re, canvas), canvas).expect("the field is canvas by canvas");386    profile387        .iter()388        .enumerate()389        .skip(1)390        .find(|(_, &v)| v < 0.0)391        .map(|(k, _)| k)392        .unwrap_or(0)393}394395/// Returns the frame's canonical bytes under the dihedral group and torus translation.396pub fn canonical(grid: &Tensor) -> Vec<u8> {397    let n = grid.shape[0];398    if grid.sum() == 0 {399        return grid.bytes().map(<[u8]>::to_vec).unwrap_or_default();400    }401    let mut best: Option<Vec<u8>> = None;402    let mut shifted = vec![0u8; n * n];403    for image in images(grid) {404        let Ok(bytes) = image.bytes() else { continue };405        for r0 in 0..n {406            for c0 in 0..n {407                if bytes[r0 * n + c0] == 0 {408                    continue;409                }410                for r in 0..n {411                    let rr = (r + n - r0) % n;412                    for c in 0..n {413                        shifted[rr * n + (c + n - c0) % n] = bytes[r * n + c];414                    }415                }416                if best.as_ref().is_none_or(|b| shifted < *b) {417                    best = Some(shifted.clone());418                }419            }420        }421    }422    best.unwrap_or_default()423}424425/// Returns the torus shift carrying frame a onto frame b, when one exists.426pub fn translate_of(a: &Tensor, b: &Tensor) -> Option<(usize, usize)> {427    let n = a.shape[0];428    if a.shape != b.shape || a.sum() != b.sum() {429        return None;430    }431    let (x, y) = (a.bytes().ok()?, b.bytes().ok()?);432    let first = x.iter().position(|&v| v != 0)?;433    let (r0, c0) = (first / n, first % n);434    for r1 in 0..n {435        for c1 in 0..n {436            if y[r1 * n + c1] == 0 {437                continue;438            }439            let (dr, dc) = ((r1 + n - r0) % n, (c1 + n - c0) % n);440            let fits = (0..n * n).all(|i| {441                let (r, c) = (i / n, i % n);442                x[i] == y[((r + dr) % n) * n + (c + dc) % n]443            });444            if fits {445                return Some((dr, dc));446            }447        }448    }449    None450}451452#[cfg(test)]453mod tests {454    use super::*;455    use mrlyrs::math::two::carpet;456    fn blinker() -> Cell2d {457        let mut t = Tensor::new(vec![5, 5]);458        t.set(&[1, 2], 1).unwrap();459        t.set(&[2, 2], 1).unwrap();460        t.set(&[3, 2], 1).unwrap();461        Cell2d::new(t).unwrap()462    }463    #[test]464    fn the_blinker_reads_one_bar_of_d2_with_no_cut() {465        let reading = read(&blinker(), None).unwrap();466        assert_eq!(467            (468                reading.fill,469                reading.components,470                reading.holes,471                reading.euler472            ),473            (3, 1, 1, 1)474        );475        assert_eq!(476            reading.symmetry,477            Symmetry {478                order: 4,479                name: "d2"480            }481        );482        assert!(reading.cuts.is_empty());483        assert_eq!(reading.factors, None);484        assert!(reading.ring > 0 && reading.share > 0.0);485        assert_eq!(reading.entropy, entropy(&blinker()));486        let turned = Cell2d::new(blinker().types().rot90(1, (0, 1)).unwrap()).unwrap();487        assert_eq!(canonical(blinker().types()), canonical(turned.types()));488        assert_eq!(read(&turned, Some(&blinker())).unwrap().churn, 4.0 / 25.0);489    }490    #[test]491    fn the_full_block_is_d4_with_a_flat_spectrum_and_cuts_at_two() {492        let block = Cell2d::new(Tensor::full(vec![4, 4], 1)).unwrap();493        let reading = read(&block, None).unwrap();494        assert_eq!(495            (496                reading.fill,497                reading.components,498                reading.holes,499                reading.euler500            ),501            (16, 1, 0, 1)502        );503        assert_eq!(504            reading.symmetry,505            Symmetry {506                order: 8,507                name: "d4"508            }509        );510        assert_eq!((reading.ring, reading.share, reading.ring_cut), (0, 0.0, 0));511        assert_eq!(reading.entropy, 0);512        assert_eq!(reading.cuts, vec![2]);513        assert_eq!(reading.factors, Some((2, 0b1111, 0b1111)));514        assert!((reading.box_slope - 2.0).abs() < 1e-9);515    }516    #[test]517    fn the_level_two_carpet_factors_at_three_into_two_carpet_tiles() {518        let frame = carpet(3, 2).unwrap();519        let reading = read(&frame, None).unwrap();520        assert_eq!(521            (522                reading.fill,523                reading.components,524                reading.holes,525                reading.euler526            ),527            (64, 1, 9, -8)528        );529        assert_eq!(reading.symmetry.order, 8);530        assert_eq!(reading.cuts, vec![3]);531        assert_eq!(reading.factors, Some((3, 495, 495)));532        let (outer, inner) = block_split(frame.types(), 3).unwrap();533        assert_eq!(534            outer.kron(&inner).bytes().unwrap(),535            frame.types().bytes().unwrap()536        );537        assert_eq!(538            block_split(&frame.types().rot90(1, (0, 1)).unwrap(), 3).map(|(a, _)| pack(&a)),539            Some(Some(495))540        );541    }542    #[test]543    fn a_shifted_frame_is_found_and_read_back() {544        let a = blinker();545        let mut t = Tensor::new(vec![5, 5]);546        t.set(&[4, 0], 1).unwrap();547        t.set(&[0, 0], 1).unwrap();548        t.set(&[1, 0], 1).unwrap();549        let b = Cell2d::new(t).unwrap();550        assert_eq!(translate_of(a.types(), b.types()), Some((3, 3)));551        assert_eq!(translate_of(a.types(), a.types()), Some((0, 0)));552        assert_eq!(canonical(a.types()), canonical(b.types()));553        let carpet = mrlyrs::math::two::carpet(3, 2).unwrap();554        let moved =555            Cell2d::new(Tensor::of(canonical(carpet.types()), vec![9, 9]).unwrap()).unwrap();556        assert_eq!(factors(moved.types()), None);557        let (_, _, d, outer, inner) = shifted_factors(moved.types()).unwrap();558        assert_eq!((d, inner), (3, 495));559        let tile = |code: u128| {560            Tensor::of(561                (0..9).map(|i| ((code >> i) & 1) as u8).collect(),562                vec![3, 3],563            )564            .unwrap()565        };566        assert!(translate_of(&tile(495), &tile(outer)).is_some());567        assert_eq!(568            first_negative_lobe(&mrlyrs::life::moore().unwrap().types().clone(), 27),569            9570        );571    }572}