geographdb_core/algorithms/

percolation.rs

//! Spatiotemporal percolation analysis for 4D graphs.
//!
//! Measures the phase transition where graph connectivity shifts from fragmented
//! to connected as the spatial query radius grows. Sweeping both radius and time
//! produces a 2D phase diagram — the percolation frontier in (radius, time) space.

use crate::algorithms::four_d::{
    reachable_4d, GraphNode4D, SpatialRegion, TemporalWindow, TraversalContext4D,
};
use glam::Vec3;

/// One measurement point in the percolation sweep.
#[derive(Debug, Clone)]
pub struct PercolationPoint {
    /// Query sphere radius at this measurement.
    pub radius: f32,
    /// Time slice centre used for the temporal window.
    pub time: u64,
    /// Nodes active under (radius, time) constraints.
    pub active: usize,
    /// Nodes reachable from the centre node under the same constraints.
    pub reachable: usize,
    /// `reachable / active`, or 0.0 when `active == 0`.
    pub fraction: f32,
}

/// Sweep a range of radii and time windows over a 4D node set.
///
/// For each (radius, time) pair, counts how many nodes are reachable from
/// `center_id` relative to the total active population. The spatial query is a
/// sphere of the given radius centred at `sphere_center`.
///
/// `time_half_width` is added/subtracted around each time value to form the
/// temporal window `[time - half, time + half]`.
pub fn percolation_sweep(
    nodes: &[GraphNode4D],
    center_id: u64,
    sphere_center: Vec3,
    radii: &[f32],
    time_slices: &[u64],
    time_half_width: u64,
) -> Vec<PercolationPoint> {
    let mut out = Vec::with_capacity(radii.len() * time_slices.len());

    for &time in time_slices {
        let t_start = time.saturating_sub(time_half_width);
        let t_end = time.saturating_add(time_half_width);
        let time_window = TemporalWindow {
            start: t_start,
            end: t_end,
        };

        for &radius in radii {
            let ctx = TraversalContext4D {
                spatial_region: Some(SpatialRegion::Sphere {
                    center: sphere_center,
                    radius,
                }),
                time_window: Some(time_window),
                ..TraversalContext4D::default()
            };

            let active = count_active(nodes, &ctx);
            let reachable = reachable_4d(nodes, center_id, &ctx).len();
            let fraction = if active == 0 {
                0.0
            } else {
                reachable as f32 / active as f32
            };

            out.push(PercolationPoint {
                radius,
                time,
                active,
                reachable,
                fraction,
            });
        }
    }

    out
}

/// Find the critical radius r* at a given time slice.
///
/// Returns the smallest radius in `points` (for the given `time`) where
/// `fraction >= threshold` AND `active >= min_active`. The `min_active` guard
/// prevents the trivial detection where a single active node yields fraction=1.0
/// at the smallest non-empty radius.
///
/// Returns `None` when no such point exists.
pub fn find_critical_radius(
    points: &[PercolationPoint],
    time: u64,
    threshold: f32,
    min_active: usize,
) -> Option<f32> {
    points
        .iter()
        .filter(|p| p.time == time && p.active >= min_active && p.fraction >= threshold)
        .map(|p| p.radius)
        .reduce(f32::min)
}

// ── helpers (not public, used by tests and the example) ─────────────────────

/// Count nodes that satisfy the traversal context filters without running BFS.
fn count_active(nodes: &[GraphNode4D], ctx: &TraversalContext4D) -> usize {
    nodes
        .iter()
        .filter(|n| {
            let time_ok = ctx
                .time_window
                .map(|w| w.overlaps(n.begin_ts, n.end_ts))
                .unwrap_or(true);
            let space_ok = match ctx.spatial_region {
                Some(SpatialRegion::Sphere { center, radius }) => {
                    n.position().distance(center) <= radius
                }
                Some(SpatialRegion::Aabb { min, max }) => {
                    let p = n.position();
                    p.x >= min.x
                        && p.x <= max.x
                        && p.y >= min.y
                        && p.y <= max.y
                        && p.z >= min.z
                        && p.z <= max.z
                }
                None => true,
            };
            time_ok && space_ok
        })
        .count()
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::algorithms::four_d::{GraphNode4D, GraphProperties, TemporalEdge};

    fn node(id: u64, x: f32, y: f32, z: f32, successors: Vec<TemporalEdge>) -> GraphNode4D {
        GraphNode4D {
            id,
            x,
            y,
            z,
            begin_ts: 0,
            end_ts: 100,
            properties: GraphProperties::new(),
            successors,
        }
    }

    fn edge(dst: u64) -> TemporalEdge {
        TemporalEdge {
            dst,
            weight: 1.0,
            begin_ts: 0,
            end_ts: 100,
        }
    }

    /// Build a line graph: 0 → 1 → 2 → … → (n-1), spaced 1.0 apart on the x-axis.
    fn line_graph(n: u64) -> Vec<GraphNode4D> {
        (0..n)
            .map(|i| {
                let succs = if i + 1 < n { vec![edge(i + 1)] } else { vec![] };
                node(i, i as f32, 0.0, 0.0, succs)
            })
            .collect()
    }

    // ── percolation_sweep ────────────────────────────────────────────────────

    #[test]
    fn test_zero_radius_gives_zero_fraction() {
        // Sphere so small it contains no nodes → active = 0, fraction = 0.
        let nodes = line_graph(5);
        let pts = percolation_sweep(
            &nodes,
            0,
            Vec3::new(-10.0, 0.0, 0.0), // centre far away
            &[0.001],
            &[50],
            10,
        );
        assert_eq!(pts.len(), 1);
        assert_eq!(
            pts[0].active, 0,
            "no node should fall inside radius 0.001 far from origin"
        );
        assert_eq!(pts[0].fraction, 0.0);
    }

    #[test]
    fn test_huge_radius_gives_full_fraction() {
        // Sphere large enough to contain everything → fraction = 1.0.
        let nodes = line_graph(5);
        let pts = percolation_sweep(&nodes, 0, Vec3::ZERO, &[1000.0], &[50], 10);
        assert_eq!(pts.len(), 1);
        assert_eq!(pts[0].active, 5);
        assert_eq!(pts[0].reachable, 5);
        assert!((pts[0].fraction - 1.0).abs() < 1e-6);
    }

    #[test]
    fn test_fraction_monotone_with_radius() {
        // Fraction must be non-decreasing as radius grows (at a fixed time).
        let nodes = line_graph(10);
        let radii: Vec<f32> = (1..=20).map(|i| i as f32 * 0.5).collect();
        let pts = percolation_sweep(&nodes, 0, Vec3::ZERO, &radii, &[50], 10);

        let fracs: Vec<f32> = pts
            .iter()
            .filter(|p| p.time == 50)
            .map(|p| p.fraction)
            .collect();
        for w in fracs.windows(2) {
            assert!(
                w[1] >= w[0] - 1e-6,
                "fraction not monotone: {} then {}",
                w[0],
                w[1]
            );
        }
    }

    #[test]
    fn test_time_filter_reduces_active_nodes() {
        // Nodes only active in [0, 100]; a window outside that range → 0 active.
        let nodes = line_graph(5);
        let pts = percolation_sweep(
            &nodes,
            0,
            Vec3::ZERO,
            &[1000.0],
            &[500], // time = 500, window [490, 510] — outside node validity [0,100]
            10,
        );
        assert_eq!(pts[0].active, 0);
        assert_eq!(pts[0].fraction, 0.0);
    }

    #[test]
    fn test_critical_radius_found_in_range() {
        // With a line of 10 nodes spaced 1.0 apart, reachability should cross
        // 0.5 somewhere between radius 1.0 and radius 10.0.
        let nodes = line_graph(10);
        let radii: Vec<f32> = (1..=20).map(|i| i as f32 * 0.5).collect();
        let pts = percolation_sweep(&nodes, 0, Vec3::ZERO, &radii, &[50], 10);
        let r_star = find_critical_radius(&pts, 50, 0.5, 1);
        assert!(r_star.is_some(), "critical radius should be found");
        let r = r_star.unwrap();
        assert!(
            (0.5..=10.0).contains(&r),
            "r* = {} out of expected range",
            r
        );
    }

    #[test]
    fn test_find_critical_radius_none_when_never_crosses() {
        // Threshold 0.5 on an isolated graph (no edges) → fraction = 1/N → None.
        let isolated: Vec<GraphNode4D> = (0..5)
            .map(|i| node(i, i as f32, 0.0, 0.0, vec![]))
            .collect();
        let pts2 = percolation_sweep(&isolated, 0, Vec3::ZERO, &[1000.0], &[50], 10);
        // Node 0 can reach only itself → fraction = 1/5 = 0.2
        let r_star = find_critical_radius(&pts2, 50, 0.5, 1);
        assert!(
            r_star.is_none(),
            "threshold 0.5 should not be crossed when fraction = 0.2"
        );
    }

    // ── count_active (internal helper) ──────────────────────────────────────

    #[test]
    fn test_count_active_no_filter() {
        let nodes = line_graph(7);
        let ctx = TraversalContext4D::default();
        assert_eq!(count_active(&nodes, &ctx), 7);
    }

    #[test]
    fn test_count_active_sphere_filter() {
        let nodes = line_graph(10); // nodes at x = 0..9, y=z=0
        let ctx = TraversalContext4D {
            spatial_region: Some(SpatialRegion::Sphere {
                center: Vec3::ZERO,
                radius: 3.5,
            }),
            ..TraversalContext4D::default()
        };
        // Nodes 0,1,2,3 are within radius 3.5 of origin (distances 0,1,2,3)
        assert_eq!(count_active(&nodes, &ctx), 4);
    }
}