walkers.rs
6.1 kB · rust · 194 lines
1use crate::design::{coords, strides};2use crate::gasket::Gasket;3use mrlycore::{Rng, Tensor};45pub const SEED: u64 = 7;6pub const WALKERS: usize = 20000;7pub const FIRST_FIT: f64 = 32.0;8pub const LAST_STEP: f64 = 1.2e5;9pub const RUNGS: usize = 110;1011pub fn ladder() -> Vec<i64> {12 let step = LAST_STEP.log10() / (RUNGS - 1) as f64;13 let mut out: Vec<i64> = (0..RUNGS)14 .map(|index| 10f64.powf(index as f64 * step) as i64)15 .collect();16 out.sort_unstable();17 out.dedup();18 out19}2021pub struct Trace {22 pub times: Vec<f64>,23 pub readings: Vec<f64>,24}2526pub fn grid_walk(grid: &Tensor, rng: &mut Rng, walkers: usize, side_cap: f64) -> Trace {27 let shape = grid.shape.clone();28 let dims = shape.len();29 let stride = strides(&shape);30 let bytes = grid.bytes();31 let mut cells: Vec<usize> = (0..grid.size()).filter(|flat| bytes[*flat] != 0).collect();32 let mut low = vec![usize::MAX; dims];33 let mut high = vec![0usize; dims];34 for flat in &cells {35 for (axis, at) in coords(*flat, &shape).iter().enumerate() {36 low[axis] = low[axis].min(*at);37 high[axis] = high[axis].max(*at);38 }39 }40 let bulk: Vec<usize> = cells41 .iter()42 .copied()43 .filter(|flat| {44 coords(*flat, &shape).iter().enumerate().all(|(axis, at)| {45 let span = high[axis] - low[axis] + 1;46 span <= 16 || (*at >= low[axis] + span / 4 && *at <= high[axis] - span / 4)47 })48 })49 .collect();50 if bulk.len() >= 200 {51 cells = bulk;52 }53 let mut start: Vec<Vec<i64>> = Vec::with_capacity(walkers);54 let mut place: Vec<Vec<i64>> = Vec::with_capacity(walkers);55 let mut seat: Vec<usize> = Vec::with_capacity(walkers);56 for _ in 0..walkers {57 let flat = cells[rng.below(cells.len())];58 let at: Vec<i64> = coords(flat, &shape).iter().map(|v| *v as i64).collect();59 start.push(at.clone());60 place.push(at);61 seat.push(flat);62 }63 let narrow = shape.iter().min().copied().unwrap_or(0) as f64;64 let cap = (narrow / side_cap).powi(2);65 let ladder = ladder();66 let mut trace = Trace {67 times: Vec::new(),68 readings: Vec::new(),69 };70 let mut rung = 0;71 let mut time = 0i64;72 while rung < ladder.len() {73 time += 1;74 for walker in 0..walkers {75 let drawn = rng.below(2 * dims);76 let axis = drawn / 2;77 let step: i64 = if drawn % 2 == 0 { 1 } else { -1 };78 let moved = place[walker][axis] + step;79 if moved < 0 || moved >= shape[axis] as i64 {80 continue;81 }82 let next = (seat[walker] as i64 + step * stride[axis] as i64) as usize;83 if bytes[next] != 0 {84 place[walker][axis] = moved;85 seat[walker] = next;86 }87 }88 if time == ladder[rung] {89 let total: i64 = (0..walkers)90 .map(|walker| {91 (0..dims)92 .map(|axis| {93 let gap = place[walker][axis] - start[walker][axis];94 gap * gap95 })96 .sum::<i64>()97 })98 .sum();99 let reading = total as f64 / walkers as f64;100 trace.times.push(time as f64);101 trace.readings.push(reading);102 rung += 1;103 if reading > cap {104 break;105 }106 }107 }108 trace109}110111pub fn gasket_walk(gasket: &Gasket, rng: &mut Rng, walkers: usize) -> Trace {112 let nodes = gasket.points.len();113 let table: Vec<[i64; 4]> = gasket114 .graph115 .adjacency116 .iter()117 .map(|row| {118 let mut slots = [-1i64; 4];119 for (slot, other) in row.iter().enumerate().take(4) {120 slots[slot] = *other as i64;121 }122 slots123 })124 .collect();125 let start: Vec<usize> = (0..walkers).map(|_| rng.below(nodes)).collect();126 let mut place = start.clone();127 let reach = gasket128 .points129 .iter()130 .flat_map(|point| [point.0, point.1])131 .max()132 .unwrap_or(0) as f64;133 let cap = (reach / 4.0).powi(2);134 let ladder = ladder();135 let mut trace = Trace {136 times: Vec::new(),137 readings: Vec::new(),138 };139 let mut rung = 0;140 let mut time = 0i64;141 while rung < ladder.len() {142 time += 1;143 for walker in 0..walkers {144 let next = table[place[walker]][rng.below(4)];145 if next >= 0 {146 place[walker] = next as usize;147 }148 }149 if time == ladder[rung] {150 let total: i64 = (0..walkers)151 .map(|walker| {152 let here = gasket.points[place[walker]];153 let there = gasket.points[start[walker]];154 (here.0 - there.0).pow(2) + (here.1 - there.1).pow(2)155 })156 .sum();157 let reading = total as f64 / walkers as f64;158 trace.times.push(time as f64);159 trace.readings.push(reading);160 rung += 1;161 if reading > cap {162 break;163 }164 }165 }166 trace167}168169pub fn slope(x: &[f64], y: &[f64]) -> f64 {170 let n = x.len() as f64;171 let mx = x.iter().sum::<f64>() / n;172 let my = y.iter().sum::<f64>() / n;173 let cov: f64 = x.iter().zip(y).map(|(a, b)| (a - mx) * (b - my)).sum();174 let var: f64 = x.iter().map(|a| (a - mx) * (a - mx)).sum();175 cov / var176}177178pub fn fit(trace: &Trace, ceiling: f64) -> (f64, f64) {179 let mut x: Vec<f64> = Vec::new();180 let mut y: Vec<f64> = Vec::new();181 for (time, reading) in trace.times.iter().zip(&trace.readings) {182 if *time >= FIRST_FIT && *reading < ceiling {183 x.push(time.ln());184 y.push(reading.ln());185 }186 }187 if x.len() < 6 {188 return (f64::NAN, f64::NAN);189 }190 let half = x.len() / 2;191 let first = 2.0 / slope(&x[..half], &y[..half]);192 let second = 2.0 / slope(&x[half..], &y[half..]);193 (2.0 / slope(&x, &y), (first - second).abs())194}