sheaf_coherence/
coherence.rs1use crate::error::SheafError;
2use crate::laplacian::SheafLaplacian;
3use crate::sheaf::CellularSheaf;
4
5#[derive(Debug, Clone, serde::Serialize, serde::Deserialize)]
11pub struct CoherenceMeasure {
12 pub alignment: f64,
14 pub disagreement: Vec<f64>,
16 pub dominant_mode: Vec<f64>,
18}
19
20impl CoherenceMeasure {
21 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 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 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 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 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 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 pub fn is_aligned(&self, threshold: f64) -> bool {
81 self.alignment >= threshold
82 }
83
84 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 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}