topology.rs

1use super::census::{corners, edges_of, fills_only};2use super::graph::slice_core_graph;3use super::models::Cell6d;4use super::{FILL, VOID};5use mrlycore::errors::{value_error, Result};6use mrlynum::graph::largest_component;7use mrlynum::graph::models::Network;8use mrlynum::spectrum::{laplacian_spectrum, spectral_exponent as slope};9use std::collections::BTreeMap;1011type Point = (i64, i64);12type Edge = (Point, Point);1314fn pieces(network: &Network) -> Vec<usize> {15    let n = network.nodes.len();16    let adjacency = network.adjacency();17    let mut seen = vec![false; n];18    let mut sizes = Vec::new();19    for start in 0..n {20        if seen[start] {21            continue;22        }23        let mut size = 0;24        let mut stack = vec![start];25        seen[start] = true;26        while let Some(current) = stack.pop() {27            size += 1;28            for &neighbor in &adjacency[&current] {29                if !seen[neighbor] {30                    seen[neighbor] = true;31                    stack.push(neighbor);32                }33            }34        }35        sizes.push(size);36    }37    sizes38}3940/// Counts the connected pieces of the fill, triangles joined across shared edges.41///42/// ```43/// let slice = mrlymath::six::cut(&mrlymath::three::carpet(5, 1).unwrap()).unwrap();44/// assert_eq!(mrlymath::six::topology::components(&slice).unwrap(), 7);45/// ```46pub fn components(cell: &Cell6d) -> Result<usize> {47    Ok(pieces(&slice_core_graph(cell)?).len())48}4950/// Returns the triangle count of the fill's largest connected piece.51pub fn giant(cell: &Cell6d) -> Result<usize> {52    Ok(pieces(&slice_core_graph(cell)?)53        .into_iter()54        .max()55        .unwrap_or(0))56}5758/// Returns the largest connected piece of the filled-triangle network as a network of its own.59pub fn giant_network(cell: &Cell6d) -> Result<Network> {60    Ok(largest_component(&slice_core_graph(cell)?))61}6263/// Reads the spectral dimension of the giant piece: twice the low-window log-log slope of the normalised Laplacian's integrated density of states.64pub fn spectral_exponent(cell: &Cell6d, window: f64) -> Result<f64> {65    let spectrum = laplacian_spectrum(&giant_network(cell)?, true)?;66    match slope(&spectrum, window) {67        Some(value) => Ok(value),68        None => value_error("The giant piece is too small to fit an exponent."),69    }70}7172/// Counts the holes of the fill, its piece count less the Euler number of the filled sub-mesh.73///74/// ```75/// let slice = mrlymath::six::cut(&mrlymath::three::carpet(3, 1).unwrap()).unwrap();76/// assert_eq!(mrlymath::six::topology::holes(&slice).unwrap(), 1);77/// ```78pub fn holes(cell: &Cell6d) -> Result<usize> {79    let count = components(cell)? as i64;80    Ok((count - fills_only(cell).euler).max(0) as usize)81}8283/// Counts the void regions the rim never reaches, the second route to the hole count.84pub fn rim_holes(cell: &Cell6d) -> Result<usize> {85    let inner = &cell.cell;86    let start = cell.start as i64;87    let (height, width) = (inner.height(), inner.width());88    let mut sites = Vec::new();89    let mut mesh: BTreeMap<Edge, usize> = BTreeMap::new();90    for y in 0..height {91        for x in 0..width {92            let v = inner.types().get(&[y, x]);93            if v != FILL && v != VOID {94                continue;95            }96            for edge in edges_of(&corners(x as i64, y as i64, start)) {97                *mesh.entry(edge).or_insert(0) += 1;98            }99            if v == VOID {100                sites.push((x as i64, y as i64));101            }102        }103    }104    let mut owners: BTreeMap<Edge, Vec<usize>> = BTreeMap::new();105    for (index, &(x, y)) in sites.iter().enumerate() {106        for edge in edges_of(&corners(x, y, start)) {107            owners.entry(edge).or_default().push(index);108        }109    }110    let mut adjacency: Vec<Vec<usize>> = vec![Vec::new(); sites.len()];111    let mut on_rim = vec![false; sites.len()];112    for (edge, shared) in &owners {113        if mesh[edge] == 1 {114            for &index in shared {115                on_rim[index] = true;116            }117        }118        if shared.len() == 2 {119            adjacency[shared[0]].push(shared[1]);120            adjacency[shared[1]].push(shared[0]);121        }122    }123    let mut seen = vec![false; sites.len()];124    let mut enclosed = 0;125    for start in 0..sites.len() {126        if seen[start] {127            continue;128        }129        let mut open = false;130        let mut stack = vec![start];131        seen[start] = true;132        while let Some(current) = stack.pop() {133            open |= on_rim[current];134            for &neighbor in &adjacency[current] {135                if !seen[neighbor] {136                    seen[neighbor] = true;137                    stack.push(neighbor);138                }139            }140        }141        if !open {142            enclosed += 1;143        }144    }145    Ok(enclosed)146}147148#[cfg(test)]149mod theorems {150    use super::*;151    use crate::bang::universe::orbit;152    use crate::formulas::six::centered_hexagonal;153    use crate::six::census::census;154    use crate::six::geometry::cut;155    use crate::six::GRID;156    use crate::three::{self, Cell3d};157    use mrlynum::graph::census::components as network_components;158159    fn slice(code: u128, number: usize, level: usize) -> Cell6d {160        cut(&three::create(code, number, level, 2).unwrap()).unwrap()161    }162163    fn section_area(cell: &Cell3d) -> u128 {164        let grid = cell.types();165        let side = grid.shape[0];166        let mut total = 0;167        for i in 0..side {168            for j in 0..side {169                for l in 0..side {170                    if grid.get(&[i, j, l]) == 0 {171                        continue;172                    }173                    let layer = 3 * side as i64 - 2 * (i + j + l) as i64;174                    total += match layer {175                        1 | 5 => 1,176                        2 | 4 => 4,177                        3 => 6,178                        _ => 0,179                    };180                }181            }182        }183        total184    }185186    #[test]187    fn the_four_families_fill_the_slice_and_name_their_classes() {188        let expected = [189            (23u128, 42usize, 23u128),190            (232, 12, 23),191            (3, 18, 3),192            (129, 12, 24),193        ];194        for (code, fill, class) in expected {195            let cut = slice(code, 3, 1);196            assert_eq!(census(&cut, false).fills, fill, "code={code}");197            assert_eq!(*orbit(code, 3).iter().next().unwrap(), class, "code={code}");198        }199    }200201    #[test]202    fn carpet_and_net_partition_the_hexagon_triangle_by_triangle() {203        let pairs = [204            (3usize, 42u128, 12u128),205            (5, 72, 78),206            (7, 204, 90),207            (9, 210, 276),208            (11, 486, 240),209        ];210        for number in (1..32).step_by(2) {211            let carpet = slice(23, number, 1);212            let net = slice(232, number, 1);213            let (left, right) = (carpet.cell.types(), net.cell.types());214            assert_eq!(left.shape, right.shape, "n={number}");215            let mut filled = 0;216            for (index, &value) in left.bytes().iter().enumerate() {217                let other = right.bytes()[index];218                if value == GRID || other == GRID {219                    assert_eq!(value, other, "n={number}");220                    continue;221                }222                assert_ne!(value, other, "n={number}");223                if value == FILL || other == FILL {224                    filled += 1;225                }226            }227            let (a, b) = (228                census(&carpet, false).fills as u128,229                census(&net, false).fills as u128,230            );231            assert_eq!(filled as u128, 6 * (number as u128).pow(2), "n={number}");232            assert_eq!(a + b, 6 * (number as u128).pow(2), "n={number}");233            if let Some(&(_, want_a, want_b)) = pairs.iter().find(|p| p.0 == number) {234                assert_eq!((a, b), (want_a, want_b), "n={number}");235            }236            if number == 31 {237                assert_eq!((a, b, a + b), (3696, 2070, 5766));238            }239        }240    }241242    #[test]243    fn the_layer_weighted_area_is_a_second_route_to_the_fill() {244        for number in 1..17usize {245            let carpet = three::create(23, number, 1, 2).unwrap();246            let net = three::create(232, number, 1, 2).unwrap();247            let whole = 6 * (number as u128).pow(2);248            assert_eq!(249                section_area(&carpet) + section_area(&net),250                whole,251                "n={number}"252            );253            if number.is_multiple_of(2) {254                assert_eq!(section_area(&carpet), whole / 2, "n={number}");255                continue;256            }257            for code in [23u128, 232, 3, 129] {258                let cell = three::create(code, number, 1, 2).unwrap();259                assert_eq!(260                    section_area(&cell),261                    census(&cut(&cell).unwrap(), false).fills as u128,262                    "code={code} n={number}"263                );264            }265        }266    }267268    #[test]269    fn the_carpet_slice_counts_its_pieces_and_holes_two_ways() {270        let pieces = [1, 1, 7, 1, 19, 1, 37, 1, 61, 1, 91, 1, 127, 1];271        let punctures = [0, 1, 0, 7, 0, 19, 0, 37, 0, 61, 0, 91, 0, 127];272        for (index, (&want_pieces, &want_holes)) in pieces.iter().zip(punctures.iter()).enumerate()273        {274            let k = index + 1;275            let cut = slice(23, 2 * k - 1, 1);276            assert_eq!(components(&cut).unwrap(), want_pieces, "k={k}");277            assert_eq!(278                network_components(&slice_core_graph(&cut).unwrap()),279                want_pieces,280                "k={k}"281            );282            assert_eq!(holes(&cut).unwrap(), want_holes, "k={k}");283            assert_eq!(rim_holes(&cut).unwrap(), want_holes, "k={k}");284            let law = centered_hexagonal(k.div_ceil(2)) as usize;285            assert_eq!(286                if k.is_multiple_of(2) {287                    want_holes288                } else {289                    want_pieces290                },291                law,292                "k={k}"293            );294        }295    }296297    #[test]298    fn the_other_families_puncture_in_opposite_phase() {299        for k in 1..11usize {300            for code in [3u128, 129] {301                let cut = slice(code, 2 * k - 1, 1);302                assert_eq!(holes(&cut).unwrap(), 0, "code={code} k={k}");303                assert_eq!(rim_holes(&cut).unwrap(), 0, "code={code} k={k}");304            }305        }306        for (k, want) in [(3usize, 1usize), (5, 7), (7, 19), (9, 37)] {307            let cut = slice(232, 2 * k - 1, 1);308            assert_eq!(holes(&cut).unwrap(), want, "k={k}");309            assert_eq!(rim_holes(&cut).unwrap(), want, "k={k}");310        }311    }312313    #[test]314    fn the_carpet_slice_percolates_at_base_three() {315        for (level, triangles) in [(1usize, 42usize), (2, 306), (3, 2250), (4, 16578)] {316            let cut = slice(23, 3, level);317            assert_eq!(census(&cut, false).fills, triangles, "l={level}");318            assert_eq!(components(&cut).unwrap(), 1, "l={level}");319            assert_eq!(giant(&cut).unwrap(), triangles, "l={level}");320            if level == 4 {321                let core = slice_core_graph(&cut).unwrap();322                assert_eq!((core.nodes.len(), core.branches.len()), (16578, 21546));323            }324        }325    }326327    #[test]328    fn the_other_slices_shatter_or_never_grow_at_base_three() {329        for level in 1..5usize {330            let net = slice(232, 3, level);331            assert_eq!(census(&net, false).fills, 12, "l={level}");332            assert_eq!(components(&net).unwrap(), 1, "l={level}");333        }334        let mut tree = 0;335        let mut anti = 0;336        for level in 1..4usize {337            let cut = slice(3, 3, level);338            assert_eq!(giant(&cut).unwrap(), 8, "l={level}");339            assert!(components(&cut).unwrap() > tree, "l={level}");340            tree = components(&cut).unwrap();341            let cut = slice(129, 3, level);342            assert_eq!(giant(&cut).unwrap(), 6, "l={level}");343            assert!(components(&cut).unwrap() > anti, "l={level}");344            anti = components(&cut).unwrap();345        }346        assert_eq!((tree, anti), (40, 55));347    }348349    #[test]350    fn the_carpet_slice_fails_to_percolate_at_base_five() {351        let cut = slice(23, 5, 2);352        assert_eq!(census(&cut, false).fills, 1164);353        assert_eq!(components(&cut).unwrap(), 20);354        assert_eq!(giant(&cut).unwrap(), 192);355    }356}357358#[cfg(test)]359mod spectra {360    use super::*;361    use crate::six::geometry::cut;362    use crate::three;363    use mrlynum::spectrum::laplacian_spectrum;364365    fn slice(code: u128, number: usize, level: usize) -> Cell6d {366        cut(&three::create(code, number, level, 2).unwrap()).unwrap()367    }368369    fn reading(code: u128, number: usize, level: usize) -> (usize, usize, f64) {370        let cell = slice(code, number, level);371        let piece = giant_network(&cell).unwrap();372        let spectrum = laplacian_spectrum(&piece, true).unwrap();373        let zeros = spectrum.iter().filter(|v| **v < 1e-10).count();374        (375            piece.nodes.len(),376            zeros,377            spectral_exponent(&cell, 0.1).unwrap(),378        )379    }380381    #[test]382    fn the_slice_exponent_climbs_with_the_size_it_is_read_at() {383        let rows = [384            (23u128, 3usize, 1usize, 42usize, 0.91),385            (23, 3, 2, 306, 1.25),386            (255, 9, 1, 486, 1.61),387            (23, 5, 2, 192, 1.05),388        ];389        for (code, number, level, nodes, want) in rows {390            let (got_nodes, zeros, exponent) = reading(code, number, level);391            assert_eq!(got_nodes, nodes, "code={code} n={number} l={level}");392            assert_eq!(zeros, 1, "code={code} n={number} l={level}");393            assert!(394                (exponent - want).abs() < 0.005,395                "code={code} n={number} l={level} {exponent}"396            );397        }398    }399400    #[test]401    #[ignore = "two thousand nodes each; run it in release"]402    fn the_deep_slices_hold_their_spectral_exponents() {403        for (code, number, level, nodes, want) in [404            (23u128, 3usize, 3usize, 2250usize, 1.44),405            (255, 19, 1, 2166, 1.79),406        ] {407            let (got_nodes, zeros, exponent) = reading(code, number, level);408            assert_eq!(got_nodes, nodes, "code={code} n={number}");409            assert_eq!(zeros, 1, "code={code} n={number}");410            assert!(411                (exponent - want).abs() < 0.005,412                "code={code} n={number} {exponent}"413            );414        }415    }416}