census.rs

1use super::models::{Branch, Network, Node};2use mrlycore::logs;3use std::collections::HashSet;4use std::f64::consts::LN_2;56const RUNGS: u32 = 32;78/// Sums the straight-line lengths of every branch.9pub fn total_length(network: &Network) -> f64 {10    network11        .branches12        .iter()13        .map(|b| {14            let a = &network.nodes[b.parent].position;15            let c = &network.nodes[b.child].position;16            a.iter()17                .zip(c)18                .map(|(x, y)| (x - y) * (x - y))19                .sum::<f64>()20                .sqrt()21        })22        .sum()23}2425/// What a node's degree makes it.26#[derive(Clone, Copy, Debug, PartialEq)]27pub enum Role {28    /// A node no branch touches.29    Alone,30    /// A node exactly one branch touches.31    Tip,32    /// A node exactly two branches pass through.33    Through,34    /// A node three or more branches meet at.35    Junction,36}3738/// Tags every node by its degree, indexed like the node list.39///40/// ```41/// let mut net = mrlynum::graph::Network::new(2);42/// net.add_node(vec![0.0, 0.0]).unwrap();43/// net.add_node(vec![1.0, 0.0]).unwrap();44/// net.add_branch(0, 1, 1.0).unwrap();45/// assert_eq!(mrlynum::graph::roles(&net), vec![mrlynum::graph::Role::Tip; 2]);46/// ```47pub fn roles(network: &Network) -> Vec<Role> {48    network49        .degree()50        .iter()51        .map(|&d| match d {52            0 => Role::Alone,53            1 => Role::Tip,54            2 => Role::Through,55            _ => Role::Junction,56        })57        .collect()58}5960/// Counts the nodes of degree one.61pub fn tips(network: &Network) -> usize {62    roles(network).iter().filter(|&&r| r == Role::Tip).count()63}6465/// Counts the nodes of degree three or more.66pub fn junctions(network: &Network) -> usize {67    roles(network)68        .iter()69        .filter(|&&r| r == Role::Junction)70        .count()71}7273/// Counts the connected components of the network.74pub fn components(network: &Network) -> usize {75    let n = network.nodes.len();76    if n == 0 {77        return 0;78    }79    let adjacency = network.adjacency();80    let mut seen = vec![false; n];81    let mut count = 0;82    for start in 0..n {83        if seen[start] {84            continue;85        }86        count += 1;87        let mut stack = vec![start];88        seen[start] = true;89        while let Some(current) = stack.pop() {90            for &neighbor in &adjacency[&current] {91                if !seen[neighbor] {92                    seen[neighbor] = true;93                    stack.push(neighbor);94                }95            }96        }97    }98    count99}100101/// Extracts the largest connected piece as a network of its own, branches re-indexed.102///103/// Ties go to the piece whose lowest node index comes first. An empty network comes back empty.104///105/// ```106/// let mut net = mrlynum::graph::Network::new(1);107/// for i in 0..3 { net.add_node(vec![i as f64]).unwrap(); }108/// net.add_branch(0, 1, 1.0).unwrap();109/// assert_eq!(mrlynum::graph::largest_component(&net).nodes.len(), 2);110/// ```111pub fn largest_component(network: &Network) -> Network {112    let n = network.nodes.len();113    let adjacency = network.adjacency();114    let mut label = vec![usize::MAX; n];115    let mut sizes: Vec<usize> = Vec::new();116    for start in 0..n {117        if label[start] != usize::MAX {118            continue;119        }120        let piece = sizes.len();121        let mut size = 0;122        let mut stack = vec![start];123        label[start] = piece;124        while let Some(current) = stack.pop() {125            size += 1;126            for &neighbor in &adjacency[&current] {127                if label[neighbor] == usize::MAX {128                    label[neighbor] = piece;129                    stack.push(neighbor);130                }131            }132        }133        sizes.push(size);134    }135    let mut best = 0;136    for (piece, &size) in sizes.iter().enumerate() {137        if size > sizes[best] {138            best = piece;139        }140    }141    let mut giant = Network::new(network.dim);142    if sizes.is_empty() {143        return giant;144    }145    let mut index_of = vec![usize::MAX; n];146    for (old, node) in network.nodes.iter().enumerate() {147        if label[old] != best {148            continue;149        }150        let index = giant.nodes.len();151        index_of[old] = index;152        giant.nodes.push(Node {153            position: node.position.clone(),154            index,155        });156    }157    for branch in &network.branches {158        if label[branch.parent] != best {159            continue;160        }161        giant.branches.push(Branch {162            parent: index_of[branch.parent],163            child: index_of[branch.child],164            radius: branch.radius,165        });166    }167    giant168}169170/// Estimates the box-counting dimension of the node cloud over a ladder of halving boxes.171pub fn fractal_dimension(network: &Network) -> f64 {172    let positions: Vec<&Vec<f64>> = network.nodes.iter().map(|n| &n.position).collect();173    if positions.len() < 2 {174        return 0.0;175    }176    let dim = network.dim;177    let mins: Vec<f64> = (0..dim)178        .map(|a| positions.iter().map(|p| p[a]).fold(f64::MAX, f64::min))179        .collect();180    let maxs: Vec<f64> = (0..dim)181        .map(|a| positions.iter().map(|p| p[a]).fold(f64::MIN, f64::max))182        .collect();183    let extent = mins184        .iter()185        .zip(&maxs)186        .map(|(lo, hi)| hi - lo)187        .fold(0.0, f64::max);188    if extent == 0.0 {189        return 0.0;190    }191    let distinct: HashSet<Vec<u64>> = positions192        .iter()193        .map(|p| p.iter().map(|v| v.to_bits()).collect())194        .collect();195    let mut scales: Vec<f64> = Vec::new();196    let mut log_count: Vec<f64> = Vec::new();197    for k in 0..=RUNGS {198        let split = f64::from_bits(((1023 + k) as u64) << 52);199        let last = (1i64 << k) - 1;200        let mut boxes: HashSet<Vec<i64>> = HashSet::new();201        for p in &positions {202            let key: Vec<i64> = (0..dim)203                .map(|a| ((((p[a] - mins[a]) * split) / extent).floor() as i64).min(last))204                .collect();205            boxes.insert(key);206        }207        scales.push(k as f64);208        log_count.push(logs::ln(boxes.len() as f64));209        if boxes.len() == distinct.len() {210            break;211        }212    }213    let n = scales.len() as f64;214    let mean_x: f64 = scales.iter().sum::<f64>() / n;215    let mean_y: f64 = log_count.iter().sum::<f64>() / n;216    let cov: f64 = scales217        .iter()218        .zip(&log_count)219        .map(|(x, y)| (x - mean_x) * (y - mean_y))220        .sum();221    let var: f64 = scales.iter().map(|x| (x - mean_x) * (x - mean_x)).sum();222    cov / (var * LN_2)223}224225/// The measurements of one network.226#[derive(Clone, Debug, PartialEq)]227pub struct Census {228    /// The node count.229    pub nodes: usize,230    /// The branch count.231    pub branches: usize,232    /// The count of degree-one nodes.233    pub tips: usize,234    /// The count of nodes of degree three or more.235    pub junctions: usize,236    /// The connected component count.237    pub components: usize,238    /// The summed branch length.239    pub total_length: f64,240    /// The box-counting dimension estimate.241    pub fractal_dimension: f64,242}243244/// Takes the full census of a network.245///246/// ```247/// let mut net = mrlynum::graph::Network::new(2);248/// net.add_node(vec![0.0, 0.0]).unwrap();249/// net.add_node(vec![3.0, 4.0]).unwrap();250/// assert_eq!(mrlynum::graph::census(&net).components, 2);251/// ```252pub fn census(network: &Network) -> Census {253    Census {254        nodes: network.nodes.len(),255        branches: network.branches.len(),256        tips: tips(network),257        junctions: junctions(network),258        components: components(network),259        total_length: total_length(network),260        fractal_dimension: fractal_dimension(network),261    }262}263264#[cfg(test)]265mod tests {266    use super::*;267    use crate::graph::extract::core_graph;268    use mrlycore::atoms;269    #[test]270    fn carpet_census() {271        let network = core_graph(&atoms::carpet_2d(3)).unwrap();272        let c = census(&network);273        assert_eq!(c.nodes, 8);274        assert_eq!(c.branches, 8);275        assert_eq!(c.components, 1);276        assert!((c.total_length - 8.0).abs() < 1e-9);277    }278    #[test]279    fn carpet_census_holds_its_pinned_counts() {280        let network = core_graph(&atoms::carpet_2d(3).fractal(4)).unwrap();281        let c = census(&network);282        assert_eq!(c.nodes, 4096);283        assert_eq!(c.branches, 6424);284        assert_eq!(c.tips, 0);285        assert_eq!(c.junctions, 3596);286        assert_eq!(c.components, 1);287        assert!((c.total_length - 6424.0).abs() < 1e-9);288    }289    #[test]290    fn the_two_node_ladder_reads_exactly_one() {291        let mut net = Network::new(2);292        net.add_node(vec![0.0, 0.0]).unwrap();293        net.add_node(vec![3.0, 4.0]).unwrap();294        assert_eq!(fractal_dimension(&net), 1.0);295        net.add_node(vec![0.0, 0.0]).unwrap();296        assert_eq!(fractal_dimension(&net), 1.0);297    }298    #[test]299    fn the_carpet_dimension_holds_its_pinned_value() {300        let network = core_graph(&atoms::carpet_2d(3).fractal(4)).unwrap();301        let d = census(&network).fractal_dimension;302        assert!((0.0..=3.0).contains(&d), "dimension {d}");303        assert_eq!(d, 1.787589465914211);304    }305}