census.rs

12.6 kB · rust · 374 lines

1use super::models::Cell6d;2use super::{FILL, GRID, VOID};3use std::collections::{BTreeMap, BTreeSet};45type Point = (i64, i64);6type Edge = (Point, Point);78fn north(x: i64, y: i64) -> [Point; 3] {9    [(x, 2 * y + 2), (x + 1, 2 * y), (x + 2, 2 * y + 2)]10}1112fn south(x: i64, y: i64) -> [Point; 3] {13    [(x, 2 * y), (x + 1, 2 * y + 2), (x + 2, 2 * y)]14}1516/// Returns the three corner points of the triangle at x, y under the start parity.17pub fn corners(x: i64, y: i64, start: i64) -> [Point; 3] {18    if (x + y + start).rem_euclid(2) == 0 {19        north(x, y)20    } else {21        south(x, y)22    }23}2425/// Returns a triangle's three edges, each sorted low corner first.26pub fn edges_of(c: &[Point; 3]) -> [Edge; 3] {27    let sorted = |a: Point, b: Point| if a <= b { (a, b) } else { (b, a) };28    [sorted(c[0], c[1]), sorted(c[1], c[2]), sorted(c[0], c[2])]29}3031/// The tally of a triangle mesh.32#[derive(Clone, Debug, PartialEq, Eq)]33pub struct Census {34    /// The count of tallied triangles.35    pub triangles: usize,36    /// The count of filled triangles.37    pub fills: usize,38    /// The count of void triangles.39    pub voids: usize,40    /// The count of backdrop triangles.41    pub grids: usize,42    /// The count of distinct corners.43    pub vertices: usize,44    /// The count of distinct edges.45    pub edges: usize,46    /// The count of edges touching one triangle.47    pub boundary_edges: usize,48    /// The count of edges shared by two triangles.49    pub interior_edges: usize,50    /// The Euler characteristic of the mesh.51    pub euler: i64,52}5354#[derive(Default)]55struct Mesh {56    vertices: BTreeSet<Point>,57    edges: BTreeMap<Edge, usize>,58}5960impl Mesh {61    fn add(&mut self, x: i64, y: i64, start: i64) {62        let c = corners(x, y, start);63        for p in c {64            self.vertices.insert(p);65        }66        for e in edges_of(&c) {67            *self.edges.entry(e).or_insert(0) += 1;68        }69    }70    fn tally(&self, triangles: usize, fills: usize, voids: usize, grids: usize) -> Census {71        let edges = self.edges.len();72        let boundary = self.edges.values().filter(|&&n| n == 1).count();73        Census {74            triangles,75            fills,76            voids,77            grids,78            vertices: self.vertices.len(),79            edges,80            boundary_edges: boundary,81            interior_edges: edges - boundary,82            euler: self.vertices.len() as i64 - edges as i64 + triangles as i64,83        }84    }85}8687/// Tallies a cell's triangles, corners and edges, counting the backdrop only on request.88pub fn census(cell: &Cell6d, include_grid: bool) -> Census {89    let inner = &cell.cell;90    let start = cell.start as i64;91    let (height, width) = (inner.height(), inner.width());92    let (mut fills, mut voids, mut grids) = (0, 0, 0);93    let mut mesh = Mesh::default();94    for y in 0..height {95        for x in 0..width {96            let v = inner.types().get(&[y, x]);97            match v {98                FILL => fills += 1,99                VOID => voids += 1,100                GRID => grids += 1,101                _ => {}102            }103            if v == GRID && !include_grid {104                continue;105            }106            mesh.add(x as i64, y as i64, start);107        }108    }109    let triangles = fills + voids + if include_grid { grids } else { 0 };110    mesh.tally(triangles, fills, voids, grids)111}112113/// Counts the filled triangles of the cell.114///115/// ```116/// use mrlymath::six::{blank, Cell6d, Orientation, Projection, FILL, VOID};117/// let hex = blank(2, Orientation::Horizontal, FILL, VOID);118/// let cell = Cell6d::new(hex, Projection::Cut, Orientation::Horizontal, 0);119/// assert_eq!(mrlymath::six::census::fills(&cell), 24);120/// ```121pub fn fills(cell: &Cell6d) -> usize {122    cell.cell123        .types()124        .bytes()125        .iter()126        .filter(|&&v| v == FILL)127        .count()128}129130/// Returns the Euler characteristic of the cell's mesh, counting the backdrop only on request.131pub fn euler(cell: &Cell6d, include_grid: bool) -> i64 {132    census(cell, include_grid).euler133}134135/// Tallies only the filled triangles, leaving the voids and the backdrop out of the mesh.136pub fn fills_only(cell: &Cell6d) -> Census {137    let inner = &cell.cell;138    let start = cell.start as i64;139    let (height, width) = (inner.height(), inner.width());140    let mut fills = 0;141    let mut mesh = Mesh::default();142    for y in 0..height {143        for x in 0..width {144            if inner.types().get(&[y, x]) != FILL {145                continue;146            }147            fills += 1;148            mesh.add(x as i64, y as i64, start);149        }150    }151    mesh.tally(fills, fills, 0, 0)152}153154#[cfg(test)]155mod tests {156    use super::*;157    use crate::six::geometry::blank;158    use crate::six::{Orientation, Projection};159    use crate::two::Cell2d;160    fn hex(radius: usize) -> Cell6d {161        Cell6d::new(162            blank(radius, Orientation::Horizontal, FILL, VOID),163            Projection::Cut,164            Orientation::Horizontal,165            0,166        )167    }168    #[test]169    fn blank_hexagon_tallies_are_pinned() {170        let expected = [171            (1, 6, 6, 0, 7, 12, 6, 1),172            (2, 28, 24, 4, 22, 49, 14, 1),173            (3, 66, 54, 12, 45, 110, 22, 1),174        ];175        for (radius, triangles, fills, voids, vertices, edges, boundary, euler) in expected {176            let c = census(&hex(radius), false);177            assert_eq!(c.triangles, triangles, "r={radius}");178            assert_eq!(c.fills, fills);179            assert_eq!(c.voids, voids);180            assert_eq!(c.vertices, vertices);181            assert_eq!(c.edges, edges);182            assert_eq!(c.boundary_edges, boundary);183            assert_eq!(c.euler, euler);184        }185    }186    #[test]187    fn single_triangle() {188        let mut t = mrlycore::Tensor::new(vec![1, 2]);189        t.set(&[0, 0], FILL);190        t.set(&[0, 1], GRID);191        let c = census(192            &Cell6d::new(Cell2d::new(t), Projection::Cut, Orientation::Horizontal, 0),193            false,194        );195        assert_eq!(c.triangles, 1);196        assert_eq!(c.vertices, 3);197        assert_eq!(c.edges, 3);198        assert_eq!(c.boundary_edges, 3);199        assert_eq!(c.euler, 1);200    }201    #[test]202    fn fills_only_ignores_the_voids() {203        let solid = hex(2);204        let whole = census(&solid, false);205        let filled = fills_only(&solid);206        assert_eq!(fills(&solid), whole.fills);207        assert_eq!(euler(&solid, false), whole.euler);208        assert_eq!(filled.triangles, whole.fills);209        assert_eq!(filled.voids, 0);210        assert!(filled.edges < whole.edges);211        assert_eq!(fills_only(&solid.clone().anti()).triangles, whole.voids);212    }213}214215#[cfg(test)]216mod theorems {217    use super::*;218    use crate::formulas::six as formulas;219    use crate::six::geometry::cut;220    use crate::six::graph::slice_core_graph;221    use crate::three;222223    fn solid_slice(number: usize) -> Cell6d {224        cut(&three::ones(number, 1).unwrap()).unwrap()225    }226227    fn readings(rec: &Census) -> [i64; 5] {228        [229            rec.triangles as i64,230            rec.boundary_edges as i64,231            rec.edges as i64,232            rec.interior_edges as i64,233            rec.vertices as i64,234        ]235    }236237    fn lagrange(seed: [i64; 3], k: i64) -> i64 {238        seed[0] * (k - 2) * (k - 3) / 2 - seed[1] * (k - 1) * (k - 3)239            + seed[2] * (k - 1) * (k - 2) / 2240    }241242    #[test]243    fn frame_is_family_invariant() {244        for number in [3, 5, 7, 9, 11] {245            let reference = census(&solid_slice(number), false);246            assert_eq!(reference.euler, 1, "n={number}");247            let families = [248                three::carpet(number, 1).unwrap(),249                three::net(number, 1).unwrap(),250                three::ztree(number, 1).unwrap(),251                three::void(number, 1).unwrap(),252            ];253            for family in families {254                let rec = census(&cut(&family).unwrap(), false);255                assert_eq!(rec.triangles, reference.triangles, "n={number}");256                assert_eq!(rec.vertices, reference.vertices, "n={number}");257                assert_eq!(rec.edges, reference.edges, "n={number}");258                assert_eq!(rec.boundary_edges, reference.boundary_edges, "n={number}");259                assert_eq!(rec.euler, reference.euler, "n={number}");260                assert_eq!(rec.fills + rec.voids, reference.fills, "n={number}");261            }262        }263    }264265    #[test]266    fn the_five_closed_forms_match_the_census_to_eight() {267        for index in 1..9usize {268            let number = 2 * index - 1;269            let rec = census(&solid_slice(number), false);270            let k = index as i64;271            assert_eq!(readings(&rec)[0], 24 * k * k - 24 * k + 6, "k={k}");272            assert_eq!(readings(&rec)[1], 12 * k - 6, "k={k}");273            assert_eq!(readings(&rec)[2], 36 * k * k - 30 * k + 6, "k={k}");274            assert_eq!(readings(&rec)[3], 36 * k * k - 42 * k + 12, "k={k}");275            assert_eq!(readings(&rec)[4], 12 * k * k - 6 * k + 1, "k={k}");276            assert_eq!(rec.euler, 1, "k={k}");277            let closed = [278                formulas::solid_slice_triangles(number).unwrap(),279                formulas::solid_slice_boundary(number).unwrap(),280                formulas::solid_slice_edges(number).unwrap(),281                formulas::solid_slice_interior(number).unwrap(),282                formulas::solid_slice_vertices(number).unwrap(),283            ];284            for (got, want) in readings(&rec).iter().zip(closed) {285                assert_eq!(*got as u128, want, "k={k}");286            }287        }288        let three = census(&solid_slice(3), false);289        assert_eq!(readings(&three), [54, 18, 90, 72, 37]);290    }291292    #[test]293    fn a_blind_quadratic_fit_reproduces_the_wider_slices() {294        let seed: Vec<[i64; 5]> = (1..4)295            .map(|k| readings(&census(&solid_slice(2 * k - 1), false)))296            .collect();297        for k in 4..11i64 {298            let rec = census(&solid_slice(2 * k as usize - 1), false);299            assert_eq!(rec.euler, 1, "k={k}");300            let got = readings(&rec);301            for j in 0..5 {302                let fitted = lagrange([seed[0][j], seed[1][j], seed[2][j]], k);303                assert_eq!(got[j], fitted, "k={k} count={j}");304            }305        }306    }307308    #[test]309    fn fresh_builds_at_the_wide_sides_hold_the_forms() {310        for k in [12usize, 16, 20] {311            let number = 2 * k - 1;312            let rec = census(&solid_slice(number), false);313            assert_eq!(314                readings(&rec).map(|v| v as u128),315                [316                    formulas::solid_slice_triangles(number).unwrap(),317                    formulas::solid_slice_boundary(number).unwrap(),318                    formulas::solid_slice_edges(number).unwrap(),319                    formulas::solid_slice_interior(number).unwrap(),320                    formulas::solid_slice_vertices(number).unwrap(),321                ],322                "k={k}"323            );324            assert_eq!(rec.euler, 1, "k={k}");325            if number == 39 {326                assert_eq!(readings(&rec), [9126, 234, 13806, 13572, 4681]);327            }328        }329    }330331    #[test]332    fn the_fill_adjacency_counts_the_sub_mesh_interior_edges() {333        let mut meshes = 0;334        for level in 1..5usize {335            for design in [None, Some(23u128), Some(232), Some(3), Some(129)] {336                let cell = match design {337                    None => three::ones(3, level).unwrap(),338                    Some(code) => three::create(code, 3, level, 2).unwrap(),339                };340                let slice = cut(&cell).unwrap();341                let sub = fills_only(&slice);342                let core = slice_core_graph(&slice).unwrap();343                assert_eq!(344                    core.nodes.len(),345                    sub.triangles,346                    "design={design:?} l={level}"347                );348                assert_eq!(349                    core.branches.len(),350                    sub.edges - sub.boundary_edges,351                    "design={design:?} l={level}"352                );353                if design == Some(23) && level == 4 {354                    assert_eq!(355                        (sub.edges, sub.boundary_edges, core.branches.len()),356                        (28188, 6642, 21546)357                    );358                }359                meshes += 1;360            }361        }362        assert_eq!(meshes, 20);363    }364365    #[test]366    fn the_lemma_needs_the_sub_mesh_and_not_the_hexagon() {367        let slice = cut(&three::carpet(3, 3).unwrap()).unwrap();368        let sub = fills_only(&slice);369        let whole = census(&slice, false);370        assert_eq!(slice_core_graph(&slice).unwrap().branches.len(), 2880);371        assert_eq!(sub.edges - sub.boundary_edges, 2880);372        assert_eq!(whole.edges - whole.boundary_edges, 6480);373    }374}