spectrum.rs
15.2 kB · rust · 466 lines
1#![allow(clippy::needless_range_loop)]23use crate::graph::models::Network;4use mrlycore::errors::{value_error, Result};56const SWEEPS: usize = 60;7const ZERO: f64 = 1e-12;89// SOLVER1011fn tridiagonalise(a: &mut [Vec<f64>]) -> (Vec<f64>, Vec<f64>) {12 let n = a.len();13 let mut d = vec![0.0; n];14 let mut e = vec![0.0; n];15 for i in (1..n).rev() {16 let l = i - 1;17 let mut h = 0.0;18 let mut scale = 0.0;19 if l > 0 {20 for k in 0..=l {21 scale += a[i][k].abs();22 }23 }24 if l == 0 || scale == 0.0 {25 e[i] = a[i][l];26 continue;27 }28 for k in 0..=l {29 a[i][k] /= scale;30 h += a[i][k] * a[i][k];31 }32 let mut f = a[i][l];33 let g = if f >= 0.0 { -h.sqrt() } else { h.sqrt() };34 e[i] = scale * g;35 h -= f * g;36 a[i][l] = f - g;37 f = 0.0;38 for j in 0..=l {39 let mut sum = 0.0;40 for k in 0..=j {41 sum += a[j][k] * a[i][k];42 }43 for k in (j + 1)..=l {44 sum += a[k][j] * a[i][k];45 }46 e[j] = sum / h;47 f += e[j] * a[i][j];48 }49 let hh = f / (h + h);50 for j in 0..=l {51 let f = a[i][j];52 let g = e[j] - hh * f;53 e[j] = g;54 for k in 0..=j {55 a[j][k] -= f * e[k] + g * a[i][k];56 }57 }58 }59 e[0] = 0.0;60 for i in 0..n {61 d[i] = a[i][i];62 }63 (d, e)64}6566fn implicit_ql(d: &mut [f64], e: &mut [f64]) -> Result<()> {67 let n = d.len();68 for i in 1..n {69 e[i - 1] = e[i];70 }71 if n > 0 {72 e[n - 1] = 0.0;73 }74 for l in 0..n {75 let mut sweeps = 0;76 loop {77 let mut m = l;78 while m + 1 < n {79 let dd = d[m].abs() + d[m + 1].abs();80 if e[m].abs() <= f64::EPSILON * dd {81 break;82 }83 m += 1;84 }85 if m == l {86 break;87 }88 sweeps += 1;89 if sweeps > SWEEPS {90 return value_error("The QL iteration did not converge.");91 }92 let mut g = (d[l + 1] - d[l]) / (2.0 * e[l]);93 let mut r = g.hypot(1.0);94 g = d[m] - d[l] + e[l] / (g + if g >= 0.0 { r.abs() } else { -r.abs() });95 let mut s = 1.0;96 let mut c = 1.0;97 let mut p = 0.0;98 let mut underflow = false;99 let mut i = m;100 while i > l {101 i -= 1;102 let f = s * e[i];103 let b = c * e[i];104 r = f.hypot(g);105 e[i + 1] = r;106 if r == 0.0 {107 d[i + 1] -= p;108 e[m] = 0.0;109 underflow = true;110 break;111 }112 s = f / r;113 c = g / r;114 g = d[i + 1] - p;115 r = (d[i] - g) * s + 2.0 * c * b;116 p = s * r;117 d[i + 1] = g + p;118 g = c * r - b;119 }120 if underflow {121 continue;122 }123 d[l] -= p;124 e[l] = g;125 e[m] = 0.0;126 }127 }128 Ok(())129}130131/// Returns the eigenvalues of a dense real symmetric matrix in ascending order.132///133/// The matrix is reduced to tridiagonal form by Householder reflections and then134/// diagonalised by implicit QL with Wilkinson shifts. Errors on a matrix that is135/// not square or whose entries disagree across the diagonal by more than `1e-12`.136///137/// ```138/// let m = vec![vec![2.0, 1.0], vec![1.0, 2.0]];139/// let values = mrlynum::spectrum::symmetric_eigenvalues(&m).unwrap();140/// assert!((values[0] - 1.0).abs() < 1e-12 && (values[1] - 3.0).abs() < 1e-12);141/// ```142pub fn symmetric_eigenvalues(matrix: &[Vec<f64>]) -> Result<Vec<f64>> {143 let n = matrix.len();144 for row in matrix {145 if row.len() != n {146 return value_error("The matrix is not square.");147 }148 }149 for i in 0..n {150 for j in (i + 1)..n {151 if (matrix[i][j] - matrix[j][i]).abs() > ZERO {152 return value_error(format!("The matrix is not symmetric at {i}, {j}."));153 }154 }155 }156 let mut work: Vec<Vec<f64>> = matrix.to_vec();157 let (mut d, mut e) = tridiagonalise(&mut work);158 implicit_ql(&mut d, &mut e)?;159 d.sort_by(|a, b| a.partial_cmp(b).unwrap());160 Ok(d)161}162163// LAPLACIAN164165/// Builds the Laplacian of a network, the combinatorial `D - A` or the normalised `I - D^-1/2 A D^-1/2`.166///167/// Self-loops are ignored and the degrees are the adjacency row sums. Errors when a168/// node carries no branch and the normalised form is asked for.169pub fn laplacian(network: &Network, normalised: bool) -> Result<Vec<Vec<f64>>> {170 let n = network.nodes.len();171 let mut matrix = vec![vec![0.0; n]; n];172 for branch in &network.branches {173 if branch.parent == branch.child {174 continue;175 }176 matrix[branch.parent][branch.child] += 1.0;177 matrix[branch.child][branch.parent] += 1.0;178 }179 let degree: Vec<f64> = matrix.iter().map(|row| row.iter().sum()).collect();180 for (index, &d) in degree.iter().enumerate() {181 if normalised && d == 0.0 {182 return value_error(format!("Node {index} carries no branch."));183 }184 }185 for i in 0..n {186 for j in 0..n {187 matrix[i][j] = if normalised {188 -matrix[i][j] / (degree[i] * degree[j]).sqrt()189 } else {190 -matrix[i][j]191 };192 }193 matrix[i][i] = if normalised { 1.0 } else { degree[i] };194 }195 Ok(matrix)196}197198/// Returns the ascending Laplacian spectrum of a network, combinatorial or normalised.199///200/// ```201/// let mut net = mrlynum::graph::Network::new(1);202/// net.add_node(vec![0.0]).unwrap();203/// net.add_node(vec![1.0]).unwrap();204/// net.add_branch(0, 1, 1.0).unwrap();205/// let values = mrlynum::spectrum::laplacian_spectrum(&net, true).unwrap();206/// assert!(values[0].abs() < 1e-12 && (values[1] - 2.0).abs() < 1e-12);207/// ```208pub fn laplacian_spectrum(network: &Network, normalised: bool) -> Result<Vec<f64>> {209 symmetric_eigenvalues(&laplacian(network, normalised)?)210}211212// READINGS213214/// Groups eigenvalues into runs split by consecutive gaps above the tolerance, each run its mean and its size.215///216/// ```217/// let groups = mrlynum::spectrum::clusters(&[0.0, 1e-15, 2.0], 1e-9);218/// assert_eq!(groups.len(), 2);219/// assert_eq!(groups[0].1, 2);220/// ```221pub fn clusters(eigenvalues: &[f64], tolerance: f64) -> Vec<(f64, usize)> {222 let mut values = eigenvalues.to_vec();223 values.sort_by(|a, b| a.partial_cmp(b).unwrap());224 let mut groups = Vec::new();225 let mut start = 0;226 for index in 1..=values.len() {227 if index == values.len() || values[index] - values[index - 1] > tolerance {228 let run = &values[start..index];229 groups.push((run.iter().sum::<f64>() / run.len() as f64, run.len()));230 start = index;231 }232 }233 groups234}235236/// Counts the eigenvalues within the tolerance of a value.237pub fn multiplicity(eigenvalues: &[f64], value: f64, tolerance: f64) -> usize {238 eigenvalues239 .iter()240 .filter(|v| (**v - value).abs() <= tolerance)241 .count()242}243244/// Builds the integrated density of states as points, each an eigenvalue and its rank fraction.245///246/// The eigenvalues are clamped at zero and sorted, exactly equal neighbours collapse to one247/// point at the run's last index over the total count, and points at or below `1e-12` drop.248///249/// ```250/// let points = mrlynum::spectrum::spectral_points(&[0.0, 0.5, 0.5, 2.0]);251/// assert_eq!(points.len(), 2);252/// assert!((points[0].1 - 0.75).abs() < 1e-12);253/// ```254pub fn spectral_points(eigenvalues: &[f64]) -> Vec<(f64, f64)> {255 let mut values: Vec<f64> = eigenvalues.iter().map(|&v| v.max(0.0)).collect();256 values.sort_by(|a, b| a.partial_cmp(b).unwrap());257 let total = values.len();258 let mut points = Vec::new();259 let mut index = 0;260 while index < total {261 let mut last = index;262 while last + 1 < total && values[last + 1] == values[index] {263 last += 1;264 }265 if values[last] > ZERO {266 points.push((values[last], (last + 1) as f64 / total as f64));267 }268 index = last + 1;269 }270 points271}272273/// Fits the low window of the integrated density of states in log-log: the intercept, the slope and the fitted count.274///275/// The points come from `spectral_points` and the first `max(floor(window * total), 3)` of276/// them are fitted by least squares. Returns `None` when fewer than two points remain.277pub fn spectral_fit(eigenvalues: &[f64], window: f64) -> Option<(f64, f64, usize)> {278 let points = spectral_points(eigenvalues);279 let xs: Vec<f64> = points.iter().map(|p| p.0.ln()).collect();280 let ys: Vec<f64> = points.iter().map(|p| p.1.ln()).collect();281 let top = ((window * eigenvalues.len() as f64).floor() as usize)282 .max(3)283 .min(xs.len());284 if top < 2 {285 return None;286 }287 let count = top as f64;288 let mean_x: f64 = xs[..top].iter().sum::<f64>() / count;289 let mean_y: f64 = ys[..top].iter().sum::<f64>() / count;290 let covariance: f64 = xs[..top]291 .iter()292 .zip(&ys[..top])293 .map(|(x, y)| (x - mean_x) * (y - mean_y))294 .sum();295 let variance: f64 = xs[..top].iter().map(|x| (x - mean_x) * (x - mean_x)).sum();296 if variance == 0.0 {297 return None;298 }299 let slope = covariance / variance;300 Some((mean_y - slope * mean_x, slope, top))301}302303/// Reads the spectral exponent: twice the log-log slope of the integrated density of states over its low window.304pub fn spectral_exponent(eigenvalues: &[f64], window: f64) -> Option<f64> {305 spectral_fit(eigenvalues, window).map(|(_, slope, _)| 2.0 * slope)306}307308#[cfg(test)]309mod tests {310 use super::*;311 use std::f64::consts::PI;312313 fn path(n: usize) -> Network {314 let mut net = Network::new(1);315 for i in 0..n {316 net.add_node(vec![i as f64]).unwrap();317 }318 for i in 1..n {319 net.add_branch(i - 1, i, 1.0).unwrap();320 }321 net322 }323324 fn complete(n: usize) -> Network {325 let mut net = Network::new(1);326 for i in 0..n {327 net.add_node(vec![i as f64]).unwrap();328 }329 for i in 0..n {330 for j in (i + 1)..n {331 net.add_branch(i, j, 1.0).unwrap();332 }333 }334 net335 }336337 fn cycle(n: usize) -> Network {338 let mut net = path(n);339 net.add_branch(n - 1, 0, 1.0).unwrap();340 net341 }342343 #[test]344 fn the_path_laplacian_reads_its_closed_form() {345 for n in 2..12usize {346 let values = laplacian_spectrum(&path(n), false).unwrap();347 let mut want: Vec<f64> = (0..n)348 .map(|k| 2.0 - 2.0 * (PI * k as f64 / n as f64).cos())349 .collect();350 want.sort_by(|a, b| a.partial_cmp(b).unwrap());351 for (got, expected) in values.iter().zip(&want) {352 assert!((got - expected).abs() < 1e-10, "n={n} {got} {expected}");353 }354 }355 }356357 #[test]358 fn the_complete_laplacian_is_zero_once_and_n_the_rest() {359 for n in 2..10usize {360 let values = laplacian_spectrum(&complete(n), false).unwrap();361 assert!(values[0].abs() < 1e-10, "n={n}");362 for value in &values[1..] {363 assert!((value - n as f64).abs() < 1e-10, "n={n} {value}");364 }365 }366 }367368 #[test]369 fn the_cycle_normalised_laplacian_reads_its_closed_form() {370 for n in 3..12usize {371 let values = laplacian_spectrum(&cycle(n), true).unwrap();372 let mut want: Vec<f64> = (0..n)373 .map(|k| 1.0 - (2.0 * PI * k as f64 / n as f64).cos())374 .collect();375 want.sort_by(|a, b| a.partial_cmp(b).unwrap());376 for (got, expected) in values.iter().zip(&want) {377 assert!((got - expected).abs() < 1e-10, "n={n} {got} {expected}");378 }379 }380 }381382 #[test]383 fn a_random_symmetric_matrix_keeps_its_trace_and_frobenius_norm() {384 let n = 40;385 let mut seed = 0x2545f4914f6cdd1du64;386 let mut next = || {387 seed ^= seed << 13;388 seed ^= seed >> 7;389 seed ^= seed << 17;390 (seed >> 11) as f64 / (1u64 << 53) as f64 - 0.5391 };392 let mut m = vec![vec![0.0; n]; n];393 for i in 0..n {394 for j in i..n {395 let v = next();396 m[i][j] = v;397 m[j][i] = v;398 }399 }400 let trace: f64 = (0..n).map(|i| m[i][i]).sum();401 let frobenius: f64 = m.iter().flatten().map(|v| v * v).sum();402 let values = symmetric_eigenvalues(&m).unwrap();403 assert_eq!(values.len(), n);404 assert!((values.iter().sum::<f64>() - trace).abs() < 1e-9);405 assert!((values.iter().map(|v| v * v).sum::<f64>() - frobenius).abs() < 1e-9);406 assert!(values.windows(2).all(|w| w[0] <= w[1]));407 }408409 #[test]410 fn a_three_by_three_matches_its_known_roots() {411 let m = vec![412 vec![2.0, 1.0, 0.0],413 vec![1.0, 2.0, 1.0],414 vec![0.0, 1.0, 2.0],415 ];416 let values = symmetric_eigenvalues(&m).unwrap();417 let root = 2.0f64.sqrt();418 for (got, want) in values.iter().zip([2.0 - root, 2.0, 2.0 + root]) {419 assert!((got - want).abs() < 1e-12, "{got} {want}");420 }421 }422423 #[test]424 fn a_crooked_matrix_is_refused() {425 assert!(symmetric_eigenvalues(&[vec![1.0, 2.0]]).is_err());426 assert!(symmetric_eigenvalues(&[vec![1.0, 2.0], vec![3.0, 1.0]]).is_err());427 let mut lonely = Network::new(1);428 lonely.add_node(vec![0.0]).unwrap();429 assert!(laplacian(&lonely, true).is_err());430 assert!(laplacian(&lonely, false).is_ok());431 }432433 #[test]434 fn the_clusters_and_multiplicities_split_a_hand_made_list() {435 let values = [0.0, 1e-15, 2e-15, 1.0, 1.0 + 1e-13, 1.0 - 1e-13, 2.0];436 let groups = clusters(&values, 1e-9);437 assert_eq!(438 groups.iter().map(|g| g.1).collect::<Vec<usize>>(),439 [3, 3, 1]440 );441 assert_eq!(groups.len(), 3);442 assert!(groups[0].0.abs() < 1e-14);443 assert!((groups[1].0 - 1.0).abs() < 1e-14);444 assert_eq!(multiplicity(&values, 1.0, 1e-12), 3);445 assert_eq!(multiplicity(&values, 1.0, 1e-14), 1);446 assert_eq!(multiplicity(&values, 0.0, 1e-9), 3);447 assert_eq!(clusters(&[], 1e-9).len(), 0);448 assert_eq!(spectral_exponent(&[0.0, 0.0], 0.1), None);449 assert_eq!(spectral_fit(&[0.0, 0.0], 0.1), None);450 }451452 #[test]453 fn a_power_law_staircase_returns_its_slope_and_intercept() {454 let total = 200;455 let power = 0.75;456 let values: Vec<f64> = (0..total)457 .map(|j| ((j + 1) as f64 / total as f64).powf(1.0 / power))458 .collect();459 let (intercept, slope, fitted) = spectral_fit(&values, 0.1).unwrap();460 assert!((slope - power).abs() < 1e-9);461 assert!(intercept.abs() < 1e-9);462 assert_eq!(fitted, 20);463 assert_eq!(spectral_points(&values).len(), total);464 assert!((spectral_exponent(&values, 0.1).unwrap() - 2.0 * power).abs() < 1e-9);465 }466}