spectrum.rs
16.2 kB · rust · 493 lines
1#![allow(clippy::needless_range_loop)]23use crate::core::error::{value_error, Result};4use crate::math::graph::models::Network;5use std::cmp::Ordering;67const SWEEPS: usize = 60;8const ZERO: f64 = 1e-12;910// SOLVER1112fn tridiagonalise(a: &mut [Vec<f64>]) -> (Vec<f64>, Vec<f64>) {13 let n = a.len();14 let mut d = vec![0.0; n];15 let mut e = vec![0.0; n];16 for i in (1..n).rev() {17 let l = i - 1;18 let mut h = 0.0;19 let mut scale = 0.0;20 if l > 0 {21 for k in 0..=l {22 scale += a[i][k].abs();23 }24 }25 if l == 0 || scale == 0.0 {26 e[i] = a[i][l];27 continue;28 }29 for k in 0..=l {30 a[i][k] /= scale;31 h += a[i][k] * a[i][k];32 }33 let mut f = a[i][l];34 let g = if f >= 0.0 { -h.sqrt() } else { h.sqrt() };35 e[i] = scale * g;36 h -= f * g;37 a[i][l] = f - g;38 f = 0.0;39 for j in 0..=l {40 let mut sum = 0.0;41 for k in 0..=j {42 sum += a[j][k] * a[i][k];43 }44 for k in (j + 1)..=l {45 sum += a[k][j] * a[i][k];46 }47 e[j] = sum / h;48 f += e[j] * a[i][j];49 }50 let hh = f / (h + h);51 for j in 0..=l {52 let f = a[i][j];53 let g = e[j] - hh * f;54 e[j] = g;55 for k in 0..=j {56 a[j][k] -= f * e[k] + g * a[i][k];57 }58 }59 }60 e[0] = 0.0;61 for i in 0..n {62 d[i] = a[i][i];63 }64 (d, e)65}6667fn implicit_ql(d: &mut [f64], e: &mut [f64]) -> Result<()> {68 let n = d.len();69 for i in 1..n {70 e[i - 1] = e[i];71 }72 if n > 0 {73 e[n - 1] = 0.0;74 }75 for l in 0..n {76 let mut sweeps = 0;77 loop {78 let mut m = l;79 while m + 1 < n {80 let dd = d[m].abs() + d[m + 1].abs();81 if e[m].abs() <= f64::EPSILON * dd {82 break;83 }84 m += 1;85 }86 if m == l {87 break;88 }89 sweeps += 1;90 if sweeps > SWEEPS {91 return value_error("The QL iteration did not converge.");92 }93 let mut g = (d[l + 1] - d[l]) / (2.0 * e[l]);94 let mut r = g.hypot(1.0);95 g = d[m] - d[l] + e[l] / (g + if g >= 0.0 { r.abs() } else { -r.abs() });96 let mut s = 1.0;97 let mut c = 1.0;98 let mut p = 0.0;99 let mut underflow = false;100 let mut i = m;101 while i > l {102 i -= 1;103 let f = s * e[i];104 let b = c * e[i];105 r = f.hypot(g);106 e[i + 1] = r;107 if r == 0.0 {108 d[i + 1] -= p;109 e[m] = 0.0;110 underflow = true;111 break;112 }113 s = f / r;114 c = g / r;115 g = d[i + 1] - p;116 r = (d[i] - g) * s + 2.0 * c * b;117 p = s * r;118 d[i + 1] = g + p;119 g = c * r - b;120 }121 if underflow {122 continue;123 }124 d[l] -= p;125 e[l] = g;126 e[m] = 0.0;127 }128 }129 Ok(())130}131132/// Returns the eigenvalues of a dense real symmetric matrix in ascending order.133///134/// The matrix is reduced to tridiagonal form by Householder reflections and then135/// diagonalised by implicit QL with Wilkinson shifts. Errors on a matrix that is136/// not square or whose entries disagree across the diagonal by more than `1e-12`.137///138/// ```139/// let m = vec![vec![2.0, 1.0], vec![1.0, 2.0]];140/// let values = mrlyrs::math::spectrum::symmetric_eigenvalues(&m).unwrap();141/// assert!((values[0] - 1.0).abs() < 1e-12 && (values[1] - 3.0).abs() < 1e-12);142/// ```143///144/// # Errors145///146/// Errors on a matrix that is not square or not symmetric, or a spectrum that is not a number.147pub fn symmetric_eigenvalues(matrix: &[Vec<f64>]) -> Result<Vec<f64>> {148 let n = matrix.len();149 for row in matrix {150 if row.len() != n {151 return value_error("The matrix is not square.");152 }153 }154 for i in 0..n {155 for j in (i + 1)..n {156 if (matrix[i][j] - matrix[j][i]).abs() > ZERO {157 return value_error(format!("The matrix is not symmetric at {i}, {j}."));158 }159 }160 }161 let mut work: Vec<Vec<f64>> = matrix.to_vec();162 let (mut d, mut e) = tridiagonalise(&mut work);163 implicit_ql(&mut d, &mut e)?;164 if d.iter().any(|v| v.is_nan()) {165 return value_error("The spectrum is not a number; the matrix carries a NaN or overflows.");166 }167 d.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));168 Ok(d)169}170171// LAPLACIAN172173/// Builds the Laplacian of a network, the combinatorial `D - A` or the normalised `I - D^-1/2 A D^-1/2`.174///175/// Self-loops are ignored and the degrees are the adjacency row sums. Errors when a176/// node carries no branch and the normalised form is asked for.177///178/// # Errors179///180/// Errors when a node carries no branch under the normalised reading.181pub fn laplacian(network: &Network, normalised: bool) -> Result<Vec<Vec<f64>>> {182 let n = network.nodes.len();183 let mut matrix = vec![vec![0.0; n]; n];184 for branch in &network.branches {185 if branch.parent == branch.child {186 continue;187 }188 matrix[branch.parent][branch.child] += 1.0;189 matrix[branch.child][branch.parent] += 1.0;190 }191 let degree: Vec<f64> = matrix.iter().map(|row| row.iter().sum()).collect();192 for (index, &d) in degree.iter().enumerate() {193 if normalised && d == 0.0 {194 return value_error(format!("Node {index} carries no branch."));195 }196 }197 for i in 0..n {198 for j in 0..n {199 matrix[i][j] = if normalised {200 -matrix[i][j] / (degree[i] * degree[j]).sqrt()201 } else {202 -matrix[i][j]203 };204 }205 matrix[i][i] = if normalised { 1.0 } else { degree[i] };206 }207 Ok(matrix)208}209210/// Returns the ascending Laplacian spectrum of a network, combinatorial or normalised.211///212/// ```213/// let mut net = mrlyrs::math::graph::Network::new(1);214/// net.add_node(vec![0.0]).unwrap();215/// net.add_node(vec![1.0]).unwrap();216/// net.add_branch(0, 1, 1.0).unwrap();217/// let values = mrlyrs::math::spectrum::laplacian_spectrum(&net, true).unwrap();218/// assert!(values[0].abs() < 1e-12 && (values[1] - 2.0).abs() < 1e-12);219/// ```220///221/// # Errors222///223/// Errors when a node carries no branch, or the spectrum is not a number.224pub fn laplacian_spectrum(network: &Network, normalised: bool) -> Result<Vec<f64>> {225 symmetric_eigenvalues(&laplacian(network, normalised)?)226}227228// READINGS229230/// Groups eigenvalues into runs split by consecutive gaps above the tolerance, each run its mean and its size.231///232/// ```233/// let groups = mrlyrs::math::spectrum::clusters(&[0.0, 1e-15, 2.0], 1e-9).unwrap();234/// assert_eq!(groups.len(), 2);235/// assert_eq!(groups[0].1, 2);236/// ```237///238/// # Errors239///240/// Errors on a NaN.241pub fn clusters(eigenvalues: &[f64], tolerance: f64) -> Result<Vec<(f64, usize)>> {242 if eigenvalues.iter().any(|v| v.is_nan()) {243 return value_error("The eigenvalues hold a NaN.");244 }245 let mut values = eigenvalues.to_vec();246 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));247 let mut groups = Vec::new();248 let mut start = 0;249 for index in 1..=values.len() {250 if index == values.len() || values[index] - values[index - 1] > tolerance {251 let run = &values[start..index];252 groups.push((run.iter().sum::<f64>() / run.len() as f64, run.len()));253 start = index;254 }255 }256 Ok(groups)257}258259/// Counts the eigenvalues within the tolerance of a value.260pub fn multiplicity(eigenvalues: &[f64], value: f64, tolerance: f64) -> usize {261 eigenvalues262 .iter()263 .filter(|v| (**v - value).abs() <= tolerance)264 .count()265}266267/// Builds the integrated density of states as points, each an eigenvalue and its rank fraction.268///269/// The eigenvalues are clamped at zero and sorted, exactly equal neighbours collapse to one270/// point at the run's last index over the total count, and points at or below `1e-12` drop.271///272/// ```273/// let points = mrlyrs::math::spectrum::spectral_points(&[0.0, 0.5, 0.5, 2.0]);274/// assert_eq!(points.len(), 2);275/// assert!((points[0].1 - 0.75).abs() < 1e-12);276/// ```277pub fn spectral_points(eigenvalues: &[f64]) -> Vec<(f64, f64)> {278 let mut values: Vec<f64> = eigenvalues.iter().map(|&v| v.max(0.0)).collect();279 values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(Ordering::Equal));280 let total = values.len();281 let mut points = Vec::new();282 let mut index = 0;283 while index < total {284 let mut last = index;285 while last + 1 < total && values[last + 1] == values[index] {286 last += 1;287 }288 if values[last] > ZERO {289 points.push((values[last], (last + 1) as f64 / total as f64));290 }291 index = last + 1;292 }293 points294}295296/// Fits the low window of the integrated density of states in log-log: the intercept, the slope and the fitted count.297///298/// The points come from `spectral_points` and the first `max(floor(window * total), 3)` of299/// them are fitted by least squares. Returns `None` when fewer than two points remain.300pub fn spectral_fit(eigenvalues: &[f64], window: f64) -> Option<(f64, f64, usize)> {301 let points = spectral_points(eigenvalues);302 let xs: Vec<f64> = points.iter().map(|p| p.0.ln()).collect();303 let ys: Vec<f64> = points.iter().map(|p| p.1.ln()).collect();304 let top = ((window * eigenvalues.len() as f64).floor() as usize)305 .max(3)306 .min(xs.len());307 if top < 2 {308 return None;309 }310 let count = top as f64;311 let mean_x: f64 = xs[..top].iter().sum::<f64>() / count;312 let mean_y: f64 = ys[..top].iter().sum::<f64>() / count;313 let covariance: f64 = xs[..top]314 .iter()315 .zip(&ys[..top])316 .map(|(x, y)| (x - mean_x) * (y - mean_y))317 .sum();318 let variance: f64 = xs[..top].iter().map(|x| (x - mean_x) * (x - mean_x)).sum();319 if variance == 0.0 {320 return None;321 }322 let slope = covariance / variance;323 Some((mean_y - slope * mean_x, slope, top))324}325326/// Reads the spectral exponent: twice the log-log slope of the integrated density of states over its low window.327pub fn spectral_exponent(eigenvalues: &[f64], window: f64) -> Option<f64> {328 spectral_fit(eigenvalues, window).map(|(_, slope, _)| 2.0 * slope)329}330331#[cfg(test)]332mod tests {333 use super::*;334 use std::f64::consts::PI;335336 fn path(n: usize) -> Network {337 let mut net = Network::new(1);338 for i in 0..n {339 net.add_node(vec![i as f64]).unwrap();340 }341 for i in 1..n {342 net.add_branch(i - 1, i, 1.0).unwrap();343 }344 net345 }346347 fn complete(n: usize) -> Network {348 let mut net = Network::new(1);349 for i in 0..n {350 net.add_node(vec![i as f64]).unwrap();351 }352 for i in 0..n {353 for j in (i + 1)..n {354 net.add_branch(i, j, 1.0).unwrap();355 }356 }357 net358 }359360 fn cycle(n: usize) -> Network {361 let mut net = path(n);362 net.add_branch(n - 1, 0, 1.0).unwrap();363 net364 }365366 #[test]367 fn the_path_laplacian_reads_its_closed_form() {368 for n in 2..12usize {369 let values = laplacian_spectrum(&path(n), false).unwrap();370 let mut want: Vec<f64> = (0..n)371 .map(|k| 2.0 - 2.0 * (PI * k as f64 / n as f64).cos())372 .collect();373 want.sort_by(|a, b| a.total_cmp(b));374 for (got, expected) in values.iter().zip(&want) {375 assert!((got - expected).abs() < 1e-10, "n={n} {got} {expected}");376 }377 }378 }379380 #[test]381 fn the_complete_laplacian_is_zero_once_and_n_the_rest() {382 for n in 2..10usize {383 let values = laplacian_spectrum(&complete(n), false).unwrap();384 assert!(values[0].abs() < 1e-10, "n={n}");385 for value in &values[1..] {386 assert!((value - n as f64).abs() < 1e-10, "n={n} {value}");387 }388 }389 }390391 #[test]392 fn the_cycle_normalised_laplacian_reads_its_closed_form() {393 for n in 3..12usize {394 let values = laplacian_spectrum(&cycle(n), true).unwrap();395 let mut want: Vec<f64> = (0..n)396 .map(|k| 1.0 - (2.0 * PI * k as f64 / n as f64).cos())397 .collect();398 want.sort_by(|a, b| a.total_cmp(b));399 for (got, expected) in values.iter().zip(&want) {400 assert!((got - expected).abs() < 1e-10, "n={n} {got} {expected}");401 }402 }403 }404405 #[test]406 fn a_random_symmetric_matrix_keeps_its_trace_and_frobenius_norm() {407 let n = 40;408 let mut seed = 0x2545f4914f6cdd1du64;409 let mut next = || {410 seed ^= seed << 13;411 seed ^= seed >> 7;412 seed ^= seed << 17;413 (seed >> 11) as f64 / (1u64 << 53) as f64 - 0.5414 };415 let mut m = vec![vec![0.0; n]; n];416 for i in 0..n {417 for j in i..n {418 let v = next();419 m[i][j] = v;420 m[j][i] = v;421 }422 }423 let trace: f64 = (0..n).map(|i| m[i][i]).sum();424 let frobenius: f64 = m.iter().flatten().map(|v| v * v).sum();425 let values = symmetric_eigenvalues(&m).unwrap();426 assert_eq!(values.len(), n);427 assert!((values.iter().sum::<f64>() - trace).abs() < 1e-9);428 assert!((values.iter().map(|v| v * v).sum::<f64>() - frobenius).abs() < 1e-9);429 assert!(values.windows(2).all(|w| w[0] <= w[1]));430 }431432 #[test]433 fn a_three_by_three_matches_its_known_roots() {434 let m = vec![435 vec![2.0, 1.0, 0.0],436 vec![1.0, 2.0, 1.0],437 vec![0.0, 1.0, 2.0],438 ];439 let values = symmetric_eigenvalues(&m).unwrap();440 let root = 2.0f64.sqrt();441 for (got, want) in values.iter().zip([2.0 - root, 2.0, 2.0 + root]) {442 assert!((got - want).abs() < 1e-12, "{got} {want}");443 }444 }445446 #[test]447 fn refuses_a_crooked_matrix_a_nan_and_a_node_without_a_branch() {448 assert!(symmetric_eigenvalues(&[vec![1.0, 2.0]]).is_err());449 assert!(symmetric_eigenvalues(&[vec![1.0, 2.0], vec![3.0, 1.0]]).is_err());450 assert!(symmetric_eigenvalues(&[vec![f64::NAN]]).is_err());451 assert!(symmetric_eigenvalues(&[vec![f64::NAN, 0.0], vec![0.0, 1.0]]).is_err());452 assert!(clusters(&[0.0, f64::NAN], 1e-9).is_err());453 let mut lonely = Network::new(1);454 lonely.add_node(vec![0.0]).unwrap();455 assert!(laplacian(&lonely, true).is_err());456 assert!(laplacian(&lonely, false).is_ok());457 assert!(laplacian_spectrum(&lonely, true).is_err());458 }459460 #[test]461 fn the_clusters_and_multiplicities_split_a_hand_made_list() {462 let values = [0.0, 1e-15, 2e-15, 1.0, 1.0 + 1e-13, 1.0 - 1e-13, 2.0];463 let groups = clusters(&values, 1e-9).unwrap();464 assert_eq!(465 groups.iter().map(|g| g.1).collect::<Vec<usize>>(),466 [3, 3, 1]467 );468 assert_eq!(groups.len(), 3);469 assert!(groups[0].0.abs() < 1e-14);470 assert!((groups[1].0 - 1.0).abs() < 1e-14);471 assert_eq!(multiplicity(&values, 1.0, 1e-12), 3);472 assert_eq!(multiplicity(&values, 1.0, 1e-14), 1);473 assert_eq!(multiplicity(&values, 0.0, 1e-9), 3);474 assert_eq!(clusters(&[], 1e-9).unwrap().len(), 0);475 assert_eq!(spectral_exponent(&[0.0, 0.0], 0.1), None);476 assert_eq!(spectral_fit(&[0.0, 0.0], 0.1), None);477 }478479 #[test]480 fn a_power_law_staircase_returns_its_slope_and_intercept() {481 let total = 200;482 let power = 0.75;483 let values: Vec<f64> = (0..total)484 .map(|j| ((j + 1) as f64 / total as f64).powf(1.0 / power))485 .collect();486 let (intercept, slope, fitted) = spectral_fit(&values, 0.1).unwrap();487 assert!((slope - power).abs() < 1e-9);488 assert!(intercept.abs() < 1e-9);489 assert_eq!(fitted, 20);490 assert_eq!(spectral_points(&values).len(), total);491 assert!((spectral_exponent(&values, 0.1).unwrap() - 2.0 * power).abs() < 1e-9);492 }493}