research-weights.rs

1use mrlycore::errors::Result;2use mrlycore::Color;3use mrlyfig::{ink, plot, save, Board, Frame};45const LEVEL: u32 = 8;6const SAMPLES: usize = 720;7const CELLS: [(i64, i64); 3] = [(0, 0), (2, 0), (0, 2)];89fn root(costs: [f64; 3]) -> f64 {10    let (mut lo, mut hi) = (0.0f64, 8.0f64);11    for _ in 0..200 {12        let mid = 0.5 * (lo + hi);13        let mass: f64 = costs.iter().map(|c| (-mid * c).exp()).sum();14        if mass > 1.0 {15            lo = mid;16        } else {17            hi = mid;18        }19    }20    0.5 * (lo + hi)21}2223fn pascal(n: usize) -> Vec<Vec<u128>> {24    let mut rows: Vec<Vec<u128>> = vec![vec![1]];25    for i in 1..=n {26        let mut row = vec![1u128; i + 1];27        for k in 1..i {28            row[k] = rows[i - 1][k - 1] + rows[i - 1][k];29        }30        rows.push(row);31    }32    rows33}3435fn stopping(budget: f64, cheap: f64, dear: f64, binom: &[Vec<u128>]) -> u128 {36    let mut total = 0u128;37    let mut p = 0usize;38    while p as f64 * cheap <= budget {39        let mut q = 0usize;40        while p as f64 * cheap + q as f64 * dear <= budget {41            total += binom[p + q][q] << p;42            q += 1;43        }44        p += 1;45    }46    total47}4849fn mass_side(dear: f64, span: (f64, f64)) -> Vec<f64> {50    let cheap = 2.0f64.ln();51    let delta = root([cheap, cheap, dear]);52    let binom = pascal(2 + (span.1 / cheap) as usize);53    (0..SAMPLES)54        .map(|i| {55            let a = span.0 + (span.1 - span.0) * i as f64 / (SAMPLES - 1) as f64;56            (stopping(a, cheap, dear, &binom) as f64).ln() - delta * a57        })58        .collect()59}6061fn gasket(nums: [u64; 3]) -> Vec<(i64, i64, u128)> {62    let mut pts = vec![(0i64, 0i64, 1u128)];63    for _ in 0..LEVEL {64        let mut next = Vec::with_capacity(pts.len() * 3);65        for &(x, y, m) in &pts {66            for (k, &(dx, dy)) in CELLS.iter().enumerate() {67                next.push((x * 3 + dx, y * 3 + dy, m * nums[k] as u128));68            }69        }70        pts = next;71    }72    pts73}7475fn length_side(nums: [u64; 3], den: u64) -> (f64, Vec<f64>) {76    let pts = gasket(nums);77    assert_eq!(pts.len(), 3usize.pow(LEVEL));78    let whole: u128 = pts.iter().map(|p| p.2).sum();79    assert_eq!(whole, (den as u128).pow(LEVEL));80    let mut rows: Vec<(u128, u128)> = pts81        .iter()82        .map(|&(x, y, m)| ((x * x + y * y) as u128, m))83        .collect();84    rows.sort_by_key(|row| row.0);85    let mut running = 0u128;86    let shells: Vec<(f64, f64)> = rows87        .iter()88        .map(|&(square, m)| {89            running += m;90            (square as f64, running as f64)91        })92        .collect();93    let base = 3.0f64.ln();94    let alpha = -(nums[0] as f64 / den as f64).ln() / base;95    let ys = (0..SAMPLES)96        .map(|i| {97            let u = 4.0 + 4.0 * i as f64 / (SAMPLES - 1) as f64;98            let radius = (u * base).exp();99            let reach = radius * radius;100            let cut = shells.partition_point(|shell| shell.0 <= reach);101            let held = if cut == 0 { 0.0 } else { shells[cut - 1].1 };102            held.ln() - LEVEL as f64 * (den as f64).ln() - alpha * u * base103        })104        .collect();105    (alpha, ys)106}107108fn center(ys: &[f64]) -> Vec<f64> {109    let mean = ys.iter().sum::<f64>() / ys.len() as f64;110    ys.iter().map(|y| y - mean).collect()111}112113fn reach(first: &[f64], second: &[f64]) -> f64 {114    first115        .iter()116        .chain(second.iter())117        .fold(0.0f64, |a, y| a.max(y.abs()))118        * 1.14119}120121fn trace(board: &mut Board, area: Frame, ys: &[f64], span: f64, thick: f64, color: Color) {122    let last = (ys.len() - 1) as f64;123    let pts: Vec<(f64, f64)> = ys124        .iter()125        .enumerate()126        .map(|(i, y)| {127            (128                area.x + area.w * i as f64 / last,129                area.y + area.h * (0.5 - 0.5 * y / span),130            )131        })132        .collect();133    board.polyline(&pts, thick, color);134}135136fn panel(board: &mut Board, area: Frame, first: &[f64], second: &[f64]) {137    plot::axis(board, area, ink::line());138    let inner = area.inset(18.0);139    let mid = inner.y + inner.h / 2.0;140    board.segment((inner.x, mid), (inner.x + inner.w, mid), 1.0, ink::line());141    let span = reach(first, second);142    trace(board, inner, second, span, 3.4, ink::orange());143    trace(board, inner, first, span, 2.1, ink::blue());144}145146fn main() -> Result<()> {147    let window = (48.0, 48.0 + 6.0 * 2.0f64.ln());148    let lattice = mass_side(4.0f64.ln(), window);149    let smooth = mass_side(3.0f64.ln(), window);150    assert!((root([2.0f64.ln(), 2.0f64.ln(), 4.0f64.ln()]) - 1.2715533032).abs() < 1e-9);151    assert!((root([2.0f64.ln(), 2.0f64.ln(), 3.0f64.ln()]) - 1.3646005647).abs() < 1e-9);152153    let (thin, ring_thin) = length_side([2, 2, 1], 5);154    let (fat, ring_fat) = length_side([3, 3, 2], 8);155    assert!((thin - 0.834043767).abs() < 1e-9);156    assert!((fat - 0.892789261).abs() < 1e-9);157158    let mut board = Board::square();159    let frame = board.frame(0.08);160    let half = frame.h / 2.0;161    let top = Frame::new(frame.x, frame.y, frame.w, half).inset(14.0);162    let foot = Frame::new(frame.x, frame.y + half, frame.w, half).inset(14.0);163    panel(&mut board, top, &center(&lattice), &center(&smooth));164    panel(&mut board, foot, &center(&ring_thin), &center(&ring_fat));165    save("research-weights", &board)?;166    Ok(())167}