surface.rs

1use crate::bang::factory;2use crate::bang::universe::Code;3use crate::formulas::counting::{fill_from_corners, positions};4use mrlycore::errors::{value_error, Result};5use mrlycore::tensor::Tensor;6use mrlynum::census::exposed;7use std::collections::HashSet;89fn strides(shape: &[usize]) -> Vec<usize> {10    (0..shape.len())11        .map(|axis| shape[axis + 1..].iter().product())12        .collect()13}1415fn occupancy(tile: &Tensor) -> u128 {16    tile.bytes().iter().filter(|&&v| v != 0).count() as u12817}1819/// Counts, per axis, the adjacent filled pairs and the cross positions whose two end cells are both filled.20///21/// A level deeper, each adjacent pair buries one face per spanning position of the block, and the22/// spanning positions of the block multiply level by level, so the exposure closes.23pub fn pairs(tile: &Tensor) -> Vec<(u128, u128)> {24    let shape = &tile.shape;25    let bytes = tile.bytes();26    let strides = strides(shape);27    (0..shape.len())28        .map(|axis| {29            let (stride, side) = (strides[axis], shape[axis]);30            let (mut adjacent, mut spanning) = (0u128, 0u128);31            for (flat, &cell) in bytes.iter().enumerate() {32                if cell == 0 {33                    continue;34                }35                let position = flat / stride % side;36                if position == 0 && bytes[flat + (side - 1) * stride] != 0 {37                    spanning += 1;38                }39                if position + 1 < side && bytes[flat + stride] != 0 {40                    adjacent += 1;41                }42            }43            (adjacent, spanning)44        })45        .collect()46}4748/// The counts the exposure recurrence runs on: the filled cells and exposed faces of the tile, and per axis its adjacent pairs and spanning positions.49///50/// With `occ` filled cells, `V(1)` exposed faces and per axis `P` adjacent pairs and `S` spanning51/// positions, `V(L + 1) = occ V(L) - 2 sum P S^L`: the perimeter in the plane, the surface in space.52#[derive(Clone, Debug, PartialEq, Eq)]53pub struct Exposure {54    /// The filled cells of the tile.55    pub occupancy: u128,56    /// The exposed faces of the tile.57    pub exposed: u128,58    /// Per axis, the adjacent filled pairs and the spanning positions.59    pub axes: Vec<(u128, u128)>,60}6162impl Exposure {63    /// Reads the counts off a rendered tile.64    pub fn of_tile(tile: &Tensor) -> Exposure {65        Exposure {66            occupancy: occupancy(tile),67            exposed: exposed(tile),68            axes: pairs(tile),69        }70    }71    /// Folds the counts from the filled residue corners at a side number, without rendering the tile.72    ///73    /// An adjacent pair sits at positions `i, i + 1` whose residues are `r, r + 1 mod q`, and a74    /// spanning position pairs residue 0 with the residue of `n - 1`.75    pub fn from_corners(76        filled: &[Vec<u8>],77        number: usize,78        dimension: usize,79        base: usize,80    ) -> Exposure {81        let set: HashSet<&Vec<u8>> = filled.iter().collect();82        let across = |corner: &[u8], axis: usize| -> u128 {83            (0..dimension)84                .filter(|&b| b != axis)85                .map(|b| positions(corner[b] as usize, number, base))86                .product()87        };88        let occupancy = fill_from_corners(filled, number, dimension, 1, base);89        let axes: Vec<(u128, u128)> = (0..dimension)90            .map(|axis| {91                let (mut adjacent, mut spanning) = (0u128, 0u128);92                for corner in filled {93                    let mut next = corner.clone();94                    next[axis] = ((corner[axis] as usize + 1) % base) as u8;95                    if set.contains(&next) {96                        adjacent +=97                            positions(corner[axis] as usize, number.saturating_sub(1), base)98                                * across(corner, axis);99                    }100                    if corner[axis] == 0 && number > 0 {101                        let mut far = corner.clone();102                        far[axis] = ((number - 1) % base) as u8;103                        if set.contains(&far) {104                            spanning += across(corner, axis);105                        }106                    }107                }108                (adjacent, spanning)109            })110            .collect();111        let buried: u128 = axes.iter().map(|&(adjacent, _)| adjacent).sum();112        Exposure {113            occupancy,114            exposed: 2 * dimension as u128 * occupancy - 2 * buried,115            axes,116        }117    }118    /// Returns the exposed faces of the level-fold Kronecker power, or none past a u128.119    pub fn at(&self, level: u32) -> Option<u128> {120        let mut value = self.exposed;121        for step in 1..level {122            let buried = self123                .axes124                .iter()125                .try_fold(0u128, |sum, &(adjacent, spanning)| {126                    sum.checked_add(adjacent.checked_mul(spanning.checked_pow(step)?)?)127                })?;128            value = self129                .occupancy130                .checked_mul(value)?131                .checked_sub(buried.checked_mul(2)?)?;132        }133        Some(value)134    }135    /// Returns the coefficients `c` of the recurrence `a(L) = c[0] a(L-1) + c[1] a(L-2) + ...` the exposure obeys.136    ///137    /// The roots are `occ` and the distinct nonzero spanning counts, `occ` doubled where a138    /// spanning count equals it.139    pub fn recurrence(&self) -> Vec<i128> {140        let mut roots: Vec<u128> = self141            .axes142            .iter()143            .map(|&(_, spanning)| spanning)144            .filter(|&spanning| spanning != 0)145            .collect();146        roots.sort_unstable();147        roots.dedup();148        roots.insert(0, self.occupancy);149        let mut poly: Vec<i128> = vec![1];150        for root in roots {151            let mut next = vec![0i128; poly.len() + 1];152            for (power, &coefficient) in poly.iter().enumerate() {153                next[power] += coefficient;154                next[power + 1] -= root as i128 * coefficient;155            }156            poly = next;157        }158        poly[1..].iter().map(|&coefficient| -coefficient).collect()159    }160}161162/// Returns the exposed face count of the tile's level-fold Kronecker power in closed form, or none past a u128.163///164/// ```165/// let carpet = mrlymath::bang::factory::create(7, 3, 2, 2, 1).unwrap();166/// let perimeter: Vec<u128> = (1..5).map(|level| mrlymath::formulas::exposure_of_tile(&carpet, level).unwrap()).collect();167/// assert_eq!(perimeter, [16, 80, 496, 3536]);168/// ```169pub fn exposure_of_tile(tile: &Tensor, level: u32) -> Option<u128> {170    Exposure::of_tile(tile).at(level)171}172173/// Returns the coefficients of the recurrence the tile's exposure obeys.174///175/// ```176/// let sponge = mrlymath::bang::factory::create(23, 3, 3, 2, 1).unwrap();177/// assert_eq!(mrlymath::formulas::exposure_recurrence(&sponge), [28, -160]);178/// ```179pub fn exposure_recurrence(tile: &Tensor) -> Vec<i128> {180    Exposure::of_tile(tile).recurrence()181}182183/// Returns the exposed face count of the code's fractal in any dimension at the given level, folded from its corners, or an error past a u128.184pub fn exposure(185    code: Code,186    number: usize,187    dimension: usize,188    level: u32,189    base: usize,190) -> Result<u128> {191    let filled = factory::code_to_corners(code, dimension, base)?;192    match Exposure::from_corners(&filled, number, dimension, base).at(level) {193        Some(value) => Ok(value),194        None => value_error("the exposure passes a hundred and twenty-eight bits."),195    }196}197198/// Returns the exposed face count of the code's 3D fractal at the given level.199pub fn surface(code: Code, number: usize, level: u32, base: usize) -> Result<u128> {200    exposure(code, number, 3, level, base)201}202203#[cfg(test)]204mod tests {205    use super::*;206    use mrlycore::atoms;207    #[test]208    fn prediction_matches_census_on_every_cube_code() {209        for code in 0..256u128 {210            for level in 1..4u32 {211                let direct = factory::create(code, 3, 3, 2, level as usize).unwrap();212                assert_eq!(213                    surface(code, 3, level, 2).unwrap(),214                    exposed(&direct),215                    "code={code} l={level}"216                );217            }218        }219    }220    #[test]221    fn prediction_matches_census_in_the_plane_and_beyond() {222        for code in 0..16u128 {223            for number in [2usize, 3, 4, 5] {224                for level in 1..4u32 {225                    let direct = factory::create(code, number, 2, 2, level as usize).unwrap();226                    assert_eq!(227                        exposure(code, number, 2, level, 2).unwrap(),228                        exposed(&direct),229                        "code={code} n={number} l={level}"230                    );231                }232            }233        }234        for code in [1u128, 23, 255, 4369, 65535, 32767] {235            for level in 1..3u32 {236                let direct = factory::create(code, 3, 4, 2, level as usize).unwrap();237                assert_eq!(exposure(code, 3, 4, level, 2).unwrap(), exposed(&direct));238            }239        }240        for code in [7u128, 100, 511] {241            let direct = factory::create(code, 3, 2, 3, 3).unwrap();242            assert_eq!(exposure(code, 3, 2, 3, 3).unwrap(), exposed(&direct));243        }244    }245    #[test]246    fn the_recurrence_holds_on_every_cube_code() {247        for code in 0..256u128 {248            let tile = factory::create(code, 3, 3, 2, 1).unwrap();249            let rule = exposure_recurrence(&tile);250            let terms: Vec<i128> = (1..8u32)251                .map(|level| exposure_of_tile(&tile, level).unwrap() as i128)252                .collect();253            for at in rule.len()..terms.len() {254                let predicted: i128 = rule255                    .iter()256                    .enumerate()257                    .map(|(back, &c)| c * terms[at - back - 1])258                    .sum();259                assert_eq!(predicted, terms[at], "code={code} at={at} rule={rule:?}");260            }261        }262    }263    #[test]264    fn the_corners_fold_what_the_tile_shows() {265        for code in 0..256u128 {266            for number in [1usize, 2, 3, 4, 5, 7] {267                let filled = factory::code_to_corners(code, 3, 2).unwrap();268                let tile = factory::create(code, number, 3, 2, 1).unwrap();269                assert_eq!(270                    Exposure::from_corners(&filled, number, 3, 2),271                    Exposure::of_tile(&tile),272                    "code={code} n={number}"273                );274            }275        }276        for (code, dimension, base) in [277            (7u128, 2usize, 3usize),278            (100, 2, 3),279            (511, 2, 3),280            (4369, 4, 2),281            (32767, 4, 2),282            (1, 1, 2),283            (2, 1, 3),284        ] {285            for number in [2usize, 3, 4, 5, 6, 9] {286                let filled = factory::code_to_corners(code, dimension, base).unwrap();287                let tile = factory::create(code, number, dimension, base, 1).unwrap();288                assert_eq!(289                    Exposure::from_corners(&filled, number, dimension, base),290                    Exposure::of_tile(&tile),291                    "code={code} d={dimension} q={base} n={number}"292                );293            }294        }295    }296    #[test]297    fn the_classics_close() {298        let sponge: Vec<u128> = (1..4).map(|l| surface(23, 3, l, 2).unwrap()).collect();299        assert_eq!(sponge, [72, 1056, 18048]);300        let carpet: Vec<u128> = (1..5).map(|l| exposure(7, 3, 2, l, 2).unwrap()).collect();301        assert_eq!(carpet, [16, 80, 496, 3536]);302        assert_eq!(303            exposure_recurrence(&factory::create(7, 3, 2, 2, 1).unwrap()),304            [11, -24]305        );306        assert_eq!(exposure_of_tile(&atoms::ones_3d(2), 3), Some(384));307        assert_eq!(exposure_of_tile(&atoms::ones_3d(1), 5), Some(6));308        assert!(exposure(23, 3, 3, 120, 2).is_err());309    }310}