surface.rs

1use crate::core::error::{value_error, Result};2use crate::core::tensor::Tensor;3use crate::math::bang::factory;4use crate::math::bang::Code;5use crate::math::counts::counting::{fill_from_corners, positions};6use serde::{Deserialize, Serialize};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.size() - tile.count(0)) 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 strides = strides(shape);26    (0..shape.len())27        .map(|axis| {28            let (stride, side) = (strides[axis], shape[axis]);29            let (mut adjacent, mut spanning) = (0u128, 0u128);30            for flat in 0..tile.size() {31                if tile.at(flat) == 0 {32                    continue;33                }34                let position = flat / stride % side;35                if position == 0 && tile.at(flat + (side - 1) * stride) != 0 {36                    spanning += 1;37                }38                if position + 1 < side && tile.at(flat + stride) != 0 {39                    adjacent += 1;40                }41            }42            (adjacent, spanning)43        })44        .collect()45}4647/// 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.48///49/// With `occ` filled cells, `V(1)` exposed faces and per axis `P` adjacent pairs and `S` spanning50/// positions, `V(L + 1) = occ V(L) - 2 sum P S^L`: the perimeter in the plane, the surface in space.51#[derive(Clone, Debug, PartialEq, Eq, Serialize, Deserialize)]52pub struct Exposure {53    /// The filled cells of the tile.54    pub occupancy: u128,55    /// The exposed faces of the tile.56    pub exposed: u128,57    /// Per axis, the adjacent filled pairs and the spanning positions.58    pub axes: Vec<(u128, u128)>,59}6061impl Exposure {62    /// Reads the counts off a rendered tile.63    pub fn of_tile(tile: &Tensor) -> Exposure {64        Exposure {65            occupancy: occupancy(tile),66            exposed: tile.exposed(),67            axes: pairs(tile),68        }69    }70    /// Folds the counts from the filled residue corners at a side number, without rendering the tile.71    ///72    /// An adjacent pair sits at positions `i, i + 1` whose residues are `r, r + 1 mod q`, and a73    /// spanning position pairs residue 0 with the residue of `n - 1`.74    pub fn from_corners(75        filled: &[Vec<u8>],76        number: usize,77        dimension: usize,78        base: usize,79    ) -> Exposure {80        let set: HashSet<&Vec<u8>> = filled.iter().collect();81        let across = |corner: &[u8], axis: usize| -> u128 {82            (0..dimension)83                .filter(|&b| b != axis)84                .map(|b| positions(corner[b] as usize, number, base))85                .product()86        };87        let occupancy = fill_from_corners(filled, number, dimension, 1, base);88        let axes: Vec<(u128, u128)> = (0..dimension)89            .map(|axis| {90                let (mut adjacent, mut spanning) = (0u128, 0u128);91                for corner in filled {92                    let mut next = corner.clone();93                    next[axis] = ((corner[axis] as usize + 1) % base) as u8;94                    if set.contains(&next) {95                        adjacent +=96                            positions(corner[axis] as usize, number.saturating_sub(1), base)97                                * across(corner, axis);98                    }99                    if corner[axis] == 0 && number > 0 {100                        let mut far = corner.clone();101                        far[axis] = ((number - 1) % base) as u8;102                        if set.contains(&far) {103                            spanning += across(corner, axis);104                        }105                    }106                }107                (adjacent, spanning)108            })109            .collect();110        let buried: u128 = axes.iter().map(|&(adjacent, _)| adjacent).sum();111        Exposure {112            occupancy,113            exposed: 2 * dimension as u128 * occupancy - 2 * buried,114            axes,115        }116    }117    /// Returns the exposed faces of the level-fold Kronecker power, or none past a u128.118    pub fn at(&self, level: u32) -> Option<u128> {119        let mut value = self.exposed;120        for step in 1..level {121            let buried = self122                .axes123                .iter()124                .try_fold(0u128, |sum, &(adjacent, spanning)| {125                    sum.checked_add(adjacent.checked_mul(spanning.checked_pow(step)?)?)126                })?;127            value = self128                .occupancy129                .checked_mul(value)?130                .checked_sub(buried.checked_mul(2)?)?;131        }132        Some(value)133    }134    /// Returns the coefficients `c` of the recurrence `a(L) = c[0] a(L-1) + c[1] a(L-2) + ...` the exposure obeys.135    ///136    /// The roots are `occ` and the distinct nonzero spanning counts, `occ` doubled where a137    /// spanning count equals it.138    pub fn recurrence(&self) -> Vec<i128> {139        let mut roots: Vec<u128> = self140            .axes141            .iter()142            .map(|&(_, spanning)| spanning)143            .filter(|&spanning| spanning != 0)144            .collect();145        roots.sort_unstable();146        roots.dedup();147        roots.insert(0, self.occupancy);148        let mut poly: Vec<i128> = vec![1];149        for root in roots {150            let mut next = vec![0i128; poly.len() + 1];151            for (power, &coefficient) in poly.iter().enumerate() {152                next[power] += coefficient;153                next[power + 1] -= root as i128 * coefficient;154            }155            poly = next;156        }157        poly[1..].iter().map(|&coefficient| -coefficient).collect()158    }159}160161/// Returns the exposed face count of the tile's level-fold Kronecker power in closed form, or none past a u128.162///163/// ```164/// use mrlyrs::math::bang::Code;165/// let carpet = mrlyrs::math::bang::factory::create(Code::from(7u64), 3, 2, 2, 1).unwrap();166/// let perimeter: Vec<u128> = (1..5).map(|level| mrlyrs::math::counts::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/// use mrlyrs::math::bang::Code;177/// let sponge = mrlyrs::math::bang::factory::create(Code::from(23u64), 3, 3, 2, 1).unwrap();178/// assert_eq!(mrlyrs::math::counts::exposure_recurrence(&sponge), [28, -160]);179/// ```180pub fn exposure_recurrence(tile: &Tensor) -> Vec<i128> {181    Exposure::of_tile(tile).recurrence()182}183184/// Returns the exposed face count of the code's fractal in any dimension at the given level, folded from its corners.185///186/// # Errors187///188/// Errors when the code is out of range, or the exposure passes a hundred and twenty-eight bits.189pub fn exposure(190    code: Code,191    number: usize,192    dimension: usize,193    level: u32,194    base: usize,195) -> Result<u128> {196    let filled = factory::code_to_corners(code, dimension, base)?;197    match Exposure::from_corners(&filled, number, dimension, base).at(level) {198        Some(value) => Ok(value),199        None => value_error("the exposure passes a hundred and twenty-eight bits."),200    }201}202203/// Returns the exposed face count of the code's 3D fractal at the given level.204///205/// # Errors206///207/// Errors when the code is out of range, or the exposure passes a hundred and twenty-eight bits.208pub fn surface(code: Code, number: usize, level: u32, base: usize) -> Result<u128> {209    exposure(code, number, 3, level, base)210}211212#[cfg(test)]213mod tests {214    use super::*;215    use crate::math::atoms;216    #[test]217    fn prediction_matches_census_on_every_cube_code() {218        for bits in 0..256u128 {219            let code = Code(bits);220            for level in 1..3u32 {221                let direct = factory::create(code, 3, 3, 2, level as usize).unwrap();222                assert_eq!(223                    surface(code, 3, level, 2).unwrap(),224                    direct.exposed(),225                    "code={code} l={level}"226                );227            }228        }229    }230    #[test]231    #[ignore = "two hundred and fifty-six level-three cubes, 2 s; run it in release"]232    fn prediction_matches_census_on_every_cube_code_at_the_third_level() {233        for bits in 0..256u128 {234            let code = Code(bits);235            let direct = factory::create(code, 3, 3, 2, 3).unwrap();236            assert_eq!(237                surface(code, 3, 3, 2).unwrap(),238                direct.exposed(),239                "code={code}"240            );241        }242    }243    #[test]244    fn prediction_matches_census_in_the_plane_and_beyond() {245        for bits in 0..16u128 {246            let code = Code(bits);247            for number in [2usize, 3, 4, 5] {248                for level in 1..4u32 {249                    let direct = factory::create(code, number, 2, 2, level as usize).unwrap();250                    assert_eq!(251                        exposure(code, number, 2, level, 2).unwrap(),252                        direct.exposed(),253                        "code={code} n={number} l={level}"254                    );255                }256            }257        }258        for code in [259            Code(1),260            Code(23),261            Code(255),262            Code(4369),263            Code(65535),264            Code(32767),265        ] {266            for level in 1..3u32 {267                let direct = factory::create(code, 3, 4, 2, level as usize).unwrap();268                assert_eq!(exposure(code, 3, 4, level, 2).unwrap(), direct.exposed());269            }270        }271        for code in [Code(7), Code(100), Code(511)] {272            let direct = factory::create(code, 3, 2, 3, 3).unwrap();273            assert_eq!(exposure(code, 3, 2, 3, 3).unwrap(), direct.exposed());274        }275    }276    #[test]277    fn the_recurrence_holds_on_every_cube_code() {278        for bits in 0..256u128 {279            let code = Code(bits);280            let tile = factory::create(code, 3, 3, 2, 1).unwrap();281            let rule = exposure_recurrence(&tile);282            let terms: Vec<i128> = (1..8u32)283                .map(|level| exposure_of_tile(&tile, level).unwrap() as i128)284                .collect();285            for at in rule.len()..terms.len() {286                let predicted: i128 = rule287                    .iter()288                    .enumerate()289                    .map(|(back, &c)| c * terms[at - back - 1])290                    .sum();291                assert_eq!(predicted, terms[at], "code={code} at={at} rule={rule:?}");292            }293        }294    }295    #[test]296    fn the_corners_fold_what_the_tile_shows() {297        for bits in 0..256u128 {298            let code = Code(bits);299            for number in [1usize, 2, 3, 4, 5, 7] {300                let filled = factory::code_to_corners(code, 3, 2).unwrap();301                let tile = factory::create(code, number, 3, 2, 1).unwrap();302                assert_eq!(303                    Exposure::from_corners(&filled, number, 3, 2),304                    Exposure::of_tile(&tile),305                    "code={code} n={number}"306                );307            }308        }309        for (code, dimension, base) in [310            (Code(7), 2usize, 3usize),311            (Code(100), 2, 3),312            (Code(511), 2, 3),313            (Code(4369), 4, 2),314            (Code(32767), 4, 2),315            (Code(1), 1, 2),316            (Code(2), 1, 3),317        ] {318            for number in [2usize, 3, 4, 5, 6, 9] {319                let filled = factory::code_to_corners(code, dimension, base).unwrap();320                let tile = factory::create(code, number, dimension, base, 1).unwrap();321                assert_eq!(322                    Exposure::from_corners(&filled, number, dimension, base),323                    Exposure::of_tile(&tile),324                    "code={code} d={dimension} q={base} n={number}"325                );326            }327        }328    }329    #[test]330    fn the_classics_close() {331        let sponge: Vec<u128> = (1..4)332            .map(|l| surface(Code(23), 3, l, 2).unwrap())333            .collect();334        assert_eq!(sponge, [72, 1056, 18048]);335        let carpet: Vec<u128> = (1..5)336            .map(|l| exposure(Code(7), 3, 2, l, 2).unwrap())337            .collect();338        assert_eq!(carpet, [16, 80, 496, 3536]);339        assert_eq!(340            exposure_recurrence(&factory::create(Code(7), 3, 2, 2, 1).unwrap()),341            [11, -24]342        );343        assert_eq!(exposure_of_tile(&atoms::ones_3d(2), 3), Some(384));344        assert_eq!(exposure_of_tile(&atoms::ones_3d(1), 5), Some(6));345        assert!(exposure(Code(23), 3, 3, 120, 2).is_err());346    }347}