census.rs

1use super::graph::edge_graph;2use super::models::Cell3d;3use crate::dim::census;4use mrlycore::errors::Result;56pub use crate::dim::census::{edges, vertices};78/// The tally of a cube's sites.9#[derive(Clone, Debug, PartialEq, Eq)]10pub struct Census {11    /// The count of filled sites.12    pub fills: usize,13    /// The count of empty sites.14    pub voids: usize,15    /// The count of exposed faces.16    pub surface: u128,17    /// The count of corners the filled sites touch.18    pub vertices: usize,19    /// The count of unit edges the filled sites touch.20    pub edges: usize,21    /// The count of unit faces the filled sites touch, shared ones counted once.22    pub faces: usize,23    /// The Euler characteristic of the filled complex.24    pub euler: i64,25}2627/// Returns the count of filled sites.28pub fn fills(cell: &Cell3d) -> usize {29    census::fills(cell)30}3132/// Returns the count of empty sites.33pub fn voids(cell: &Cell3d) -> usize {34    census::voids(cell)35}3637/// Returns the filled-site count, the cube's volume.38pub fn volume(cell: &Cell3d) -> usize {39    fills(cell)40}4142/// Returns the count of filled faces exposed to void or the outside.43pub fn surface(cell: &Cell3d) -> u128 {44    census::exposure(cell)45}4647/// Returns the count of faces buried between two filled sites, six per site less the exposed surface.48///49/// ```50/// let block = mrlymath::three::ones(2, 1).unwrap();51/// assert_eq!(mrlymath::three::census::hidden(&block), 24);52/// ```53pub fn hidden(cell: &Cell3d) -> u128 {54    6 * fills(cell) as u128 - surface(cell)55}5657/// Returns the count of unit faces the filled sites touch, a face shared by two sites counted once.58///59/// Surface counts only the faces open to void; this counts every face of the complex.60///61/// ```62/// let block = mrlymath::three::ones(2, 1).unwrap();63/// assert_eq!(mrlymath::three::census::faces(&block), 36);64/// assert_eq!(mrlymath::three::census::surface(&block), 24);65/// ```66pub fn faces(cell: &Cell3d) -> usize {67    let grid = cell.types();68    let (dx, dy, dz) = (grid.shape[0], grid.shape[1], grid.shape[2]);69    let mut planes = [70        vec![false; (dx + 1) * dy * dz],71        vec![false; dx * (dy + 1) * dz],72        vec![false; dx * dy * (dz + 1)],73    ];74    for i in 0..dx {75        for j in 0..dy {76            for k in 0..dz {77                if grid.get(&[i, j, k]) == 0 {78                    continue;79                }80                for step in 0..2 {81                    planes[0][((i + step) * dy + j) * dz + k] = true;82                    planes[1][(i * (dy + 1) + j + step) * dz + k] = true;83                    planes[2][(i * dy + j) * (dz + 1) + k + step] = true;84                }85            }86        }87    }88    planes.iter().flatten().filter(|&&on| on).count()89}9091/// Returns the Euler characteristic of the filled complex, vertices less edges plus faces less sites.92///93/// One for a solid block, and one less than twice the genus below zero for a tunnelled one.94///95/// ```96/// assert_eq!(mrlymath::three::census::euler(&mrlymath::three::ones(2, 1).unwrap()).unwrap(), 1);97/// assert_eq!(mrlymath::three::census::euler(&mrlymath::three::carpet(3, 1).unwrap()).unwrap(), -4);98/// ```99pub fn euler(cell: &Cell3d) -> Result<i64> {100    let net = edge_graph(cell)?;101    let (v, e) = (net.nodes.len() as i64, net.branches.len() as i64);102    Ok(v - e + faces(cell) as i64 - fills(cell) as i64)103}104105/// Tallies a cell's sites, its exposed surface and its Euler characteristic in one reading.106///107/// ```108/// let tally = mrlymath::three::census::census(&mrlymath::three::carpet(3, 1).unwrap()).unwrap();109/// assert_eq!((tally.fills, tally.voids, tally.surface), (20, 7, 72));110/// assert_eq!((tally.vertices, tally.edges, tally.faces, tally.euler), (64, 144, 96, -4));111/// ```112pub fn census(cell: &Cell3d) -> Result<Census> {113    let net = edge_graph(cell)?;114    let (vertices, edges) = (net.nodes.len(), net.branches.len());115    let (fills, faces) = (fills(cell), faces(cell));116    Ok(Census {117        fills,118        voids: voids(cell),119        surface: surface(cell),120        vertices,121        edges,122        faces,123        euler: vertices as i64 - edges as i64 + faces as i64 - fills as i64,124    })125}126127#[cfg(test)]128mod tests {129    use super::*;130    use crate::formulas;131    use crate::three::designs;132    #[test]133    fn census_matches_formulas() {134        for code in [23u128, 129, 17, 232] {135            for level in 1..3u32 {136                let cell = designs::create(code, 3, level as usize, 2).unwrap();137                assert_eq!(138                    fills(&cell) as u128,139                    formulas::fill(code, 3, 3, level, 2).unwrap()140                );141                assert_eq!(142                    surface(&cell),143                    formulas::surface(code, 3, level, 2).unwrap(),144                    "code={code} l={level}"145                );146            }147        }148    }149    #[test]150    fn menger_census() {151        let c = designs::carpet(3, 1).unwrap();152        let result = census(&c).unwrap();153        assert_eq!(result.fills, 20);154        assert_eq!(result.voids, 7);155        assert_eq!(result.surface, 72);156        assert_eq!(result.vertices, 64);157        assert_eq!(result.edges, 144);158        assert_eq!(result.faces, 96);159        assert_eq!(result.euler, -4);160    }161    #[test]162    fn a_solid_block_is_contractible() {163        for n in [1, 2, 3] {164            let block = designs::ones(n, 1).unwrap();165            assert_eq!(euler(&block).unwrap(), 1, "n={n}");166        }167        let one = census(&designs::ones(1, 1).unwrap()).unwrap();168        assert_eq!((one.vertices, one.edges, one.faces), (8, 12, 6));169    }170    #[test]171    fn the_sponge_deepens_its_genus() {172        let level_two = census(&designs::carpet(3, 2).unwrap()).unwrap();173        assert_eq!(level_two.fills, 400);174        assert_eq!(level_two.vertices, 896);175        assert_eq!(level_two.edges, 2304);176        assert_eq!(level_two.faces, 1728);177        assert_eq!(level_two.euler, -80);178    }179    #[test]180    fn the_parts_agree_with_the_whole() {181        for cell in [182            designs::net(3, 1).unwrap(),183            designs::void(4, 1).unwrap(),184            designs::xtree(3, 2).unwrap(),185        ] {186            let tally = census(&cell).unwrap();187            assert_eq!(tally.vertices, vertices(&cell).unwrap());188            assert_eq!(tally.edges, edges(&cell).unwrap());189            assert_eq!(tally.faces, faces(&cell));190            assert_eq!(tally.euler, euler(&cell).unwrap());191            assert_eq!(tally.fills, volume(&cell));192        }193    }194}195196#[cfg(test)]197mod theorems {198    use super::*;199    use crate::bang::universe::{orbit, total_exposure, touches_every_corner};200    use crate::three::designs;201    use crate::two;202    use std::collections::BTreeSet;203204    fn tile_fill(code: u128, number: usize) -> u128 {205        fills(&designs::create(code, number, 1, 2).unwrap()) as u128206    }207208    fn state(code: u128, number: usize, level: usize) -> (i128, i128) {209        let cell = designs::create(code, number, level, 2).unwrap();210        (surface(&cell) as i128, hidden(&cell) as i128)211    }212213    fn second_eigenvalue(code: u128, number: usize) -> i128 {214        let half = (number / 2) as i128;215        let side = number as i128;216        match code {217            23 => side * side - half * half,218            232 => half * half,219            _ => (side - half) * (side - half),220        }221    }222223    #[test]224    fn the_face_ledger_prints_the_family_closed_forms() {225        let visible = [226            (23u128, [72u128, 1056, 18048, 336384]),227            (232, [30, 198, 1374, 9606]),228            (3, [56, 608, 7040, 83456]),229            (129, [54, 486, 4374, 39366]),230        ];231        for (code, faces) in visible {232            let fc = tile_fill(code, 3);233            for (index, &want) in faces.iter().enumerate() {234                let level = index + 1;235                let cell = designs::create(code, 3, level, 2).unwrap();236                assert_eq!(surface(&cell), want, "code={code} l={level}");237                assert_eq!(238                    hidden(&cell),239                    6 * fc.pow(level as u32) - want,240                    "code={code} l={level}"241                );242            }243        }244        let carpet: Vec<u128> = (1..5)245            .map(|level| hidden(&designs::create(23, 3, level, 2).unwrap()))246            .collect();247        assert_eq!(carpet, [48, 1344, 29952, 623616]);248    }249250    #[test]251    fn the_face_matrix_fits_its_eigenvalues_and_predicts() {252        let mut fitted = 0;253        for (number, top) in [(3usize, 4usize), (5, 3), (7, 2)] {254            for code in [23u128, 232, 3, 129] {255                let fc = tile_fill(code, number) as i128;256                let l2 = second_eigenvalue(code, number);257                let states: Vec<(i128, i128)> =258                    (1..=top).map(|level| state(code, number, level)).collect();259                let work = states[0].1 / 2;260                for (index, &(v, h)) in states.iter().enumerate() {261                    let level = index as u32 + 1;262                    assert_eq!(v + h, 6 * fc.pow(level), "code={code} n={number} l={level}");263                    let want = if work == 0 {264                        0265                    } else {266                        2 * work * (fc.pow(level) - l2.pow(level)) / (fc - l2)267                    };268                    assert_eq!(h, want, "code={code} n={number} l={level}");269                }270                if top < 3 || work == 0 {271                    continue;272                }273                let ((v1, h1), (v2, h2), (v3, h3)) = (states[0], states[1], states[2]);274                let base = v1 * h2 - v2 * h1;275                assert_eq!(276                    v1 * h3 - h1 * v3,277                    (fc + l2) * base,278                    "code={code} n={number}"279                );280                assert_eq!(v2 * h3 - v3 * h2, fc * l2 * base, "code={code} n={number}");281                if code == 23 && number == 3 {282                    let entries = [283                        v2 * h2 - v3 * h1,284                        v1 * v3 - v2 * v2,285                        h2 * h2 - h3 * h1,286                        v1 * h3 - v2 * h2,287                    ];288                    for entry in entries {289                        assert_eq!(entry % base, 0, "code={code} n={number}");290                    }291                    let matrix = [292                        [entries[0] / base, entries[1] / base],293                        [entries[2] / base, entries[3] / base],294                    ];295                    assert_eq!(matrix, [[12, 4], [8, 16]]);296                }297                let numerator = (v2 * h2 - v3 * h1) * v3 + (v3 * v1 - v2 * v2) * h3;298                assert_eq!(numerator % base, 0, "code={code} n={number}");299                let next = numerator / base;300                let closed = 6 * fc.pow(4) - 2 * work * (fc.pow(4) - l2.pow(4)) / (fc - l2);301                assert_eq!(next, closed, "code={code} n={number}");302                if number == 3 {303                    assert_eq!(next, states[3].0, "code={code}");304                }305                fitted += 1;306            }307        }308        assert_eq!(fitted, 6);309    }310311    #[test]312    fn the_void_buries_no_face() {313        for k in 1..13u128 {314            let number = 2 * k as usize - 1;315            let cell = designs::create(129, number, 1, 2).unwrap();316            let cells = k * k * k + (k - 1) * (k - 1) * (k - 1);317            assert_eq!(fills(&cell) as u128, cells, "k={k}");318            assert_eq!(surface(&cell), 6 * cells, "k={k}");319            assert_eq!(hidden(&cell), 0, "k={k}");320            let flat = two::create(9, number, 1, 0, 2).unwrap();321            let tally = two::census::census(&flat).unwrap();322            assert_eq!(tally.edges as u128, tally.perimeter, "k={k}");323        }324        for level in 1..5u32 {325            let cell = designs::create(129, 3, level as usize, 2).unwrap();326            assert_eq!(surface(&cell), 6 * 9u128.pow(level), "l={level}");327        }328    }329330    #[test]331    fn total_exposure_holds_for_the_independent_corner_sets() {332        for number in [3usize, 5, 7] {333            for code in 0..256u128 {334                let cell = designs::create(code, number, 1, 2).unwrap();335                let open = surface(&cell) == 6 * fills(&cell) as u128;336                assert_eq!(open, total_exposure(code, 3), "code={code} n={number}");337            }338        }339        let exposed: Vec<u128> = (0..256).filter(|&c| total_exposure(c, 3)).collect();340        let classes: BTreeSet<u128> = exposed341            .iter()342            .map(|&c| *orbit(c, 3).iter().next().unwrap())343            .collect();344        assert_eq!(exposed.len(), 35);345        assert_eq!(346            classes.into_iter().collect::<Vec<u128>>(),347            [0, 1, 6, 22, 24, 105]348        );349    }350351    #[test]352    fn the_all_even_rule_touches_every_grid_corner() {353        for k in 1..4usize {354            let number = 2 * k - 1;355            let grid = (number + 1).pow(3);356            for code in 0..256u128 {357                let cell = designs::create(code, number, 1, 2).unwrap();358                let whole = vertices(&cell).unwrap() == grid;359                assert_eq!(whole, touches_every_corner(code, 3), "code={code} k={k}");360            }361        }362        assert_eq!(363            (0..256u128).filter(|&c| touches_every_corner(c, 3)).count(),364            128365        );366    }367368    #[test]369    fn the_net_falls_short_of_the_grid_corners() {370        for k in 1..21usize {371            let cell = designs::create(232, 2 * k - 1, 1, 2).unwrap();372            let m = k - 1;373            assert_eq!(vertices(&cell).unwrap(), 8 * m * m * (k + 2), "k={k}");374            assert_eq!(375                8 * k * k * k - vertices(&cell).unwrap(),376                24 * k - 16,377                "k={k}"378            );379        }380        for k in 1..25usize {381            let flat = two::net(2 * k - 1, 1).unwrap();382            assert_eq!(383                two::census::vertices(&flat).unwrap(),384                4 * k * k - 4,385                "k={k}"386            );387        }388    }389}