readings.rs

1use mrlycore::errors::Result;2use mrlycore::tensor::Tensor;3use mrlymath::dim::graph::{core_graph, tunnel_graph};4use mrlymath::life::{churn, entropy};5use mrlymath::two::{census, Cell2d};6use mrlynum::fft::{peak_ring, radial_profile};7use mrlynum::graph::census::components;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]) != 0 {107                        block.set(&[p, q], 1);108                        live = true;109                    }110                }111            }112            if !live {113                continue;114            }115            outer.set(&[i, j], 1);116            match &inner {117                None => inner = Some(block),118                Some(first) => {119                    if first.bytes() != block.bytes() {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, &b) in tile.bytes().iter().enumerate() {144        if b != 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 bytes = grid.bytes();170    let mut shifted = Tensor::new(vec![n, n]);171    for dr in 0..n {172        for dc in 0..n {173            for r in 0..n {174                for c in 0..n {175                    let v = bytes[((r + dr) % n) * n + (c + dc) % n];176                    shifted.bytes_mut()[r * n + c] = v;177                }178            }179            if let Some((d, a, b)) = factors(&shifted) {180                return Some((dr, dc, d, a, b));181            }182        }183    }184    None185}186187fn images(grid: &Tensor) -> Vec<Tensor> {188    let mut out = Vec::with_capacity(8);189    for k in 0..4 {190        let turned = grid.rot90(k, (0, 1));191        out.push(turned.flip(1));192        out.push(turned);193    }194    out195}196197/// Returns the dihedral subgroup fixing the frame.198pub fn symmetry(grid: &Tensor) -> Symmetry {199    let same = |t: &Tensor| t.bytes() == grid.bytes();200    let r90 = same(&grid.rot90(1, (0, 1)));201    let r180 = same(&grid.rot90(2, (0, 1)));202    let fh = same(&grid.flip(1));203    let fv = same(&grid.flip(0));204    let t = same(&grid.transpose(0, 1));205    let at = same(&grid.rot90(2, (0, 1)).transpose(0, 1));206    let order = 1 + [r90, r180, r90, fh, fv, t, at]207        .iter()208        .filter(|&&b| b)209        .count();210    let name = match (order, r90, fh, t, r180) {211        (8, ..) => "d4",212        (4, true, ..) => "c4",213        (4, _, true, ..) => "d2",214        (4, ..) => "d2d",215        (2, _, _, _, true) => "c2",216        (2, ..) if fh || fv => "m",217        (2, ..) => "md",218        _ => "1",219    };220    Symmetry { order, name }221}222223/// Returns the least-squares slope of ln(box count) against ln(side over box) over dyadic boxes, zero on an empty frame.224pub fn box_slope(grid: &Tensor) -> f64 {225    let side = grid.shape[0];226    if grid.sum() == 0 || side < 2 {227        return 0.0;228    }229    let mut points = Vec::new();230    let mut b = 1;231    while b < side {232        let m = side.div_ceil(b);233        let mut boxes = vec![false; m * m];234        for r in 0..side {235            for c in 0..side {236                if grid.get(&[r, c]) != 0 {237                    boxes[(r / b) * m + c / b] = true;238                }239            }240        }241        let count = boxes.iter().filter(|&&x| x).count();242        points.push(((side as f64 / b as f64).ln(), (count as f64).ln()));243        b *= 2;244    }245    let n = points.len() as f64;246    let (sx, sy): (f64, f64) = points247        .iter()248        .fold((0.0, 0.0), |(a, b), (x, y)| (a + x, b + y));249    let (sxx, sxy): (f64, f64) = points250        .iter()251        .fold((0.0, 0.0), |(a, b), (x, y)| (a + x * x, b + x * y));252    let denominator = n * sxx - sx * sx;253    if denominator.abs() < 1e-12 {254        return 0.0;255    }256    (n * sxy - sx * sy) / denominator257}258259/// Transforms a real n-square field exactly, rows then columns, returning the real and imaginary parts.260pub fn dft2(field: &[f64], n: usize) -> (Vec<f64>, Vec<f64>) {261    let cos: Vec<f64> = (0..n)262        .map(|k| (2.0 * PI * k as f64 / n as f64).cos())263        .collect();264    let sin: Vec<f64> = (0..n)265        .map(|k| (2.0 * PI * k as f64 / n as f64).sin())266        .collect();267    let mut re = vec![0.0; n * n];268    let mut im = vec![0.0; n * n];269    for r in 0..n {270        for kc in 0..n {271            let (mut a, mut b) = (0.0, 0.0);272            for c in 0..n {273                let v = field[r * n + c];274                if v != 0.0 {275                    let phase = (kc * c) % n;276                    a += v * cos[phase];277                    b -= v * sin[phase];278                }279            }280            re[r * n + kc] = a;281            im[r * n + kc] = b;282        }283    }284    let mut out_re = vec![0.0; n * n];285    let mut out_im = vec![0.0; n * n];286    for kc in 0..n {287        for kr in 0..n {288            let (mut a, mut b) = (0.0, 0.0);289            for r in 0..n {290                let phase = (kr * r) % n;291                let (x, y) = (re[r * n + kc], im[r * n + kc]);292                a += x * cos[phase] + y * sin[phase];293                b += y * cos[phase] - x * sin[phase];294            }295            out_re[kr * n + kc] = a;296            out_im[kr * n + kc] = b;297        }298    }299    (out_re, out_im)300}301302fn centred(values: &[f64], n: usize) -> Vec<f64> {303    let half = n / 2;304    let mut out = vec![0.0; n * n];305    for r in 0..n {306        for c in 0..n {307            out[((r + half) % n) * n + (c + half) % n] = values[r * n + c];308        }309    }310    out311}312313/// Returns the centred power square of a frame on its own torus.314pub fn power_square(grid: &Tensor) -> Vec<f64> {315    let n = grid.shape[0];316    let field: Vec<f64> = grid.bytes().iter().map(|&b| b as f64).collect();317    let (re, im) = dft2(&field, n);318    let power: Vec<f64> = re.iter().zip(&im).map(|(a, b)| a * a + b * b).collect();319    centred(&power, n)320}321322fn spectrum_peaks(grid: &Tensor) -> (usize, f64, usize) {323    let n = grid.shape[0];324    let power = power_square(grid);325    let half = n / 2;326    let ring_of = |r: usize, c: usize| {327        let dr = r as f64 - half as f64;328        let dc = c as f64 - half as f64;329        dr.hypot(dc).round() as usize330    };331    let floor = 1e-9 * power[half * n + half];332    let mut best = 0.0;333    let mut ring = 0;334    let mut total = 0.0;335    for r in 0..n {336        for c in 0..n {337            if r == half && c == half {338                continue;339            }340            let p = power[r * n + c];341            if p <= floor {342                continue;343            }344            total += p;345            if p > best {346                best = p;347                ring = ring_of(r, c);348            }349        }350    }351    if ring == 0 || total <= 0.0 {352        return (0, 0.0, 0);353    }354    let mut on_ring = 0.0;355    for r in 0..n {356        for c in 0..n {357            if !(r == half && c == half) && ring_of(r, c) == ring && power[r * n + c] > floor {358                on_ring += power[r * n + c];359            }360        }361    }362    let cut = peak_ring(&radial_profile(&power, n));363    (ring, on_ring / total, cut)364}365366/// 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.367pub fn first_negative_lobe(mask: &Tensor, canvas: usize) -> usize {368    let side = mask.shape[0];369    let centre = (side - 1) / 2;370    let mut field = vec![0.0; canvas * canvas];371    for r in 0..side {372        for c in 0..side {373            if mask.get(&[r, c]) != 0 {374                let rr = (r + canvas - centre % canvas) % canvas;375                let cc = (c + canvas - centre % canvas) % canvas;376                field[rr * canvas + cc] += 1.0;377            }378        }379    }380    let (re, _) = dft2(&field, canvas);381    let profile = radial_profile(&centred(&re, canvas), canvas);382    profile383        .iter()384        .enumerate()385        .skip(1)386        .find(|(_, &v)| v < 0.0)387        .map(|(k, _)| k)388        .unwrap_or(0)389}390391/// Returns the frame's canonical bytes under the dihedral group and torus translation.392pub fn canonical(grid: &Tensor) -> Vec<u8> {393    let n = grid.shape[0];394    if grid.sum() == 0 {395        return grid.bytes().to_vec();396    }397    let mut best: Option<Vec<u8>> = None;398    let mut shifted = vec![0u8; n * n];399    for image in images(grid) {400        let bytes = image.bytes();401        for r0 in 0..n {402            for c0 in 0..n {403                if bytes[r0 * n + c0] == 0 {404                    continue;405                }406                for r in 0..n {407                    let rr = (r + n - r0) % n;408                    for c in 0..n {409                        shifted[rr * n + (c + n - c0) % n] = bytes[r * n + c];410                    }411                }412                if best.as_ref().is_none_or(|b| shifted < *b) {413                    best = Some(shifted.clone());414                }415            }416        }417    }418    best.unwrap_or_default()419}420421/// Returns the torus shift carrying frame a onto frame b, when one exists.422pub fn translate_of(a: &Tensor, b: &Tensor) -> Option<(usize, usize)> {423    let n = a.shape[0];424    if a.shape != b.shape || a.sum() != b.sum() {425        return None;426    }427    let (x, y) = (a.bytes(), b.bytes());428    let first = x.iter().position(|&v| v != 0)?;429    let (r0, c0) = (first / n, first % n);430    for r1 in 0..n {431        for c1 in 0..n {432            if y[r1 * n + c1] == 0 {433                continue;434            }435            let (dr, dc) = ((r1 + n - r0) % n, (c1 + n - c0) % n);436            let fits = (0..n * n).all(|i| {437                let (r, c) = (i / n, i % n);438                x[i] == y[((r + dr) % n) * n + (c + dc) % n]439            });440            if fits {441                return Some((dr, dc));442            }443        }444    }445    None446}447448#[cfg(test)]449mod tests {450    use super::*;451    use mrlymath::two::carpet;452    fn blinker() -> Cell2d {453        let mut t = Tensor::new(vec![5, 5]);454        t.set(&[1, 2], 1);455        t.set(&[2, 2], 1);456        t.set(&[3, 2], 1);457        Cell2d::new(t)458    }459    #[test]460    fn the_blinker_reads_one_bar_of_d2_with_no_cut() {461        let reading = read(&blinker(), None).unwrap();462        assert_eq!(463            (464                reading.fill,465                reading.components,466                reading.holes,467                reading.euler468            ),469            (3, 1, 1, 1)470        );471        assert_eq!(472            reading.symmetry,473            Symmetry {474                order: 4,475                name: "d2"476            }477        );478        assert!(reading.cuts.is_empty());479        assert_eq!(reading.factors, None);480        assert!(reading.ring > 0 && reading.share > 0.0);481        assert_eq!(reading.entropy, entropy(&blinker()));482        let turned = Cell2d::new(blinker().types().rot90(1, (0, 1)));483        assert_eq!(canonical(blinker().types()), canonical(turned.types()));484        assert_eq!(read(&turned, Some(&blinker())).unwrap().churn, 4.0 / 25.0);485    }486    #[test]487    fn the_full_block_is_d4_with_a_flat_spectrum_and_cuts_at_two() {488        let block = Cell2d::new(Tensor::full(vec![4, 4], 1));489        let reading = read(&block, None).unwrap();490        assert_eq!(491            (492                reading.fill,493                reading.components,494                reading.holes,495                reading.euler496            ),497            (16, 1, 0, 1)498        );499        assert_eq!(500            reading.symmetry,501            Symmetry {502                order: 8,503                name: "d4"504            }505        );506        assert_eq!((reading.ring, reading.share, reading.ring_cut), (0, 0.0, 0));507        assert_eq!(reading.entropy, 0);508        assert_eq!(reading.cuts, vec![2]);509        assert_eq!(reading.factors, Some((2, 0b1111, 0b1111)));510        assert!((reading.box_slope - 2.0).abs() < 1e-9);511    }512    #[test]513    fn the_level_two_carpet_factors_at_three_into_two_carpet_tiles() {514        let frame = carpet(3, 2).unwrap();515        let reading = read(&frame, None).unwrap();516        assert_eq!(517            (518                reading.fill,519                reading.components,520                reading.holes,521                reading.euler522            ),523            (64, 1, 9, -8)524        );525        assert_eq!(reading.symmetry.order, 8);526        assert_eq!(reading.cuts, vec![3]);527        assert_eq!(reading.factors, Some((3, 495, 495)));528        let (outer, inner) = block_split(frame.types(), 3).unwrap();529        assert_eq!(outer.kron(&inner).bytes(), frame.types().bytes());530        assert_eq!(531            block_split(&frame.types().rot90(1, (0, 1)), 3).map(|(a, _)| pack(&a)),532            Some(Some(495))533        );534    }535    #[test]536    fn a_shifted_frame_is_found_and_read_back() {537        let a = blinker();538        let mut t = Tensor::new(vec![5, 5]);539        t.set(&[4, 0], 1);540        t.set(&[0, 0], 1);541        t.set(&[1, 0], 1);542        let b = Cell2d::new(t);543        assert_eq!(translate_of(a.types(), b.types()), Some((3, 3)));544        assert_eq!(translate_of(a.types(), a.types()), Some((0, 0)));545        assert_eq!(canonical(a.types()), canonical(b.types()));546        let carpet = mrlymath::two::carpet(3, 2).unwrap();547        let moved = Cell2d::new(Tensor::of(canonical(carpet.types()), vec![9, 9]));548        assert_eq!(factors(moved.types()), None);549        let (_, _, d, outer, inner) = shifted_factors(moved.types()).unwrap();550        assert_eq!((d, inner), (3, 495));551        let tile = |code: u128| {552            Tensor::of(553                (0..9).map(|i| ((code >> i) & 1) as u8).collect(),554                vec![3, 3],555            )556        };557        assert!(translate_of(&tile(495), &tile(outer)).is_some());558        assert_eq!(559            first_negative_lobe(&mrlymath::life::moore().types().clone(), 27),560            9561        );562    }563}