demo-spectra.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board};3use mrlymath::{six, three};4use mrlynum::graph::{census, largest_component};5use mrlynum::spectrum::{laplacian_spectrum, spectral_fit, spectral_points};67const WINDOW: f64 = 0.10;89fn main() -> Result<()> {10    let cell = six::cut(&three::create(23, 3, 2, 2)?)?;11    let whole = six::graph::slice_core_graph(&cell)?;12    let pieces = census(&whole).components;13    let network = largest_component(&whole);14    let values = laplacian_spectrum(&network, true)?;15    let points = spectral_points(&values);16    let (intercept, slope, fitted) = spectral_fit(&values, WINDOW).expect("the low window fits");17    assert_eq!(18        (network.nodes.len(), network.branches.len(), pieces),19        (306, 378, 1)20    );21    assert_eq!((points.len(), fitted), (286, 30));22    assert!((2.0 * slope - 1.253284).abs() < 1e-5);2324    let mut board = Board::square();25    let frame = board.frame(0.08);26    let xs: Vec<f64> = points.iter().map(|p| p.0.ln()).collect();27    let ys: Vec<f64> = points.iter().map(|p| p.1.ln()).collect();28    let (left, right) = (xs[0], *xs.last().unwrap());29    let (foot, roof) = (ys[0], 0.0);30    let at = |x: f64, y: f64| {31        (32            frame.x + frame.w * (x - left) / (right - left),33            frame.y + frame.h * (1.0 - (y - foot) / (roof - foot)),34        )35    };36    let ray = |x: f64| intercept + slope * x;37    let reach = |from: f64, to: f64| {38        let low = ((foot - intercept) / slope).max(from);39        let high = ((roof - intercept) / slope).min(to);40        (at(low, ray(low)), at(high, ray(high)))41    };4243    let (edge, _) = at(xs[fitted - 1], roof);44    board.rect(frame.x, frame.y, edge - frame.x, frame.h, ink::panel());45    plot::axis(&mut board, frame, ink::line());46    let mut stair = vec![at(xs[0], ys[0])];47    for index in 1..points.len() {48        stair.push(at(xs[index], ys[index - 1]));49        stair.push(at(xs[index], ys[index]));50    }51    board.polyline(&stair, 4.0, ink::blue());52    let (a, b) = reach(left, right);53    board.segment(a, b, 2.5, ink::fade(ink::yellow(), 0.6));54    let (a, b) = reach(left, xs[fitted - 1]);55    board.segment(a, b, 6.0, ink::yellow());56    save("demo-spectra", &board)?;57    Ok(())58}