spectral.rs
5.8 kB · rust · 200 lines
1use crate::design::Graph;2use faer::{Mat, Side};3use mrlycore::Rng;45pub const MODES: usize = 4;6pub const DENSE_LIMIT: usize = 2000;7const BLOCK: usize = 6;8const ROUNDS: usize = 14;9const SOLVE_TOL: f64 = 1e-9;10const SOLVE_CAP: usize = 200000;11const RITZ_TOL: f64 = 1e-11;1213pub fn apply(graph: &Graph, x: &[f64], out: &mut [f64]) {14 for (node, row) in graph.adjacency.iter().enumerate() {15 let mut total = row.len() as f64 * x[node];16 for other in row {17 total -= x[*other as usize];18 }19 out[node] = total;20 }21}2223pub fn dense(graph: &Graph) -> Mat<f64> {24 let n = graph.nodes();25 let mut out = Mat::<f64>::zeros(n, n);26 for (node, row) in graph.adjacency.iter().enumerate() {27 out[(node, node)] = row.len() as f64;28 for other in row {29 out[(node, *other as usize)] = -1.0;30 }31 }32 out33}3435pub fn low(graph: &Graph) -> Vec<f64> {36 let n = graph.nodes();37 let want = MODES.min(n.saturating_sub(1));38 if want == 0 {39 return vec![0.0];40 }41 if n <= DENSE_LIMIT {42 let mut values = dense(graph)43 .as_ref()44 .self_adjoint_eigenvalues(Side::Lower)45 .expect("the dense eigensolver converges");46 values.truncate(want + 1);47 return values;48 }49 let mut values = vec![0.0];50 values.extend(krylov(graph, want));51 values52}5354pub fn slow_mode(graph: &Graph) -> Vec<f64> {55 let eigen = dense(graph)56 .as_ref()57 .self_adjoint_eigen(Side::Lower)58 .expect("the dense eigensolver converges");59 let u = eigen.U();60 (0..graph.nodes()).map(|node| u[(node, 1)]).collect()61}6263fn krylov(graph: &Graph, want: usize) -> Vec<f64> {64 let n = graph.nodes();65 let mut rng = Rng::new(12345);66 let mut block: Vec<Vec<f64>> = (0..BLOCK)67 .map(|_| (0..n).map(|_| rng.unit() * 2.0 - 1.0).collect())68 .collect();69 orthonormalise(&mut block, &[]);70 let mut basis: Vec<Vec<f64>> = Vec::new();71 let mut gram: Vec<Vec<f64>> = Vec::new();72 let mut work = vec![0.0; n];73 let mut latest: Vec<f64> = Vec::new();74 for _ in 0..ROUNDS {75 let mut grown: Vec<Vec<f64>> = block76 .iter()77 .map(|column| solve(graph, column, &mut work))78 .collect();79 orthonormalise(&mut grown, &basis);80 if grown.is_empty() {81 break;82 }83 for column in &grown {84 apply(graph, column, &mut work);85 let mut row: Vec<f64> = basis.iter().map(|seat| dot(seat, &work)).collect();86 row.push(dot(column, &work));87 for (seat, entry) in row.iter().enumerate().take(gram.len()) {88 gram[seat].push(*entry);89 }90 gram.push(row);91 basis.push(column.clone());92 }93 block = grown;94 if basis.len() < want + 2 {95 continue;96 }97 let candidate = ritz(&gram, want);98 let settled = latest.len() == candidate.len()99 && latest100 .iter()101 .zip(&candidate)102 .all(|(old, new)| (old - new).abs() <= RITZ_TOL * new.abs());103 latest = candidate;104 if settled {105 break;106 }107 }108 latest109}110111fn ritz(gram: &[Vec<f64>], want: usize) -> Vec<f64> {112 let size = gram.len();113 let small = Mat::<f64>::from_fn(size, size, |row, column| {114 0.5 * (gram[row][column] + gram[column][row])115 });116 let mut values = small117 .as_ref()118 .self_adjoint_eigenvalues(Side::Lower)119 .expect("the Ritz eigensolver converges");120 values.truncate(want);121 values122}123124fn solve(graph: &Graph, rhs: &[f64], work: &mut [f64]) -> Vec<f64> {125 let n = graph.nodes();126 let mut residual = rhs.to_vec();127 centre(&mut residual);128 let target = dot(&residual, &residual).sqrt();129 let mut solution = vec![0.0; n];130 if target == 0.0 {131 return solution;132 }133 let mut direction = residual.clone();134 let mut square = dot(&residual, &residual);135 for _ in 0..SOLVE_CAP {136 apply(graph, &direction, work);137 let curvature = dot(&direction, work);138 if curvature <= 0.0 {139 break;140 }141 let stride = square / curvature;142 for seat in 0..n {143 solution[seat] += stride * direction[seat];144 residual[seat] -= stride * work[seat];145 }146 centre(&mut residual);147 let next = dot(&residual, &residual);148 if next.sqrt() <= SOLVE_TOL * target {149 break;150 }151 let blend = next / square;152 for seat in 0..n {153 direction[seat] = residual[seat] + blend * direction[seat];154 }155 square = next;156 }157 centre(&mut solution);158 solution159}160161fn dot(left: &[f64], right: &[f64]) -> f64 {162 left.iter().zip(right).map(|(a, b)| a * b).sum()163}164165fn centre(column: &mut [f64]) {166 let mean = column.iter().sum::<f64>() / column.len() as f64;167 for entry in column.iter_mut() {168 *entry -= mean;169 }170}171172fn orthonormalise(block: &mut Vec<Vec<f64>>, against: &[Vec<f64>]) {173 let mut kept: Vec<Vec<f64>> = Vec::new();174 for column in block.iter_mut() {175 for _ in 0..2 {176 centre(column);177 for seat in against.iter().chain(kept.iter()) {178 let overlap = dot(seat, column);179 for (entry, base) in column.iter_mut().zip(seat) {180 *entry -= overlap * base;181 }182 }183 }184 let norm = dot(column, column).sqrt();185 if norm > 1e-8 {186 for entry in column.iter_mut() {187 *entry /= norm;188 }189 kept.push(column.clone());190 }191 }192 *block = kept;193}194195pub fn exponents(coarse: &[f64], fine: &[f64], scale: f64) -> Vec<f64> {196 let modes = MODES.min(coarse.len() - 1).min(fine.len() - 1);197 (1..=modes)198 .map(|mode| (coarse[mode] / fine[mode]).ln() / scale.ln())199 .collect()200}