census.rs

1use super::models::{Branch, Network, Node};2use crate::core::error::Result;3use serde::{Deserialize, Serialize};4use std::collections::HashSet;5use std::f64::consts::LN_2;67const RUNGS: u32 = 32;89/// Sums the straight-line lengths of every branch.10pub fn total_length(network: &Network) -> f64 {11    network12        .branches13        .iter()14        .map(|b| {15            let a = &network.nodes[b.parent].position;16            let c = &network.nodes[b.child].position;17            a.iter()18                .zip(c)19                .map(|(x, y)| (x - y) * (x - y))20                .sum::<f64>()21                .sqrt()22        })23        .sum()24}2526/// What a node's degree makes it.27#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]28pub enum Role {29    /// A node no branch touches.30    Alone,31    /// A node exactly one branch touches.32    Tip,33    /// A node exactly two branches pass through.34    Through,35    /// A node three or more branches meet at.36    Junction,37}3839/// Tags every node by its degree, indexed like the node list.40///41/// ```42/// let mut net = mrlyrs::math::graph::Network::new(2);43/// net.add_node(vec![0.0, 0.0]).unwrap();44/// net.add_node(vec![1.0, 0.0]).unwrap();45/// net.add_branch(0, 1, 1.0).unwrap();46/// assert_eq!(mrlyrs::math::graph::roles(&net).unwrap(), vec![mrlyrs::math::graph::Role::Tip; 2]);47/// ```48///49/// # Errors50///51/// Errors when a branch names a node the network does not hold.52pub fn roles(network: &Network) -> Result<Vec<Role>> {53    Ok(network54        .degree()?55        .iter()56        .map(|&d| match d {57            0 => Role::Alone,58            1 => Role::Tip,59            2 => Role::Through,60            _ => Role::Junction,61        })62        .collect())63}6465/// Counts the nodes of degree one.66///67/// # Errors68///69/// Errors when a branch names a node the network does not hold.70pub fn tips(network: &Network) -> Result<usize> {71    Ok(roles(network)?.iter().filter(|&&r| r == Role::Tip).count())72}7374/// Counts the nodes of degree three or more.75///76/// # Errors77///78/// Errors when a branch names a node the network does not hold.79pub fn junctions(network: &Network) -> Result<usize> {80    Ok(roles(network)?81        .iter()82        .filter(|&&r| r == Role::Junction)83        .count())84}8586/// Counts the connected components of the network.87///88/// # Errors89///90/// Errors when a branch names a node the network does not hold.91pub fn components(network: &Network) -> Result<usize> {92    let n = network.nodes.len();93    if n == 0 {94        return Ok(0);95    }96    let adjacency = network.adjacency()?;97    let mut seen = vec![false; n];98    let mut count = 0;99    for start in 0..n {100        if seen[start] {101            continue;102        }103        count += 1;104        let mut stack = vec![start];105        seen[start] = true;106        while let Some(current) = stack.pop() {107            for &neighbor in &adjacency[&current] {108                if !seen[neighbor] {109                    seen[neighbor] = true;110                    stack.push(neighbor);111                }112            }113        }114    }115    Ok(count)116}117118/// Extracts the largest connected piece as a network of its own, branches re-indexed.119///120/// Ties go to the piece whose lowest node index comes first. An empty network comes back empty.121///122/// ```123/// let mut net = mrlyrs::math::graph::Network::new(1);124/// for i in 0..3 { net.add_node(vec![i as f64]).unwrap(); }125/// net.add_branch(0, 1, 1.0).unwrap();126/// assert_eq!(mrlyrs::math::graph::largest_component(&net).unwrap().nodes.len(), 2);127/// ```128///129/// # Errors130///131/// Errors when a branch names a node the network does not hold.132pub fn largest_component(network: &Network) -> Result<Network> {133    let n = network.nodes.len();134    let adjacency = network.adjacency()?;135    let mut label = vec![usize::MAX; n];136    let mut sizes: Vec<usize> = Vec::new();137    for start in 0..n {138        if label[start] != usize::MAX {139            continue;140        }141        let piece = sizes.len();142        let mut size = 0;143        let mut stack = vec![start];144        label[start] = piece;145        while let Some(current) = stack.pop() {146            size += 1;147            for &neighbor in &adjacency[&current] {148                if label[neighbor] == usize::MAX {149                    label[neighbor] = piece;150                    stack.push(neighbor);151                }152            }153        }154        sizes.push(size);155    }156    let mut best = 0;157    for (piece, &size) in sizes.iter().enumerate() {158        if size > sizes[best] {159            best = piece;160        }161    }162    let mut giant = Network::new(network.dim);163    if sizes.is_empty() {164        return Ok(giant);165    }166    let mut index_of = vec![usize::MAX; n];167    for (old, node) in network.nodes.iter().enumerate() {168        if label[old] != best {169            continue;170        }171        let index = giant.nodes.len();172        index_of[old] = index;173        giant.nodes.push(Node {174            position: node.position.clone(),175            index,176        });177    }178    for branch in &network.branches {179        if label[branch.parent] != best {180            continue;181        }182        giant.branches.push(Branch {183            parent: index_of[branch.parent],184            child: index_of[branch.child],185            radius: branch.radius,186        });187    }188    Ok(giant)189}190191/// Estimates the box-counting dimension of the node cloud over a ladder of halving boxes.192pub fn fractal_dimension(network: &Network) -> f64 {193    let positions: Vec<&Vec<f64>> = network.nodes.iter().map(|n| &n.position).collect();194    if positions.len() < 2 {195        return 0.0;196    }197    let dim = network.dim;198    let mins: Vec<f64> = (0..dim)199        .map(|a| positions.iter().map(|p| p[a]).fold(f64::MAX, f64::min))200        .collect();201    let maxs: Vec<f64> = (0..dim)202        .map(|a| positions.iter().map(|p| p[a]).fold(f64::MIN, f64::max))203        .collect();204    let extent = mins205        .iter()206        .zip(&maxs)207        .map(|(lo, hi)| hi - lo)208        .fold(0.0, f64::max);209    if extent == 0.0 {210        return 0.0;211    }212    let distinct: HashSet<Vec<u64>> = positions213        .iter()214        .map(|p| p.iter().map(|v| v.to_bits()).collect())215        .collect();216    let mut scales: Vec<f64> = Vec::new();217    let mut log_count: Vec<f64> = Vec::new();218    for k in 0..=RUNGS {219        let split = f64::from_bits(((1023 + k) as u64) << 52);220        let last = (1i64 << k) - 1;221        let mut boxes: HashSet<Vec<i64>> = HashSet::new();222        for p in &positions {223            let key: Vec<i64> = (0..dim)224                .map(|a| ((((p[a] - mins[a]) * split) / extent).floor() as i64).min(last))225                .collect();226            boxes.insert(key);227        }228        scales.push(k as f64);229        log_count.push(f64::ln(boxes.len() as f64));230        if boxes.len() == distinct.len() {231            break;232        }233    }234    let n = scales.len() as f64;235    let mean_x: f64 = scales.iter().sum::<f64>() / n;236    let mean_y: f64 = log_count.iter().sum::<f64>() / n;237    let cov: f64 = scales238        .iter()239        .zip(&log_count)240        .map(|(x, y)| (x - mean_x) * (y - mean_y))241        .sum();242    let var: f64 = scales.iter().map(|x| (x - mean_x) * (x - mean_x)).sum();243    cov / (var * LN_2)244}245246/// The measurements of one network.247#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]248pub struct Census {249    /// The node count.250    pub nodes: usize,251    /// The branch count.252    pub branches: usize,253    /// The count of degree-one nodes.254    pub tips: usize,255    /// The count of nodes of degree three or more.256    pub junctions: usize,257    /// The connected component count.258    pub components: usize,259    /// The summed branch length.260    pub total_length: f64,261    /// The box-counting dimension estimate.262    pub fractal_dimension: f64,263}264265/// Takes the full census of a network.266///267/// ```268/// let mut net = mrlyrs::math::graph::Network::new(2);269/// net.add_node(vec![0.0, 0.0]).unwrap();270/// net.add_node(vec![3.0, 4.0]).unwrap();271/// assert_eq!(mrlyrs::math::graph::census(&net).unwrap().components, 2);272/// ```273///274/// # Errors275///276/// Errors when a branch names a node the network does not hold.277pub fn census(network: &Network) -> Result<Census> {278    Ok(Census {279        nodes: network.nodes.len(),280        branches: network.branches.len(),281        tips: tips(network)?,282        junctions: junctions(network)?,283        components: components(network)?,284        total_length: total_length(network),285        fractal_dimension: fractal_dimension(network),286    })287}288289#[cfg(test)]290mod tests {291    use super::*;292    use crate::math::atoms;293    use crate::math::graph::extract::core_graph;294    #[test]295    fn carpet_census() {296        let network = core_graph(&atoms::carpet_2d(3)).unwrap();297        let c = census(&network).unwrap();298        assert_eq!(c.nodes, 8);299        assert_eq!(c.branches, 8);300        assert_eq!(c.components, 1);301        assert!((c.total_length - 8.0).abs() < 1e-9);302    }303    #[test]304    fn carpet_census_holds_its_pinned_counts() {305        let network = core_graph(&atoms::carpet_2d(3).fractal(4)).unwrap();306        let c = census(&network).unwrap();307        assert_eq!(c.nodes, 4096);308        assert_eq!(c.branches, 6424);309        assert_eq!(c.tips, 0);310        assert_eq!(c.junctions, 3596);311        assert_eq!(c.components, 1);312        assert!((c.total_length - 6424.0).abs() < 1e-9);313    }314    #[test]315    fn the_two_node_ladder_reads_exactly_one() {316        let mut net = Network::new(2);317        net.add_node(vec![0.0, 0.0]).unwrap();318        net.add_node(vec![3.0, 4.0]).unwrap();319        assert_eq!(fractal_dimension(&net), 1.0);320        net.add_node(vec![0.0, 0.0]).unwrap();321        assert_eq!(fractal_dimension(&net), 1.0);322    }323    #[test]324    fn the_carpet_dimension_holds_its_pinned_value() {325        let network = core_graph(&atoms::carpet_2d(3).fractal(4)).unwrap();326        let d = census(&network).unwrap().fractal_dimension;327        assert!((0.0..=3.0).contains(&d), "dimension {d}");328        assert!((d - 1.787589465914211).abs() < 1e-12, "dimension {d}");329    }330331    #[test]332    fn refuses_a_branch_past_the_nodes() {333        let mut net = Network::new(1);334        net.add_node(vec![0.0]).unwrap();335        net.branches.push(Branch {336            parent: 0,337            child: 4,338            radius: 1.0,339        });340        assert!(roles(&net).is_err());341        assert!(tips(&net).is_err());342        assert!(junctions(&net).is_err());343        assert!(components(&net).is_err());344        assert!(largest_component(&net).is_err());345        assert!(census(&net).is_err());346    }347}