roulette.rs

1use crate::core::error::{value_error, Result};2use crate::math::spirograph::{trace, Pencil, Track};3use crate::num::factor::gcd;4use serde::{Deserialize, Serialize};5use std::collections::HashMap;67const GRID: (usize, usize) = (8, 1024);8const JOIN: f64 = 1e-5;910// THE NODES1112/// Every crossing of a traced roulette: the curves against themselves, the curves against one another, and how crowded the worst node is.13#[derive(Clone, Debug, Default, PartialEq, Eq, Serialize, Deserialize)]14pub struct Nodes {15    /// The curves counted, in the order the pencils came in.16    pub curves: usize,17    /// How often each curve crosses itself, curve by curve.18    pub selves: Vec<usize>,19    /// How often each pair of curves crosses, the lower curve first, in lexicographic order.20    pub pairs: Vec<usize>,21    /// The most crossings one node carries: one at a plain double point, and `n(n - 1)/2` where `n` branches meet.22    pub most: usize,23    /// The nodes more than one crossing clusters at.24    pub crowded: usize,25    /// The distinct points the crossings sit at, one for every cluster.26    pub points: usize,27    /// 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.28    pub branches: usize,29    /// The segment pairs that meet without crossing: collinear or end to end.30    pub touches: usize,31}3233impl Nodes {34    /// How often the curves `i` and `j` cross, either order, and zero when they are one curve.35    pub fn pair(&self, i: usize, j: usize) -> usize {36        if i == j || i >= self.curves || j >= self.curves {37            return 0;38        }39        let (i, j) = (i.min(j), i.max(j));40        self.pairs[i * self.curves - i * (i + 1) / 2 + j - i - 1]41    }4243    /// Every self crossing.44    pub fn selved(&self) -> usize {45        self.selves.iter().sum()46    }4748    /// Every crossing of two curves.49    pub fn paired(&self) -> usize {50        self.pairs.iter().sum()51    }5253    /// 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.54    pub fn total(&self) -> usize {55        self.selved() + self.paired()56    }57}5859/// Counts the nodes of the roulette the pencils draw on the track: `mrlyrs::math::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.60///61/// ```62/// use mrlyrs::math::spirograph::{pencils, track};63/// let seats = pencils(&[1, 0, 0, 1], 2, 2, "fill", 0.5, 0.0, 1).unwrap();64/// let path = track("in", 3, 1, 4, 1).unwrap();65/// assert_eq!(mrlyrs::math::roulette::nodes(&path, &seats, 600, 4e-4).unwrap().total(), 6);66/// ```67///68/// # Errors69///70/// Errors on an empty pencil list, under three samples, a tolerance outside zero to one, or a trace with no extent.71pub fn nodes(track: &Track, pencils: &[Pencil], samples: usize, tol: f64) -> Result<Nodes> {72    if pencils.is_empty() {73        return value_error("a roulette needs a pencil.");74    }75    if samples < 3 {76        return value_error("a node count needs at least three samples.");77    }78    if !(0.0..1.0).contains(&tol) {79        return value_error("the tolerance is a fraction of the picture's longer side.");80    }81    let trail = trace(track, pencils, samples)?;82    let (curves, steps) = (pencils.len(), samples - 1);83    let at = |k: usize, i: usize| {84        let base = (k * samples + i) * 2;85        [f64::from(trail[base]), f64::from(trail[base + 1])]86    };87    let (low, high) = box_of(&trail);88    let span = (high[0] - low[0]).max(high[1] - low[1]);89    if span <= 0.0 {90        return value_error("the trace has no extent.");91    }92    let size = [93        (high[0] - low[0]).max(span * 1e-9),94        (high[1] - low[1]).max(span * 1e-9),95    ];96    let g = ((curves * steps) as f64).sqrt() as usize;97    let g = g.clamp(GRID.0, GRID.1);98    let cell = |p: [f64; 2]| {99        let place = |v: f64, k: usize| (((v - low[k]) / size[k] * g as f64) as usize).min(g - 1);100        (place(p[0], 0), place(p[1], 1))101    };102    let joined: Vec<bool> = (0..curves)103        .map(|k| {104            let (first, last) = (at(k, 0), at(k, steps));105            (first[0] - last[0]).hypot(first[1] - last[1]) < span * JOIN106        })107        .collect();108    let mut buckets: Vec<Vec<u32>> = vec![Vec::new(); g * g];109    for k in 0..curves {110        for i in 0..steps {111            let (lo, hi) = ends(at(k, i), at(k, i + 1));112            let ((x0, y0), (x1, y1)) = (cell(lo), cell(hi));113            for x in x0..=x1 {114                for y in y0..=y1 {115                    buckets[x * g + y].push((k * steps + i) as u32);116                }117            }118        }119    }120    let mut out = Nodes {121        curves,122        selves: vec![0; curves],123        pairs: vec![0; curves * (curves - 1) / 2],124        ..Nodes::default()125    };126    let mut found: Vec<[f64; 2]> = Vec::new();127    for x in 0..g {128        for y in 0..g {129            let list = &buckets[x * g + y];130            for (u, &left) in list.iter().enumerate() {131                for &right in &list[u + 1..] {132                    let (p, q) = (left as usize, right as usize);133                    let ((kp, ip), (kq, iq)) = ((p / steps, p % steps), (q / steps, q % steps));134                    if kp == kq && neighbours(ip, iq, steps, joined[kp]) {135                        continue;136                    }137                    let (a, b) = (at(kp, ip), at(kp, ip + 1));138                    let (c, d) = (at(kq, iq), at(kq, iq + 1));139                    let (first, second) = (ends(a, b), ends(c, d));140                    let (one, two) = (cell(first.0), cell(second.0));141                    if (one.0.max(two.0), one.1.max(two.1)) != (x, y) {142                        continue;143                    }144                    let (s1, s2) = (side(a, b, c), side(a, b, d));145                    let (s3, s4) = (side(c, d, a), side(c, d, b));146                    if s1 * s2 < 0 && s3 * s4 < 0 {147                        if kp == kq {148                            out.selves[kp] += 1;149                        } else {150                            let (i, j) = (kp.min(kq), kp.max(kq));151                            out.pairs[i * curves - i * (i + 1) / 2 + j - i - 1] += 1;152                        }153                        found.push(meet(a, b, c, d));154                    } else if s1 * s2 <= 0 && s3 * s4 <= 0 {155                        out.touches += 1;156                    }157                }158            }159        }160    }161    let (most, crowded, points, branches) = crowd(&found, tol * span);162    out.most = most;163    out.crowded = crowded;164    out.points = points;165    out.branches = branches;166    Ok(out)167}168169/// 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 `mrlyrs::math::spirograph::distinct` and `mrlyrs::math::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.170pub fn spread(track: &Track, pencils: &[Pencil], exact: bool) -> Vec<Pencil> {171    if !exact {172        return pencils.to_vec();173    }174    let mut seen: HashMap<(i64, i64), usize> = HashMap::new();175    let mut out = Vec::new();176    for pencil in pencils {177        let key = match track.kind.as_str() {178            "in" | "out" => least(pencil.seat, gcd(track.ratio.1 as u128, 4) as usize),179            _ => (out.len() as i64, 1),180        };181        if seen.insert(key, out.len()).is_none() {182            out.push(*pencil);183        }184    }185    out186}187188fn least(seat: (i64, i64), order: usize) -> (i64, i64) {189    let quarter = |(u, v): (i64, i64)| (-v, u);190    let (mut best, mut turned) = (seat, seat);191    for _ in 1..order {192        for _ in 0..4 / order {193            turned = quarter(turned);194        }195        best = best.min(turned);196    }197    best198}199200// THE GEOMETRY201202fn box_of(trail: &[f32]) -> ([f64; 2], [f64; 2]) {203    let mut low = [f64::MAX; 2];204    let mut high = [f64::MIN; 2];205    for pair in trail.chunks_exact(2) {206        for k in 0..2 {207            let value = f64::from(pair[k]);208            low[k] = low[k].min(value);209            high[k] = high[k].max(value);210        }211    }212    (low, high)213}214215fn ends(a: [f64; 2], b: [f64; 2]) -> ([f64; 2], [f64; 2]) {216    (217        [a[0].min(b[0]), a[1].min(b[1])],218        [a[0].max(b[0]), a[1].max(b[1])],219    )220}221222fn neighbours(i: usize, j: usize, steps: usize, joined: bool) -> bool {223    let (i, j) = (i.min(j), i.max(j));224    j == i + 1 || (joined && i == 0 && j == steps - 1)225}226227/// 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.228pub fn side(a: [f64; 2], b: [f64; 2], c: [f64; 2]) -> i32 {229    let (ux, uy) = (b[0] - a[0], b[1] - a[1]);230    let (vx, vy) = (c[0] - a[0], c[1] - a[1]);231    let cross = uy * vx;232    let slip = f64::mul_add(uy, vx, -cross);233    let turn = f64::mul_add(ux, vy, -cross) - slip;234    match turn.partial_cmp(&0.0) {235        Some(std::cmp::Ordering::Greater) => 1,236        Some(std::cmp::Ordering::Less) => -1,237        _ => 0,238    }239}240241fn meet(a: [f64; 2], b: [f64; 2], c: [f64; 2], d: [f64; 2]) -> [f64; 2] {242    let run = [b[0] - a[0], b[1] - a[1]];243    let other = [d[0] - c[0], d[1] - c[1]];244    let denominator = run[0] * other[1] - run[1] * other[0];245    let step = ((c[0] - a[0]) * other[1] - (c[1] - a[1]) * other[0]) / denominator;246    [a[0] + step * run[0], a[1] + step * run[1]]247}248249// THE CROWD250251fn crowd(found: &[[f64; 2]], tol: f64) -> (usize, usize, usize, usize) {252    if found.is_empty() {253        return (0, 0, 0, 0);254    }255    if tol <= 0.0 {256        return (1, 0, found.len(), 2 * found.len());257    }258    let key = |p: [f64; 2]| ((p[0] / tol).floor() as i64, (p[1] / tol).floor() as i64);259    let mut cells: HashMap<(i64, i64), Vec<usize>> = HashMap::new();260    for (index, &point) in found.iter().enumerate() {261        cells.entry(key(point)).or_default().push(index);262    }263    let mut parent: Vec<usize> = (0..found.len()).collect();264    for (index, &point) in found.iter().enumerate() {265        let (x, y) = key(point);266        for dx in -1..=1 {267            for dy in -1..=1 {268                for &other in cells.get(&(x + dx, y + dy)).into_iter().flatten() {269                    let far = (found[other][0] - point[0]).hypot(found[other][1] - point[1]);270                    if other != index && far <= tol {271                        join(&mut parent, index, other);272                    }273                }274            }275        }276    }277    let mut sizes: HashMap<usize, usize> = HashMap::new();278    for index in 0..found.len() {279        *sizes.entry(root(&mut parent, index)).or_default() += 1;280    }281    let meeting = |count: usize| {282        let mut n = 2;283        while n * (n - 1) / 2 < count {284            n += 1;285        }286        n287    };288    (289        sizes.values().copied().max().unwrap_or(0),290        sizes.values().filter(|&&count| count > 1).count(),291        sizes.len(),292        sizes.values().map(|&count| meeting(count)).sum(),293    )294}295296fn root(parent: &mut [usize], mut index: usize) -> usize {297    while parent[index] != index {298        parent[index] = parent[parent[index]];299        index = parent[index];300    }301    index302}303304fn join(parent: &mut [usize], a: usize, b: usize) {305    let (a, b) = (root(parent, a), root(parent, b));306    if a != b {307        parent[a.max(b)] = a.min(b);308    }309}310311#[cfg(test)]312mod tests {313    use super::*;314    use crate::math::spirograph::{distinct, pencils, track};315316    fn carpet() -> Vec<u8> {317        vec![1, 1, 1, 1, 0, 1, 1, 1, 1]318    }319320    #[test]321    fn two_curves_on_one_orbit_cross_twice_the_ratio_whatever_their_seats() {322        let seats = pencils(&carpet(), 3, 3, "fill", 0.3, 0.0, 1).unwrap();323        for a in [3, 4, 5, 6, 7] {324            let path = track("in", a, 1, 4, 1).unwrap();325            let count = nodes(&path, &seats, 3000, 4e-4).unwrap();326            assert_eq!(count.paired(), 56 * a);327            assert_eq!(count.selved(), 0);328            assert_eq!((count.most, count.crowded), (1, 0));329        }330    }331332    #[test]333    fn a_curve_below_the_threshold_crosses_itself_a_times_b_less_one() {334        let seats = pencils(&[1, 0, 0, 0, 0, 0, 0, 1, 0], 3, 3, "fill", 0.9, 0.0, 1).unwrap();335        let path = track("in", 7, 4, 4, 1).unwrap();336        let count = nodes(&path, &seats, 6000, 4e-4).unwrap();337        assert_eq!(count.selves, vec![21, 21]);338        assert_eq!(count.pair(0, 1), 56);339        assert_eq!(count.total(), 98);340    }341342    #[test]343    fn the_orientation_sign_holds_where_the_plain_determinant_rounds_to_nothing() {344        let far = 134_217_728.0;345        let (b, c) = ([far + 1.0, far], [far, far - 1.0]);346        assert_eq!(side([0.0, 0.0], b, c), -1);347        assert_eq!(348            ((b[0] * c[1]) - (b[1] * c[0])).partial_cmp(&0.0),349            Some(std::cmp::Ordering::Equal)350        );351        assert_eq!(side([0.0, 0.0], [1.0, 1.0], [3.0, 3.0]), 0);352        assert_eq!(side([0.0, 0.0], c, b), 1);353    }354355    #[test]356    fn an_alignment_reach_crowds_the_nodes_and_the_points_fall_short() {357        let seats = pencils(&carpet(), 3, 3, "fill", 0.79, 0.0, 1).unwrap();358        let path = track("in", 7, 3, 4, 1).unwrap();359        let count = nodes(&path, &seats, 6000, 4e-4).unwrap();360        assert_eq!((count.total(), count.crowded, count.most), (1288, 28, 6));361        assert_eq!((count.points, count.branches), (1148, 2352));362    }363364    #[test]365    fn the_coincident_pencils_collapse_to_one_pencil_a_curve() {366        let seats = pencils(&carpet(), 3, 3, "fill", 0.9, 0.0, 1).unwrap();367        for (kind, ring, wheel) in [("in", 7, 4), ("in", 7, 2), ("in", 7, 3), ("out", 5, 8)] {368            let path = track(kind, ring, wheel, 4, 1).unwrap();369            assert_eq!(370                spread(&path, &seats, true).len(),371                distinct(&path, &seats, true)372            );373        }374    }375}