reach.rs

1use super::models::Cell3d;2use mrlycore::errors::{value_error, Result};3use mrlycore::rng::Rng;4use std::collections::HashSet;56const DIRS: [[isize; 3]; 6] = [7    [1, 0, 0],8    [-1, 0, 0],9    [0, 1, 0],10    [0, -1, 0],11    [0, 0, 1],12    [0, 0, -1],13];1415/// The tally of one diffusion run over a cube's wall faces.16#[derive(Clone, Debug, PartialEq)]17pub struct Reach {18    /// The count of wall faces between an inner void and a fill.19    pub faces: usize,20    /// The count of wall faces some walker reacted on.21    pub used: usize,22    /// The used share of the wall faces.23    pub utilisation: f64,24    /// The share of walkers that reacted before exiting or stalling.25    pub reacted: f64,26    /// The count of open entry sites on the first face.27    pub entries: usize,28}2930fn shape_of(cell: &Cell3d) -> [usize; 3] {31    let grid = cell.types();32    [grid.shape[0], grid.shape[1], grid.shape[2]]33}3435fn shift(site: [usize; 3], step: &[isize; 3], shape: &[usize; 3]) -> Option<[usize; 3]> {36    let mut next = [0usize; 3];37    for axis in 0..3 {38        let moved = site[axis] as isize + step[axis];39        if moved < 0 || moved >= shape[axis] as isize {40            return None;41        }42        next[axis] = moved as usize;43    }44    Some(next)45}4647/// Returns every wall face as a void site and the direction index of its filled neighbor.48pub fn walls(cell: &Cell3d) -> Vec<[usize; 4]> {49    let grid = cell.types();50    let shape = shape_of(cell);51    let mut out = Vec::new();52    for x in 0..shape[0] {53        for y in 0..shape[1] {54            for z in 0..shape[2] {55                if grid.get(&[x, y, z]) != 0 {56                    continue;57                }58                for (d, step) in DIRS.iter().enumerate() {59                    if let Some(site) = shift([x, y, z], step, &shape) {60                        if grid.get(&site) != 0 {61                            out.push([x, y, z, d]);62                        }63                    }64                }65            }66        }67    }68    out69}7071fn doors(cell: &Cell3d) -> Vec<[usize; 3]> {72    let grid = cell.types();73    let shape = shape_of(cell);74    let mut out = Vec::new();75    for y in 0..shape[1] {76        for z in 0..shape[2] {77            if grid.get(&[0, y, z]) == 0 {78                out.push([0, y, z]);79            }80        }81    }82    out83}8485/// Walks seeded reactant in from the open first face and tallies the wall faces it reaches.86///87/// Each walker enters at a random void site of the first face and steps at random. A step88/// back out through that face exits, any other step off the cube reflects, and a step into89/// a fill reacts with the given probability, marking the wall face, or reflects otherwise.90pub fn explore(91    cell: &Cell3d,92    p_react: f64,93    walkers: usize,94    max_steps: usize,95    seed: u64,96) -> Result<Reach> {97    if walkers == 0 {98        return value_error("the walk needs at least one walker.");99    }100    let grid = cell.types();101    let shape = shape_of(cell);102    let faces = walls(cell);103    if faces.is_empty() {104        return value_error("the design offers no wall faces.");105    }106    let entry = doors(cell);107    if entry.is_empty() {108        return value_error("the first face offers no entry site.");109    }110    let mut rng = Rng::new(seed);111    let mut hit: HashSet<[usize; 4]> = HashSet::new();112    let mut reacted = 0usize;113    for _ in 0..walkers {114        let mut site = entry[rng.below(entry.len())];115        for _ in 0..max_steps {116            let d = rng.below(6);117            let step = &DIRS[d];118            let Some(next) = shift(site, step, &shape) else {119                if site[0] as isize + step[0] < 0 {120                    break;121                }122                continue;123            };124            if grid.get(&next) != 0 {125                if rng.chance(p_react) {126                    hit.insert([site[0], site[1], site[2], d]);127                    reacted += 1;128                    break;129                }130                continue;131            }132            site = next;133        }134    }135    Ok(Reach {136        faces: faces.len(),137        used: hit.len(),138        utilisation: hit.len() as f64 / faces.len() as f64,139        reacted: reacted as f64 / walkers as f64,140        entries: entry.len(),141    })142}143144/// Sweeps the reaction probability and pairs each value with its utilisation, same seed each run.145pub fn sweep(146    cell: &Cell3d,147    probabilities: &[f64],148    walkers: usize,149    max_steps: usize,150    seed: u64,151) -> Result<Vec<(f64, f64)>> {152    probabilities153        .iter()154        .map(|&p| Ok((p, explore(cell, p, walkers, max_steps, seed)?.utilisation)))155        .collect()156}157158#[cfg(test)]159mod tests {160    use super::*;161    use crate::three::{carpet, census, net, ones, void};162163    fn hull(cell: &Cell3d) -> u128 {164        let grid = cell.types();165        let shape = shape_of(cell);166        let mut count = 0u128;167        for x in 0..shape[0] {168            for y in 0..shape[1] {169                for z in 0..shape[2] {170                    if grid.get(&[x, y, z]) == 0 {171                        continue;172                    }173                    for (site, size) in [(x, shape[0]), (y, shape[1]), (z, shape[2])] {174                        if site == 0 {175                            count += 1;176                        }177                        if site + 1 == size {178                            count += 1;179                        }180                    }181                }182            }183        }184        count185    }186187    #[test]188    fn walls_complete_the_censused_surface() {189        for cell in [190            carpet(3, 1).unwrap(),191            carpet(3, 2).unwrap(),192            net(3, 1).unwrap(),193            void(2, 1).unwrap(),194            void(2, 2).unwrap(),195        ] {196            let inner = walls(&cell).len() as u128;197            assert_eq!(inner + hull(&cell), census::surface(&cell));198        }199    }200201    #[test]202    fn a_solid_block_of_eight_shows_only_its_hull() {203        let cell = ones(2, 1).unwrap();204        assert_eq!(census::fills(&cell), 8);205        assert_eq!(census::surface(&cell), 24);206        assert!(walls(&cell).is_empty());207        assert!(explore(&cell, 0.5, 8, 8, 0).is_err());208    }209210    #[test]211    fn the_carpet_cube_offers_twenty_four_walls() {212        let cell = carpet(3, 1).unwrap();213        assert_eq!(walls(&cell).len(), 24);214        assert_eq!(census::surface(&cell), 72);215        assert_eq!(explore(&cell, 0.5, 16, 32, 0).unwrap().entries, 1);216    }217218    #[test]219    fn a_fixed_seed_repeats_exactly() {220        let cell = carpet(3, 2).unwrap();221        let a = explore(&cell, 0.5, 64, 128, 9).unwrap();222        let b = explore(&cell, 0.5, 64, 128, 9).unwrap();223        assert_eq!(a, b);224        assert!(a.used <= a.faces);225        assert!((0.0..=1.0).contains(&a.utilisation));226        assert!((0.0..=1.0).contains(&a.reacted));227    }228229    #[test]230    fn faster_reactions_screen_the_interior() {231        let cell = carpet(3, 1).unwrap();232        let curve = sweep(&cell, &[0.05, 0.5, 0.9], 800, 400, 0).unwrap();233        for pair in curve.windows(2) {234            assert!(pair[0].1 >= pair[1].1, "{curve:?}");235        }236        assert!(curve[0].1 > curve[2].1, "{curve:?}");237    }238}