roulette.rs

1use mrlycore::errors::{value_error, Result};2use mrlynum::factor::gcd;3use mrlynum::spirograph::{trace, Pencil, Track};4use std::collections::HashMap;56const GRID: (usize, usize) = (8, 1024);7const JOIN: f64 = 1e-5;89// THE NODES1011/// Every crossing of a traced roulette: the curves against themselves, the curves against one another, and how crowded the worst node is.12#[derive(Clone, Debug, Default, PartialEq, Eq)]13pub struct Nodes {14    /// The curves counted, in the order the pencils came in.15    pub curves: usize,16    /// How often each curve crosses itself, curve by curve.17    pub selves: Vec<usize>,18    /// How often each pair of curves crosses, the lower curve first, in lexicographic order.19    pub pairs: Vec<usize>,20    /// The most crossings one node carries: one at a plain double point, and `n(n - 1)/2` where `n` branches meet.21    pub most: usize,22    /// The nodes more than one crossing clusters at.23    pub crowded: usize,24    /// The distinct points the crossings sit at, one for every cluster.25    pub points: usize,26    /// The branches through every node added up, which is the edge count of the picture as a plane graph, `n` at a node where `n` branches meet and `2` times `points` when no node is crowded.27    pub branches: usize,28    /// The segment pairs that meet without crossing: collinear or end to end.29    pub touches: usize,30}3132impl Nodes {33    /// How often the curves `i` and `j` cross, either order, and zero when they are one curve.34    pub fn pair(&self, i: usize, j: usize) -> usize {35        if i == j || i >= self.curves || j >= self.curves {36            return 0;37        }38        let (i, j) = (i.min(j), i.max(j));39        self.pairs[i * self.curves - i * (i + 1) / 2 + j - i - 1]40    }4142    /// Every self crossing.43    pub fn selved(&self) -> usize {44        self.selves.iter().sum()45    }4647    /// Every crossing of two curves.48    pub fn paired(&self) -> usize {49        self.pairs.iter().sum()50    }5152    /// Every crossing, self and pair together, which counts a node where `n` branches meet `n(n - 1)/2` times; `points` is the count of distinct nodes and the two agree exactly when `crowded` is zero.53    pub fn total(&self) -> usize {54        self.selved() + self.paired()55    }56}5758/// Counts the nodes of the roulette the pencils draw on the track: `mrlynum::spirograph::trace` at `samples` points a pencil, every pair of polyline segments tested for a proper crossing by orientation signs on a grid of buckets, and crossings within `tol` of the picture's longer side read as one node. A pair of segments is counted in one bucket alone, the first they share, so no crossing is counted twice; the sign of an orientation is `side`, exact for any endpoints whose two differences are exact, which two `f32` endpoints are while the picture's coordinates keep their exponents within 29 of one another, as these pictures do. A seat at the wheel's centre draws one circle `b` times over and the count is meaningless there, the passes crossing one another as the sampling wanders.59pub fn nodes(track: &Track, pencils: &[Pencil], samples: usize, tol: f64) -> Result<Nodes> {60    if pencils.is_empty() {61        return value_error("a roulette needs a pencil.");62    }63    if samples < 3 {64        return value_error("a node count needs at least three samples.");65    }66    if !(0.0..1.0).contains(&tol) {67        return value_error("the tolerance is a fraction of the picture's longer side.");68    }69    let trail = trace(track, pencils, samples)?;70    let (curves, steps) = (pencils.len(), samples - 1);71    let at = |k: usize, i: usize| {72        let base = (k * samples + i) * 2;73        [f64::from(trail[base]), f64::from(trail[base + 1])]74    };75    let (low, high) = box_of(&trail);76    let span = (high[0] - low[0]).max(high[1] - low[1]);77    if span <= 0.0 {78        return value_error("the trace has no extent.");79    }80    let size = [81        (high[0] - low[0]).max(span * 1e-9),82        (high[1] - low[1]).max(span * 1e-9),83    ];84    let g = ((curves * steps) as f64).sqrt() as usize;85    let g = g.clamp(GRID.0, GRID.1);86    let cell = |p: [f64; 2]| {87        let place = |v: f64, k: usize| (((v - low[k]) / size[k] * g as f64) as usize).min(g - 1);88        (place(p[0], 0), place(p[1], 1))89    };90    let joined: Vec<bool> = (0..curves)91        .map(|k| {92            let (first, last) = (at(k, 0), at(k, steps));93            (first[0] - last[0]).hypot(first[1] - last[1]) < span * JOIN94        })95        .collect();96    let mut buckets: Vec<Vec<u32>> = vec![Vec::new(); g * g];97    for k in 0..curves {98        for i in 0..steps {99            let (lo, hi) = ends(at(k, i), at(k, i + 1));100            let ((x0, y0), (x1, y1)) = (cell(lo), cell(hi));101            for x in x0..=x1 {102                for y in y0..=y1 {103                    buckets[x * g + y].push((k * steps + i) as u32);104                }105            }106        }107    }108    let mut out = Nodes {109        curves,110        selves: vec![0; curves],111        pairs: vec![0; curves * (curves - 1) / 2],112        ..Nodes::default()113    };114    let mut found: Vec<[f64; 2]> = Vec::new();115    for x in 0..g {116        for y in 0..g {117            let list = &buckets[x * g + y];118            for (u, &left) in list.iter().enumerate() {119                for &right in &list[u + 1..] {120                    let (p, q) = (left as usize, right as usize);121                    let ((kp, ip), (kq, iq)) = ((p / steps, p % steps), (q / steps, q % steps));122                    if kp == kq && neighbours(ip, iq, steps, joined[kp]) {123                        continue;124                    }125                    let (a, b) = (at(kp, ip), at(kp, ip + 1));126                    let (c, d) = (at(kq, iq), at(kq, iq + 1));127                    let (first, second) = (ends(a, b), ends(c, d));128                    let (one, two) = (cell(first.0), cell(second.0));129                    if (one.0.max(two.0), one.1.max(two.1)) != (x, y) {130                        continue;131                    }132                    let (s1, s2) = (side(a, b, c), side(a, b, d));133                    let (s3, s4) = (side(c, d, a), side(c, d, b));134                    if s1 * s2 < 0 && s3 * s4 < 0 {135                        if kp == kq {136                            out.selves[kp] += 1;137                        } else {138                            let (i, j) = (kp.min(kq), kp.max(kq));139                            out.pairs[i * curves - i * (i + 1) / 2 + j - i - 1] += 1;140                        }141                        found.push(meet(a, b, c, d));142                    } else if s1 * s2 <= 0 && s3 * s4 <= 0 {143                        out.touches += 1;144                    }145                }146            }147        }148    }149    let (most, crowded, points, branches) = crowd(&found, tol * span);150    out.most = most;151    out.crowded = crowded;152    out.points = points;153    out.branches = branches;154    Ok(out)155}156157/// One pencil for every distinct curve, the coincidence law read on the exact seats when `exact` says the seats carry no jitter: the first pencil of each family, in the order they came in. On a circle the seats fall into classes under the rotation group of order `gcd(b, 4)`, which is the clause `mrlynum::spirograph::distinct` and `mrlynum::spirograph::representatives` read; on a line and on a polygon every distinct seat draws its own curve, two seats of one radius on a line drawing translates of one shape and never one curve.158pub fn spread(track: &Track, pencils: &[Pencil], exact: bool) -> Vec<Pencil> {159    if !exact {160        return pencils.to_vec();161    }162    let mut seen: HashMap<(i64, i64), usize> = HashMap::new();163    let mut out = Vec::new();164    for pencil in pencils {165        let key = match track.kind.as_str() {166            "in" | "out" => least(pencil.seat, gcd(track.ratio.1, 4)),167            _ => (out.len() as i64, 1),168        };169        if seen.insert(key, out.len()).is_none() {170            out.push(*pencil);171        }172    }173    out174}175176fn least(seat: (i64, i64), order: usize) -> (i64, i64) {177    let quarter = |(u, v): (i64, i64)| (-v, u);178    let (mut best, mut turned) = (seat, seat);179    for _ in 1..order {180        for _ in 0..4 / order {181            turned = quarter(turned);182        }183        best = best.min(turned);184    }185    best186}187188// THE GEOMETRY189190fn box_of(trail: &[f32]) -> ([f64; 2], [f64; 2]) {191    let mut low = [f64::MAX; 2];192    let mut high = [f64::MIN; 2];193    for pair in trail.chunks_exact(2) {194        for k in 0..2 {195            let value = f64::from(pair[k]);196            low[k] = low[k].min(value);197            high[k] = high[k].max(value);198        }199    }200    (low, high)201}202203fn ends(a: [f64; 2], b: [f64; 2]) -> ([f64; 2], [f64; 2]) {204    (205        [a[0].min(b[0]), a[1].min(b[1])],206        [a[0].max(b[0]), a[1].max(b[1])],207    )208}209210fn neighbours(i: usize, j: usize, steps: usize, joined: bool) -> bool {211    let (i, j) = (i.min(j), i.max(j));212    j == i + 1 || (joined && i == 0 && j == steps - 1)213}214215/// Which side of the line from `a` to `b` the point `c` lies: plus one to the left, minus one to the right, zero on it. The sign is exact whenever the two differences `b - a` and `c - a` are exact, whatever the size of the products: the determinant is taken by the fused multiply-add identity of Kahan, whose error is at most twice the rounding unit times the determinant itself, so it can neither flip a sign nor invent one.216pub fn side(a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> i32 {217    let (ux, uy) = (b[0] - a[0], b[1] - a[1]);218    let (vx, vy) = (c[0] - a[0], c[1] - a[1]);219    let cross = uy * vx;220    let slip = f64::mul_add(uy, vx, -cross);221    let turn = f64::mul_add(ux, vy, -cross) - slip;222    match turn.partial_cmp(&0.0) {223        Some(std::cmp::Ordering::Greater) => 1,224        Some(std::cmp::Ordering::Less) => -1,225        _ => 0,226    }227}228229fn meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> [f64; 2] {230    let run = [b[0] - a[0], b[1] - a[1]];231    let other = [d[0] - c[0], d[1] - c[1]];232    let denominator = run[0] * other[1] - run[1] * other[0];233    let step = ((c[0] - a[0]) * other[1] - (c[1] - a[1]) * other[0]) / denominator;234    [a[0] + step * run[0], a[1] + step * run[1]]235}236237// THE CROWD238239fn crowd(found: &[[f64; 2]], tol: f64) -> (usize, usize, usize, usize) {240    if found.is_empty() {241        return (0, 0, 0, 0);242    }243    if tol <= 0.0 {244        return (1, 0, found.len(), 2 * found.len());245    }246    let key = |p: [f64; 2]| ((p[0] / tol).floor() as i64, (p[1] / tol).floor() as i64);247    let mut cells: HashMap<(i64, i64), Vec<usize>> = HashMap::new();248    for (index, &point) in found.iter().enumerate() {249        cells.entry(key(point)).or_default().push(index);250    }251    let mut parent: Vec<usize> = (0..found.len()).collect();252    for (index, &point) in found.iter().enumerate() {253        let (x, y) = key(point);254        for dx in -1..=1 {255            for dy in -1..=1 {256                for &other in cells.get(&(x + dx, y + dy)).into_iter().flatten() {257                    let far = (found[other][0] - point[0]).hypot(found[other][1] - point[1]);258                    if other != index && far <= tol {259                        join(&mut parent, index, other);260                    }261                }262            }263        }264    }265    let mut sizes: HashMap<usize, usize> = HashMap::new();266    for index in 0..found.len() {267        *sizes.entry(root(&mut parent, index)).or_default() += 1;268    }269    let meeting = |count: usize| {270        let mut n = 2;271        while n * (n - 1) / 2 < count {272            n += 1;273        }274        n275    };276    (277        sizes.values().copied().max().unwrap_or(0),278        sizes.values().filter(|&&count| count > 1).count(),279        sizes.len(),280        sizes.values().map(|&count| meeting(count)).sum(),281    )282}283284fn root(parent: &mut [usize], mut index: usize) -> usize {285    while parent[index] != index {286        parent[index] = parent[parent[index]];287        index = parent[index];288    }289    index290}291292fn join(parent: &mut [usize], a: usize, b: usize) {293    let (a, b) = (root(parent, a), root(parent, b));294    if a != b {295        parent[a.max(b)] = a.min(b);296    }297}298299#[cfg(test)]300mod tests {301    use super::*;302    use mrlynum::spirograph::{distinct, pencils, track};303304    fn carpet() -> Vec<u8> {305        vec![1, 1, 1, 1, 0, 1, 1, 1, 1]306    }307308    #[test]309    fn two_curves_on_one_orbit_cross_twice_the_ratio_whatever_their_seats() {310        let seats = pencils(&carpet(), 3, 3, "fill", 0.3, 0.0, 1).unwrap();311        for a in [3, 4, 5, 6, 7] {312            let path = track("in", a, 1, 4, 1).unwrap();313            let count = nodes(&path, &seats, 3000, 4e-4).unwrap();314            assert_eq!(count.paired(), 56 * a);315            assert_eq!(count.selved(), 0);316            assert_eq!((count.most, count.crowded), (1, 0));317        }318    }319320    #[test]321    fn a_curve_below_the_threshold_crosses_itself_a_times_b_less_one() {322        let seats = pencils(&[1, 0, 0, 0, 0, 0, 0, 1, 0], 3, 3, "fill", 0.9, 0.0, 1).unwrap();323        let path = track("in", 7, 4, 4, 1).unwrap();324        let count = nodes(&path, &seats, 6000, 4e-4).unwrap();325        assert_eq!(count.selves, vec![21, 21]);326        assert_eq!(count.pair(0, 1), 56);327        assert_eq!(count.total(), 98);328    }329330    #[test]331    fn the_orientation_sign_holds_where_the_plain_determinant_rounds_to_nothing() {332        let far = 134_217_728.0;333        let (b, c) = ([far + 1.0, far], [far, far - 1.0]);334        assert_eq!(side([0.0, 0.0], b, c), -1);335        assert_eq!(336            ((b[0] * c[1]) - (b[1] * c[0])).partial_cmp(&0.0),337            Some(std::cmp::Ordering::Equal)338        );339        assert_eq!(side([0.0, 0.0], [1.0, 1.0], [3.0, 3.0]), 0);340        assert_eq!(side([0.0, 0.0], c, b), 1);341    }342343    #[test]344    fn an_alignment_reach_crowds_the_nodes_and_the_points_fall_short() {345        let seats = pencils(&carpet(), 3, 3, "fill", 0.79, 0.0, 1).unwrap();346        let path = track("in", 7, 3, 4, 1).unwrap();347        let count = nodes(&path, &seats, 6000, 4e-4).unwrap();348        assert_eq!((count.total(), count.crowded, count.most), (1288, 28, 6));349        assert_eq!((count.points, count.branches), (1148, 2352));350    }351352    #[test]353    fn the_coincident_pencils_collapse_to_one_pencil_a_curve() {354        let seats = pencils(&carpet(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();355        for (kind, ring, wheel) in [("in", 7, 4), ("in", 7, 2), ("in", 7, 3), ("out", 5, 8)] {356            let path = track(kind, ring, wheel, 4, 1).unwrap();357            assert_eq!(358                spread(&path, &seats, true).len(),359                distinct(&path, &seats, true)360            );361        }362    }363}