Skip to main content

sheaf_coherence/
coherence.rs

1use crate::error::SheafError;
2use crate::laplacian::SheafLaplacian;
3use crate::sheaf::CellularSheaf;
4
5/// Measures how well a belief assignment aligns across the sheaf.
6///
7/// Alignment = 1 - ||L_F x|| / ||x|| (higher = more coherent).
8/// Per-edge disagreement = ||F_{ij} x_i - F_{ij} x_j|| where F_{ij} is the restriction map.
9/// Dominant mode = eigenvector of the smallest nonzero eigenvalue (most coherent non-trivial pattern).
10#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
11pub struct CoherenceMeasure {
12    /// Overall alignment score: 1 - ||L_F x|| / ||x||.
13    pub alignment: f64,
14    /// Per-edge disagreement values.
15    pub disagreement: Vec<f64>,
16    /// Eigenvector of smallest nonzero eigenvalue.
17    pub dominant_mode: Vec<f64>,
18}
19
20impl CoherenceMeasure {
21    /// Compute coherence measures for a flat belief vector.
22    pub fn from_flat(sheaf: &CellularSheaf, x: &[f64], max_iter: usize, tol: f64) -> Result<Self, SheafError> {
23        sheaf.validate()?;
24        let total_dim: usize = sheaf.stalk_dims.iter().sum();
25        if x.len() != total_dim {
26            return Err(SheafError::BeliefDimensionMismatch {
27                agent: "global".into(),
28                expected: total_dim,
29                got: x.len(),
30            });
31        }
32
33        let lap = SheafLaplacian::from_sheaf(sheaf)?;
34
35        // Alignment
36        let x_norm = norm(x);
37        let lx_norm = lap.residual_norm(x);
38        let alignment = if x_norm > 1e-15 {
39            1.0 - lx_norm / x_norm
40        } else {
41            1.0
42        };
43
44        // Per-edge disagreement
45        let mut offsets = Vec::with_capacity(sheaf.node_count());
46        let mut off = 0;
47        for &d in &sheaf.stalk_dims {
48            offsets.push(off);
49            off += d;
50        }
51
52        let mut disagreement = Vec::new();
53        for (src, tgt, f_map) in &sheaf.restriction_maps {
54            let x_src = &x[offsets[*src]..offsets[*src] + sheaf.stalk_dims[*src]];
55            let x_tgt = &x[offsets[*tgt]..offsets[*tgt] + sheaf.stalk_dims[*tgt]];
56            // disagreement = ||F x_src - F x_tgt|| (for undirected, both directions)
57            let f_src = mat_vec(f_map, x_src);
58            let f_tgt = mat_vec(f_map, x_tgt);
59            let diff: Vec<f64> = f_src.iter().zip(&f_tgt).map(|(a, b)| a - b).collect();
60            disagreement.push(norm(&diff));
61        }
62
63        // Dominant mode: eigenvector of smallest eigenvalue
64        let (_, dominant_mode) = lap.smallest_eigenvalue(max_iter, tol);
65
66        Ok(Self {
67            alignment: alignment.clamp(0.0, 1.0),
68            disagreement,
69            dominant_mode,
70        })
71    }
72
73    /// Compute coherence from per-node values.
74    pub fn from_values(sheaf: &CellularSheaf, values: &[Vec<f64>], max_iter: usize, tol: f64) -> Result<Self, SheafError> {
75        let flat: Vec<f64> = values.iter().flatten().copied().collect();
76        Self::from_flat(sheaf, &flat, max_iter, tol)
77    }
78
79    /// True if alignment is above a threshold.
80    pub fn is_aligned(&self, threshold: f64) -> bool {
81        self.alignment >= threshold
82    }
83
84    /// Average per-edge disagreement.
85    pub fn avg_disagreement(&self) -> f64 {
86        if self.disagreement.is_empty() {
87            return 0.0;
88        }
89        self.disagreement.iter().sum::<f64>() / self.disagreement.len() as f64
90    }
91
92    /// Maximum per-edge disagreement.
93    pub fn max_disagreement(&self) -> f64 {
94        self.disagreement.iter().cloned().fold(0.0_f64, f64::max)
95    }
96}
97
98fn mat_vec(m: &[Vec<f64>], v: &[f64]) -> Vec<f64> {
99    m.iter()
100        .map(|row| row.iter().zip(v.iter()).map(|(a, b)| a * b).sum())
101        .collect()
102}
103
104fn norm(v: &[f64]) -> f64 {
105    v.iter().map(|x| x * x).sum::<f64>().sqrt()
106}