demo-novelty.rs

1use mrlycore::errors::Result;2use mrlyfig::{ink, plot, save, Board, Frame};3use mrlynum::lattice::totients;4use mrlynum::zeta::{novelty_main, novelty_wave, smoothed_novelty, Line};56const LOW: f64 = 8.0;7const HIGH: f64 = 14.0;8const PER_OCTAVE: usize = 16;9const ONE: usize = 1;10const MANY: usize = 30;11const CURVE: usize = 2400;1213fn trace(area: Frame, span: f64, values: &[f64]) -> Vec<(f64, f64)> {14    let last = (values.len() - 1) as f64;15    values16        .iter()17        .enumerate()18        .map(|(i, v)| {19            (20                area.x + area.w * i as f64 / last,21                area.y + area.h * (0.5 - 0.5 * v / span),22            )23        })24        .collect()25}2627fn main() -> Result<()> {28    let reach = 2.0f64.powf(HIGH + 1.0) as usize;29    let phi = totients(reach);30    assert_eq!(reach, 32_768);31    assert_eq!(&phi[1..=10], &[1, 1, 2, 2, 4, 2, 6, 4, 6, 4]);32    let main = novelty_main();33    let samples = ((HIGH - LOW) * PER_OCTAVE as f64).round() as usize + 1;34    assert_eq!(samples, 97);35    let js: Vec<f64> = (0..samples)36        .map(|k| LOW + k as f64 / PER_OCTAVE as f64)37        .collect();38    let dots: Vec<f64> = js39        .iter()40        .map(|&j| {41            let y = 2.0f64.powf(-j);42            smoothed_novelty(&phi, y, main) / y.powf(1.5)43        })44        .collect();45    let line = Line::new();46    let gammas = line.zeros(MANY);47    assert!((gammas[0] - 14.134_725).abs() < 1e-6);48    let coef = line.novelty_coefficients(&gammas);49    assert!((coef[0].abs() - 0.1879).abs() < 5e-4);50    let wave =51        |count: usize, j: f64| novelty_wave(&gammas[..count], &coef[..count], -j * 2.0f64.ln());52    let peak = dots.iter().fold(0.0f64, |a, v| a.max(v.abs()));53    let miss = |count: usize| {54        js.iter()55            .zip(&dots)56            .map(|(&j, &v)| (v - wave(count, j)).abs())57            .fold(0.0f64, f64::max)58            / peak59    };60    assert!(miss(MANY) < 1e-2);61    assert!(miss(ONE) > 0.3 && miss(ONE) < 0.7);62    let sweep = |count: usize| -> Vec<f64> {63        (0..=CURVE)64            .map(|k| wave(count, LOW + (HIGH - LOW) * k as f64 / CURVE as f64))65            .collect()66    };67    let (one, many) = (sweep(ONE), sweep(MANY));6869    let mut board = Board::square();70    let frame = board.frame(0.08);71    plot::axis(&mut board, frame, ink::line());72    let mid = frame.y + frame.h / 2.0;73    board.segment((frame.x, mid), (frame.x + frame.w, mid), 1.5, ink::line());74    let inner = frame.inset(22.0);75    let span = peak * 1.12;76    board.polyline(&trace(inner, span, &one), 4.0, ink::dim());77    board.polyline(&trace(inner, span, &many), 2.0, ink::orange());78    plot::dots(&mut board, &trace(inner, span, &dots), 5.0, ink::blue());79    save("demo-novelty", &board)?;80    Ok(())81}