facett-graphview 0.1.12

//! **Force-Directed Edge Bundling (FDEB)** — the answer to the *hairball* curse (#2).
//!
//! Straight edges in a dense directed graph smear into an unreadable hairball. FDEB
//! (Holten & van Wijk, 2009) turns each edge into a flexible wire subdivided into
//! control points, then lets **compatible** edges — ones travelling the same way,
//! at a similar scale and close together — **attract each other's control points**.
//! Edges going the same direction collapse into shared "information highways" that
//! split again near their endpoints, so the eye follows trunks instead of spaghetti.
//!
//! # Reused force pipeline
//! The bundling force is the **same boids/n-body shape** as
//! [`facett_core::render::gpu::particles`] and [`crate::gpu_layout`]: a per-point
//! **electrostatic attraction** toward the matching control points of compatible
//! edges (cohesion) balanced by a **spring** holding each wire's own points evenly
//! spaced (the rubberband). It is a pure, deterministic function of the input
//! positions + params (FC-7): fixed iteration order, no RNG.
//!
//! # What lands here vs the roadmap
//! This is the **CPU reference** FDEB: compatibility (angle · scale · position) is
//! precomputed once, then a fixed-subdivision force loop bundles the wires. The
//! cycle-refinement schedule (double the subdivisions each cycle) and the **GPU**
//! port (the control points are particles; compatibility is the neighbour kernel —
//! it drops straight onto the `particles.wgsl` storage-ping-pong) are documented in
//! `.nornir/elite-graph-engine.md`. The output is bundled polylines a host paints as
//! glowing splines.

/// A 2D point.
pub type P2 = [f32; 2];

/// FDEB tuning. [`Default`] bundles moderately on a unit-ish coordinate scale.
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct BundleParams {
    /// Interior control points per edge (the wire's flexibility). Endpoints are extra.
    pub subdivisions: usize,
    /// Force-application iterations.
    pub iterations: usize,
    /// Per-iteration step size (how far a control point moves per step).
    pub step: f32,
    /// Global spring stiffness holding a wire's own points evenly spaced.
    pub stiffness: f32,
    /// Compatibility floor — edge pairs below this don't attract (the precompute prune).
    pub compat_threshold: f32,
}

impl Default for BundleParams {
    fn default() -> Self {
        Self { subdivisions: 6, iterations: 60, step: 0.04, stiffness: 0.1, compat_threshold: 0.05 }
    }
}

/// One bundled edge: its endpoints + the subdivided control polyline (endpoints
/// included, so `path.len() == subdivisions + 2`). Feed `path` to a spline painter.
#[derive(Clone, Debug, PartialEq)]
pub struct BundledEdge {
    pub from: usize,
    pub to: usize,
    pub path: Vec<P2>,
}

fn sub(a: P2, b: P2) -> P2 {
    [a[0] - b[0], a[1] - b[1]]
}
fn len(a: P2) -> f32 {
    (a[0] * a[0] + a[1] * a[1]).sqrt()
}

/// Edge-pair **compatibility** `Ca · Cs · Cp` (angle, scale, position) — the classic
/// FDEB measure in `[0,1]`. High when two edges run parallel, at a similar length,
/// close together (i.e. should bundle).
fn compatibility(p0: P2, p1: P2, q0: P2, q1: P2) -> f32 {
    let pv = sub(p1, p0);
    let qv = sub(q1, q0);
    let lp = len(pv).max(1e-6);
    let lq = len(qv).max(1e-6);
    // angle: |cos θ| between the two edge vectors.
    let ca = ((pv[0] * qv[0] + pv[1] * qv[1]) / (lp * lq)).abs();
    // scale: penalise very different lengths.
    let lavg = (lp + lq) * 0.5;
    let cs = 2.0 / (lavg * (lp.min(lq)).recip() + (lp.max(lq)) / lavg);
    // position: penalise distant midpoints (relative to lavg).
    let pm = [(p0[0] + p1[0]) * 0.5, (p0[1] + p1[1]) * 0.5];
    let qm = [(q0[0] + q1[0]) * 0.5, (q0[1] + q1[1]) * 0.5];
    let cp = lavg / (lavg + len(sub(pm, qm)));
    (ca * cs * cp).clamp(0.0, 1.0)
}

/// **Bundle** the `edges` (index pairs into `nodes`) with FDEB. Returns one
/// [`BundledEdge`] per input edge whose `path` is the bundled control polyline.
/// Deterministic in `(nodes, edges, p)`. Self-loops / out-of-range edges become a
/// straight degenerate path (left unbundled).
#[must_use]
pub fn bundle(nodes: &[P2], edges: &[(usize, usize)], p: &BundleParams) -> Vec<BundledEdge> {
    let m = p.subdivisions.max(1);
    let n = edges.len();
    // ── initialise each wire as evenly-spaced points on its straight chord ──
    let mut paths: Vec<Vec<P2>> = Vec::with_capacity(n);
    let mut endpoints: Vec<(P2, P2)> = Vec::with_capacity(n);
    for &(a, b) in edges {
        let (pa, pb) = match (nodes.get(a), nodes.get(b)) {
            (Some(&pa), Some(&pb)) => (pa, pb),
            _ => ([0.0, 0.0], [0.0, 0.0]),
        };
        endpoints.push((pa, pb));
        let path = (0..=m + 1)
            .map(|i| {
                let t = i as f32 / (m + 1) as f32;
                [pa[0] + (pb[0] - pa[0]) * t, pa[1] + (pb[1] - pa[1]) * t]
            })
            .collect();
        paths.push(path);
    }

    // ── precompute compatibility (symmetric, pruned at the threshold) ──
    let mut compat: Vec<Vec<(usize, f32)>> = vec![Vec::new(); n];
    for i in 0..n {
        let (pi0, pi1) = endpoints[i];
        for j in (i + 1)..n {
            let (qj0, qj1) = endpoints[j];
            let c = compatibility(pi0, pi1, qj0, qj1);
            if c >= p.compat_threshold {
                compat[i].push((j, c));
                compat[j].push((i, c));
            }
        }
    }

    // ── force iterations: spring (own wire) + electrostatic (compatible wires) ──
    for _ in 0..p.iterations {
        let snapshot = paths.clone();
        for e in 0..n {
            let rest = len(sub(endpoints[e].1, endpoints[e].0)).max(1e-4);
            let kp = p.stiffness / (rest * (m + 1) as f32);
            for i in 1..=m {
                let cur = snapshot[e][i];
                // spring toward the midpoint of this wire's own neighbours.
                let nb = snapshot[e][i - 1];
                let nf = snapshot[e][i + 1];
                let mut force = [
                    kp * ((nb[0] - cur[0]) + (nf[0] - cur[0])),
                    kp * ((nb[1] - cur[1]) + (nf[1] - cur[1])),
                ];
                // electrostatic attraction to compatible wires' matching points.
                for &(q, c) in &compat[e] {
                    let d = sub(snapshot[q][i], cur);
                    let dist = len(d).max(1e-3);
                    let w = c / dist; // inverse-distance, compatibility-weighted
                    force[0] += d[0] * w;
                    force[1] += d[1] * w;
                }
                paths[e][i] = [cur[0] + p.step * force[0], cur[1] + p.step * force[1]];
            }
        }
    }

    edges
        .iter()
        .zip(paths)
        .map(|(&(from, to), path)| BundledEdge { from, to, path })
        .collect()
}

/// The straight-edge baseline polylines (same subdivision count, no bundling) — the
/// control comparison for the occlusion metric.
#[must_use]
pub fn straight(nodes: &[P2], edges: &[(usize, usize)], subdivisions: usize) -> Vec<BundledEdge> {
    let p = BundleParams { subdivisions, iterations: 0, ..BundleParams::default() };
    bundle(nodes, edges, &p)
}

/// **Occlusion / ink-coverage metric**: how many distinct cells of a `cell`-sized grid
/// the edge polylines touch (sampled densely along each path). Bundling concentrates
/// edges onto shared trunks → fewer distinct cells → less occlusion. Read as DATA: a
/// drop vs the straight baseline proves the bundle de-hairballs the view.
#[must_use]
pub fn occupancy(edges: &[BundledEdge], cell: f32) -> usize {
    use std::collections::BTreeSet;
    let cell = cell.max(1e-4);
    let mut cells: BTreeSet<(i32, i32)> = BTreeSet::new();
    for e in edges {
        for w in e.path.windows(2) {
            let segs = 8;
            for s in 0..=segs {
                let t = s as f32 / segs as f32;
                let x = w[0][0] + (w[1][0] - w[0][0]) * t;
                let y = w[0][1] + (w[1][1] - w[0][1]) * t;
                cells.insert(((x / cell).floor() as i32, (y / cell).floor() as i32));
            }
        }
    }
    cells.len()
}

/// Total polyline length over all edges (bundling bows wires, so this is a secondary
/// diagnostic — the primary win is the [`occupancy`] drop).
#[must_use]
pub fn total_length(edges: &[BundledEdge]) -> f32 {
    edges.iter().map(|e| e.path.windows(2).map(|w| len(sub(w[1], w[0]))).sum::<f32>()).sum()
}

#[cfg(test)]
mod tests {
    use super::*;

    /// A bipartite hairball: 6 nodes on the left (x=0) fully connected to 6 nodes on
    /// the right (x=10). Straight, the 36 near-parallel edges smear across the whole
    /// middle band; they're highly compatible (parallel, same scale, overlapping) so
    /// FDEB should pull them into a trunk.
    fn hairball() -> (Vec<P2>, Vec<(usize, usize)>) {
        let mut nodes = Vec::new();
        for k in 0..6 {
            nodes.push([0.0, k as f32]); // left column 0..5
        }
        for k in 0..6 {
            nodes.push([10.0, k as f32]); // right column 6..11
        }
        let mut edges = Vec::new();
        for l in 0..6 {
            for r in 6..12 {
                edges.push((l, r));
            }
        }
        (nodes, edges)
    }

    #[test]
    fn endpoints_stay_pinned() {
        let (nodes, edges) = hairball();
        let bundled = bundle(&nodes, &edges, &BundleParams::default());
        for (b, &(a, c)) in bundled.iter().zip(&edges) {
            assert!(len(sub(*b.path.first().unwrap(), nodes[a])) < 1e-4, "start pinned to node");
            assert!(len(sub(*b.path.last().unwrap(), nodes[c])) < 1e-4, "end pinned to node");
            assert_eq!(b.path.len(), BundleParams::default().subdivisions + 2);
        }
    }

    #[test]
    fn bundling_reduces_occlusion_vs_straight() {
        let (nodes, edges) = hairball();
        let p = BundleParams::default();
        let straight_edges = straight(&nodes, &edges, p.subdivisions);
        let bundled = bundle(&nodes, &edges, &p);

        let cell = 0.5;
        let occ_straight = occupancy(&straight_edges, cell);
        let occ_bundled = occupancy(&bundled, cell);
        // The hairball collapses onto shared trunks → markedly fewer occupied cells.
        assert!(
            (occ_bundled as f32) < (occ_straight as f32) * 0.8,
            "bundling cut occlusion by >20%: {occ_straight} → {occ_bundled} cells"
        );
    }

    #[test]
    fn bundling_is_deterministic() {
        let (nodes, edges) = hairball();
        let p = BundleParams::default();
        assert_eq!(bundle(&nodes, &edges, &p), bundle(&nodes, &edges, &p), "FDEB is reproducible (FC-7)");
    }

    #[test]
    fn compatible_parallel_edges_converge_in_the_middle() {
        // Two parallel horizontal edges one unit apart. Their mid control points
        // should be pulled closer together than the straight 1.0 separation.
        let nodes = vec![[0.0, 0.0], [10.0, 0.0], [0.0, 1.0], [10.0, 1.0]];
        let edges = vec![(0, 1), (2, 3)];
        let p = BundleParams { iterations: 120, ..BundleParams::default() };
        let b = bundle(&nodes, &edges, &p);
        let mid = b[0].path.len() / 2;
        let gap = len(sub(b[1].path[mid], b[0].path[mid]));
        assert!(gap < 0.9, "parallel wires bundled closer than their 1.0 endpoint gap (got {gap})");
    }

    #[test]
    fn perpendicular_edges_are_incompatible() {
        // A horizontal and a vertical edge crossing: low angle-compatibility ⇒ they
        // barely attract, so the horizontal wire stays ~straight.
        let nodes = vec![[0.0, 0.0], [10.0, 0.0], [5.0, -5.0], [5.0, 5.0]];
        let edges = vec![(0, 1), (2, 3)];
        let b = bundle(&nodes, &edges, &BundleParams::default());
        let mid = b[0].path.len() / 2;
        // The horizontal wire's midpoint barely leaves the y=0 chord.
        assert!(b[0].path[mid][1].abs() < 0.6, "perpendicular edge didn't bundle (y={})", b[0].path[mid][1]);
    }
}