tourbillon.rs

1use crate::core::error::{value_error, Result};2use crate::core::Rng;3use crate::math::spin::Blend;4use crate::num::factor::{factorize, gcd, mobius, squarefree};5use crate::num::prime::{is_prime, primes, squares};6use serde::{Deserialize, Serialize};78const SCALE_CAP: usize = 199;9const SIZE_FLOOR: usize = 16;10const SIZE_CAP: usize = 1024;11const SHOWN: usize = 8;12const PEAKS: usize = 3;13const LADDER: usize = 1000;14const PLANE_CAP: usize = 64_000_000;15const PERIOD_CAP: usize = 360;16const TIGHT: f64 = 1e-9;1718// THE LAYERS1920/// One carpet of a stack: the odd scale it is drawn at, the weight the linear blends carry it at and the turn it takes about the centre, in degrees.21#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]22pub struct Layer {23    /// The odd scale, the number of cells across the carpet.24    pub scale: usize,25    /// The weight, scaled so the magnitudes average one.26    pub weight: f64,27    /// The turn about the centre, in degrees.28    pub degrees: f64,29}3031/// The readings of a spun stack: its layers, the first eight scales and angles, the mean and RMS contrast over the disc, that contrast times the root of the layer count, the exact centre value, whether the blend carries the weights, the span the raster covers, the sites inside the disc and the brightest three.32#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]33pub struct Stats {34    /// The count of layers in the stack.35    pub layers: usize,36    /// The first eight scales.37    pub scales: Vec<usize>,38    /// The first eight angles, in degrees.39    pub angles: Vec<f64>,40    /// The mean over the disc.41    pub mean: f64,42    /// The RMS contrast over the disc.43    pub rms: f64,44    /// The RMS contrast times the root of the layer count.45    pub faded: f64,46    /// The exact value at the centre, computed and never sampled.47    pub centre: f64,48    /// Whether the blend carries the weights.49    pub weighted: bool,50    /// The smallest value over the disc.51    pub low: f64,52    /// The largest value over the disc.53    pub high: f64,54    /// The count of sites inside the disc.55    pub inside: usize,56    /// The brightest three sites, each as its unit coordinates and its value.57    pub peaks: Vec<[f64; 3]>,58    /// The least whole number of increments that closes a quarter turn, none past the cap.59    pub period: Option<usize>,60    /// The count of distinct angle classes the layers fall in, read a quarter turn apart.61    pub classes: usize,62    /// The count of layer pairs sharing an angle class.63    pub pairs: usize,64}6566/// One angle of the quarter-turn lattice: the turn in degrees and the ninety a over q that names it.67#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]68pub struct Eye {69    /// The turn in degrees.70    pub angle: f64,71    /// The numerator, ninety times a.72    pub numer: usize,73    /// The denominator q.74    pub denom: usize,75}7677fn golden() -> f64 {78    180.0 * (3.0 - 5f64.sqrt())79}8081// THE QUARTER TURN8283/// The angles a quarter turn shares with itself: ninety a over q for every q up to the cap and every a from zero to four q coprime to it, sorted by angle.84///85/// The odd carpet is unchanged by a quarter turn, so two layers whose angles agree modulo ninety stand on one exact lattice of nodes and the stack stops moving there.86///87/// ```88/// let list = mrlyrs::math::tourbillon::eyes(2);89/// assert_eq!(list.len(), 9);90/// assert_eq!(list[1].angle, 45.0);91/// ```92pub fn eyes(qmax: usize) -> Vec<Eye> {93    let mut out = Vec::new();94    for denom in 1..=qmax {95        for step in 0..=4 * denom {96            if gcd(step as u128, denom as u128) != 1 {97                continue;98            }99            out.push(Eye {100                angle: 90.0 * step as f64 / denom as f64,101                numer: 90 * step,102                denom,103            });104        }105    }106    out.sort_by(|a, b| a.angle.total_cmp(&b.angle));107    out108}109110/// The least whole number of increments that closes a quarter turn, none once the count passes the cap.111///112/// ```113/// assert_eq!(mrlyrs::math::tourbillon::period(18.0), Some(5));114/// assert_eq!(mrlyrs::math::tourbillon::period(1.0), Some(90));115/// ```116pub fn period(increment: f64) -> Option<usize> {117    (1..=PERIOD_CAP).find(|count| {118        let turns = *count as f64 * increment / 90.0;119        (turns - turns.round()).abs() <= TIGHT120    })121}122123fn quarter(degrees: f64) -> f64 {124    let angle = degrees.rem_euclid(90.0);125    if 90.0 - angle <= TIGHT {126        0.0127    } else {128        angle129    }130}131132/// The angle classes of a stack read a quarter turn apart: how many the layers fall in, and how many layer pairs share one.133pub fn sharing(list: &[Layer]) -> (usize, usize) {134    let mut angles: Vec<f64> = list.iter().map(|layer| quarter(layer.degrees)).collect();135    angles.sort_by(f64::total_cmp);136    let mut counts: Vec<usize> = Vec::new();137    let mut held = f64::NEG_INFINITY;138    for angle in angles {139        match counts.last_mut() {140            Some(count) if angle - held <= TIGHT => *count += 1,141            _ => {142                counts.push(1);143                held = angle;144            }145        }146    }147    let pairs = counts.iter().map(|count| count * (count - 1) / 2).sum();148    (counts.len(), pairs)149}150151fn odd_top(top: usize) -> Result<usize> {152    if !(3..=SCALE_CAP).contains(&top) || top.is_multiple_of(2) {153        return value_error(format!(154            "the top scale must be odd between 3 and {SCALE_CAP}."155        ));156    }157    Ok(top)158}159160fn side(size: usize) -> Result<usize> {161    if !(SIZE_FLOOR..=SIZE_CAP).contains(&size) {162        return value_error(format!(163            "the raster must be between {SIZE_FLOOR} and {SIZE_CAP} pixels."164        ));165    }166    Ok(size)167}168169fn blend_of(name: &str) -> Result<Blend> {170    match Blend::named(name) {171        Some(blend) => Ok(blend),172        None => value_error(format!(173            "blend {name:?} is not mean, sum, union, meet, parity or difference."174        )),175    }176}177178fn linear(blend: Blend) -> bool {179    matches!(blend, Blend::Mean | Blend::Sum)180}181182fn chosen(scale: usize, set: &str) -> Result<bool> {183    match set {184        "odd" => Ok(true),185        "primes" => Ok(is_prime(scale)),186        "squarefree" => Ok(squarefree(scale)),187        "prime powers" => Ok(factorize(scale).len() == 1),188        _ => value_error(format!(189            "layer set {set:?} is not odd, primes, squarefree or prime powers."190        )),191    }192}193194fn weight_of(scale: usize, weights: &str) -> Result<f64> {195    match weights {196        "plain" => Ok(1.0),197        "mobius" => Ok(f64::from(mobius(scale))),198        "harmonic" => Ok(1.0 / scale as f64),199        _ => value_error(format!(200            "weights {weights:?} are not plain, mobius or harmonic."201        )),202    }203}204205/// The layers of a stack: every scale one, three, five up to the top the set keeps, each with its weight and its angle.206///207/// The schedule is unspun, degrees, golden, primes, random or gaussian, the layer set is odd, primes, squarefree or prime powers and the weights are plain, mobius or harmonic. The weights are rescaled so their magnitudes sum to the layer count, so the mean blend stays paper coverage and the sum blend stays the inked layer count.208///209/// # Errors210///211/// Errors at an even top or one out of range, or for a schedule, set or weights name the stack does not know.212pub fn layers(213    top: usize,214    schedule: &str,215    increment: f64,216    set: &str,217    weights: &str,218    seed: u32,219) -> Result<Vec<Layer>> {220    let top = odd_top(top)?;221    let ladder = primes(LADDER);222    let mut rng = Rng::new(u64::from(seed));223    let mut out: Vec<Layer> = Vec::new();224    let mut index = 0usize;225    for scale in (1..=top).step_by(2) {226        if !chosen(scale, set)? {227            continue;228        }229        index += 1;230        let degrees = match schedule {231            "unspun" => 0.0,232            "degrees" => index as f64 * increment,233            "golden" => index as f64 * golden(),234            "primes" => match index {235                1 => 0.0,236                _ => ladder.get(index - 2).copied().unwrap_or(0) as f64,237            },238            "random" => rng.unit() * 360.0,239            "gaussian" => {240                squares(scale).map_or(0.0, |(a, b)| (b as f64).atan2(a as f64).to_degrees())241            }242            _ => {243                return value_error(format!(244                "schedule {schedule:?} is not unspun, degrees, golden, primes, random or gaussian."245            ))246            }247        };248        out.push(Layer {249            scale,250            weight: weight_of(scale, weights)?,251            degrees,252        });253    }254    let mass: f64 = out.iter().map(|layer| layer.weight.abs()).sum();255    if mass > 0.0 {256        let factor = out.len() as f64 / mass;257        for layer in &mut out {258            layer.weight *= factor;259        }260    }261    Ok(out)262}263264fn lit(scale: usize) -> f32 {265    f32::from(u8::from(((scale as f64 * 0.5).floor() as i64) & 1 == 1))266}267268fn sample(list: &[Layer], blend: Blend) -> Vec<f32> {269    list.iter()270        .map(|layer| {271            let ink = lit(layer.scale);272            if linear(blend) {273                layer.weight as f32 * ink274            } else {275                ink276            }277        })278        .collect()279}280281// THE RASTER282283fn disc(size: usize) -> Vec<bool> {284    let mut out = vec![false; size * size];285    for row in 0..size {286        let v = (row as f64 + 0.5) / size as f64 - 0.5;287        for column in 0..size {288            let u = (column as f64 + 0.5) / size as f64 - 0.5;289            out[row * size + column] = u * u + v * v <= 0.25;290        }291    }292    out293}294295fn cells(layer: &Layer, size: usize, inside: &[bool], ink: &mut [u8]) {296    let (sin, cos) = layer.degrees.to_radians().sin_cos();297    let scale = layer.scale as f64;298    let step = 1.0 / size as f64;299    let (dx, dy) = (cos * step, -sin * step);300    let edge = 0.5 * step - 0.5;301    for row in 0..size {302        let v = (row as f64 + 0.5) * step - 0.5;303        let mut x = cos * edge + sin * v + 0.5;304        let mut y = cos * v - sin * edge + 0.5;305        for column in 0..size {306            let at = row * size + column;307            let held = inside[at] && (0.0..1.0).contains(&x) && (0.0..1.0).contains(&y);308            ink[at] = u8::from(309                held && (scale * x).floor() as i64 & 1 == 1 && (scale * y).floor() as i64 & 1 == 1,310            );311            x += dx;312            y += dy;313        }314    }315}316317fn edges(size: usize, inside: &[bool], ink: &[u8], mark: &mut [u8]) {318    for row in 0..size {319        for column in 0..size {320            let at = row * size + column;321            if !inside[at] {322                mark[at] = 0;323                continue;324            }325            let right = column + 1 < size && inside[at + 1] && ink[at + 1] != ink[at];326            let below = row + 1 < size && inside[at + size] && ink[at + size] != ink[at];327            mark[at] = u8::from(right || below);328        }329    }330}331332fn corners(layer: &Layer, size: usize, inside: &[bool], mark: &mut [u8]) {333    mark.fill(0);334    let (sin, cos) = layer.degrees.to_radians().sin_cos();335    let scale = layer.scale;336    for down in (1..scale).step_by(2) {337        for across in (1..scale).step_by(2) {338            for (vx, vy) in [339                (across, down),340                (across + 1, down),341                (across, down + 1),342                (across + 1, down + 1),343            ] {344                let x = vx as f64 / scale as f64 - 0.5;345                let y = vy as f64 / scale as f64 - 0.5;346                let u = cos * x - sin * y;347                let v = sin * x + cos * y;348                if u * u + v * v > 0.25 {349                    continue;350                }351                let column = ((u + 0.5) * size as f64 - 0.5).round() as i64;352                let row = ((v + 0.5) * size as f64 - 0.5).round() as i64;353                for dy in -1..=1i64 {354                    for dx in -1..=1i64 {355                        let (a, b) = (column + dx, row + dy);356                        if a < 0 || b < 0 || a >= size as i64 || b >= size as i64 {357                            continue;358                        }359                        let at = b as usize * size + a as usize;360                        if inside[at] {361                            mark[at] = 1;362                        }363                    }364                }365            }366        }367    }368}369370/// Rasters the layers onto a square of the size, every one turned about the centre by its own angle and masked to the inscribed disc, then merged site by site.371///372/// The render mode is cells, the lit squares of the carpet, edges, the boundaries between its cells, or corners, the vertices of its lit squares. The weights ride on the layers under the linear blends only. Sites outside the disc are NaN.373///374/// # Errors375///376/// Errors for a mode that is not cells, edges or corners, or a raster past the site budget.377pub fn stack(list: &[Layer], size: usize, mode: &str, blend: Blend) -> Result<Vec<f32>> {378    if !["cells", "edges", "corners"].contains(&mode) {379        return value_error(format!(380            "render mode {mode:?} is not cells, edges or corners."381        ));382    }383    if list.len() * size * size > PLANE_CAP {384        return value_error(format!(385            "{} layers on a {size} raster passes the {PLANE_CAP} site budget.",386            list.len()387        ));388    }389    let inside = disc(size);390    let mut out = vec![f32::NAN; size * size];391    if list.is_empty() {392        return Ok(out);393    }394    let mut ink = vec![0u8; size * size];395    let mut planes: Vec<Vec<u8>> = Vec::with_capacity(list.len());396    for layer in list {397        let mut mark = vec![0u8; size * size];398        match mode {399            "corners" => corners(layer, size, &inside, &mut mark),400            "edges" => {401                cells(layer, size, &inside, &mut ink);402                edges(size, &inside, &ink, &mut mark);403            }404            _ => cells(layer, size, &inside, &mut mark),405        }406        planes.push(mark);407    }408    let scaled: Vec<f32> = list409        .iter()410        .map(|layer| {411            if linear(blend) {412                layer.weight as f32413            } else {414                1.0415            }416        })417        .collect();418    let mut buffer = vec![0f32; list.len()];419    for at in 0..size * size {420        if !inside[at] {421            continue;422        }423        for (slot, plane) in planes.iter().enumerate() {424            buffer[slot] = scaled[slot] * f32::from(plane[at]);425        }426        out[at] = blend.fold(&buffer);427    }428    Ok(out)429}430431fn peaks(field: &[f32], size: usize, top: usize) -> Vec<[f64; 3]> {432    let block = (size / top.max(1)).max(1);433    let wide = size.div_ceil(block);434    let mut best: Vec<(f64, usize, usize)> = Vec::with_capacity(wide * wide);435    for tile in 0..wide * wide {436        let (bx, by) = (tile % wide * block, tile / wide * block);437        let mut hit: Option<(f64, usize, usize)> = None;438        for row in by..(by + block).min(size) {439            for column in bx..(bx + block).min(size) {440                let value = field[row * size + column];441                if value.is_nan() {442                    continue;443                }444                let value = f64::from(value);445                if hit.is_none_or(|(seen, _, _)| value > seen) {446                    hit = Some((value, column, row));447                }448            }449        }450        if let Some(found) = hit {451            best.push(found);452        }453    }454    best.sort_by(|a, b| b.0.total_cmp(&a.0));455    let apart = block as f64 / size as f64;456    let mut out: Vec<[f64; 3]> = Vec::new();457    for (value, column, row) in best {458        let x = (column as f64 + 0.5) / size as f64;459        let y = (row as f64 + 0.5) / size as f64;460        if out461            .iter()462            .all(|site| (site[0] - x).hypot(site[1] - y) >= apart)463        {464            out.push([x, y, value]);465        }466        if out.len() == PEAKS {467            break;468        }469    }470    out471}472473// THE STACK474475/// Spins the odd parity carpets at the scales one, three, five up to the top into one stack on a square of the size, every layer turned about the centre by its own angle and masked to the inscribed disc, so every pixel sees every layer.476///477/// The schedule is unspun, degrees, golden, primes, random or gaussian, the layer set is odd, primes, squarefree or prime powers, the weights are plain, mobius or harmonic, the render mode is cells, edges or corners and the blend is mean, sum, union, meet, parity or difference. The weights are scaled to average magnitude one and apply only under the linear blends, so mean is paper coverage and sum is the inked layer count. Sites outside the disc are NaN.478///479/// # Errors480///481/// Errors at a raster or top out of range, or for a schedule, set, weights or blend name the stack does not know.482#[allow(clippy::too_many_arguments)]483pub fn field(484    top: usize,485    size: usize,486    schedule: &str,487    increment: f64,488    set: &str,489    weights: &str,490    mode: &str,491    blend: &str,492    seed: u32,493) -> Result<Vec<f32>> {494    let size = side(size)?;495    let blend = blend_of(blend)?;496    let list = layers(top, schedule, increment, set, weights, seed)?;497    stack(&list, size, mode, blend)498}499500/// Reads a spun stack against the schedule that made it: the layer count, the first eight scales and angles, the mean and RMS contrast over the disc, that contrast times the root of the layer count, the exact centre value, whether the blend carries the weights, the span the raster covers and the brightest three sites.501///502/// The centre is computed from the layers at the exact half-half point, never sampled off the raster, so no turn of the layers moves it.503///504/// # Errors505///506/// Errors when the field is not size by size, or a name the stack does not know.507#[allow(clippy::too_many_arguments)]508pub fn stats(509    field: &[f32],510    size: usize,511    top: usize,512    schedule: &str,513    increment: f64,514    set: &str,515    weights: &str,516    blend: &str,517    seed: u32,518) -> Result<Stats> {519    let size = side(size)?;520    if field.len() != size * size {521        return value_error("the field must be size by size with size at least 1.");522    }523    let blend = blend_of(blend)?;524    let list = layers(top, schedule, increment, set, weights, seed)?;525    let seen: Vec<f64> = field526        .iter()527        .filter(|value| !value.is_nan())528        .map(|value| f64::from(*value))529        .collect();530    let count = seen.len().max(1) as f64;531    let mean = seen.iter().sum::<f64>() / count;532    let rms = (seen.iter().map(|value| (value - mean).powi(2)).sum::<f64>() / count).sqrt();533    let low = seen.iter().cloned().fold(f64::INFINITY, f64::min);534    let high = seen.iter().cloned().fold(f64::NEG_INFINITY, f64::max);535    let centre = if list.is_empty() {536        0.0537    } else {538        f64::from(blend.fold(&sample(&list, blend)))539    };540    let (classes, pairs) = sharing(&list);541    Ok(Stats {542        layers: list.len(),543        scales: list.iter().take(SHOWN).map(|layer| layer.scale).collect(),544        angles: list.iter().take(SHOWN).map(|layer| layer.degrees).collect(),545        mean,546        rms,547        faded: rms * (list.len() as f64).sqrt(),548        centre,549        weighted: linear(blend),550        low: if low.is_finite() { low } else { 0.0 },551        high: if high.is_finite() { high } else { 0.0 },552        inside: seen.len(),553        peaks: peaks(field, size, top),554        period: period(increment),555        classes,556        pairs,557    })558}559560#[cfg(test)]561mod tests {562    use super::*;563564    #[test]565    fn the_schedule_keeps_the_odd_scales_and_turns_them_by_the_primes() {566        let odd = layers(55, "primes", 0.0, "odd", "plain", 1).unwrap();567        assert_eq!(odd.len(), 28);568        let scales: Vec<usize> = odd.iter().take(SHOWN).map(|layer| layer.scale).collect();569        assert_eq!(scales, vec![1, 3, 5, 7, 9, 11, 13, 15]);570        let angles: Vec<f64> = odd.iter().take(SHOWN).map(|layer| layer.degrees).collect();571        assert_eq!(angles, vec![0.0, 2.0, 3.0, 5.0, 7.0, 11.0, 13.0, 17.0]);572        assert_eq!(573            layers(199, "unspun", 0.0, "primes", "plain", 1)574                .unwrap()575                .len(),576            45577        );578        let counts: Vec<usize> = ["odd", "primes", "squarefree", "prime powers"]579            .iter()580            .map(|set| layers(55, "unspun", 0.0, set, "plain", 1).unwrap().len())581            .collect();582        assert_eq!(counts, vec![28, 15, 23, 19]);583        assert!(layers(54, "unspun", 0.0, "odd", "plain", 1).is_err());584        assert!(layers(201, "unspun", 0.0, "odd", "plain", 1).is_err());585        assert!(layers(55, "spiral", 0.0, "odd", "plain", 1).is_err());586        assert!(layers(55, "unspun", 0.0, "even", "plain", 1).is_err());587        assert!(layers(55, "unspun", 0.0, "odd", "zeta", 1).is_err());588    }589590    #[test]591    fn the_unspun_stack_counts_eighteen_layers_at_its_brightest() {592        let list = layers(55, "unspun", 0.0, "odd", "plain", 1).unwrap();593        let raster = stack(&list, 512, "cells", Blend::Sum).unwrap();594        assert_eq!(raster.len(), 512 * 512);595        let high = raster596            .iter()597            .filter(|value| !value.is_nan())598            .cloned()599            .fold(f32::NEG_INFINITY, f32::max);600        assert_eq!(high, 18.0);601        assert!(raster[0].is_nan());602        let read = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "sum", 1).unwrap();603        assert_eq!((read.high, read.low), (18.0, 0.0));604        assert!(stack(&list, 512, "dots", Blend::Sum).is_err());605    }606607    #[test]608    fn the_field_and_the_stats_hold_the_centre_at_half_the_layers() {609        let raster = field(55, 512, "unspun", 0.0, "odd", "plain", "cells", "mean", 1).unwrap();610        let read = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "mean", 1).unwrap();611        assert_eq!(read.layers, 28);612        assert_eq!(read.centre, 14.0 / 28.0);613        assert_eq!(format!("{:.6}", read.high), format!("{:.6}", 18.0 / 28.0));614        assert!(read.weighted);615        let counted = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "sum", 1).unwrap();616        assert_eq!(counted.centre, 14.0);617        let folded = stats(&raster, 512, 55, "unspun", 0.0, "odd", "plain", "parity", 1).unwrap();618        assert_eq!(folded.centre, 0.0);619        assert!(!folded.weighted);620        assert!(field(55, 4096, "unspun", 0.0, "odd", "plain", "cells", "mean", 1).is_err());621        assert!(field(55, 512, "unspun", 0.0, "odd", "plain", "cells", "blur", 1).is_err());622        assert!(stats(&raster, 256, 55, "unspun", 0.0, "odd", "plain", "mean", 1).is_err());623        let odd = field(51, 64, "unspun", 0.0, "odd", "plain", "cells", "sum", 1).unwrap();624        assert_eq!(625            stats(&odd, 64, 51, "unspun", 0.0, "odd", "plain", "sum", 1)626                .unwrap()627                .centre,628            13.0629        );630        assert_eq!(631            stats(&odd, 64, 51, "unspun", 0.0, "odd", "plain", "parity", 1)632                .unwrap()633                .centre,634            1.0635        );636    }637638    #[test]639    fn no_turn_of_the_layers_moves_the_centre() {640        let centre = |schedule: &str, increment: f64| {641            let raster = field(642                55, 64, schedule, increment, "odd", "plain", "cells", "mean", 1,643            )644            .unwrap();645            stats(646                &raster, 64, 55, schedule, increment, "odd", "plain", "mean", 1,647            )648            .unwrap()649            .centre650        };651        for step in 0..8 {652            assert_eq!(centre("degrees", f64::from(step)), 14.0 / 28.0);653        }654        for schedule in ["golden", "primes", "random", "gaussian"] {655            assert_eq!(centre(schedule, 0.0), 14.0 / 28.0);656        }657    }658659    #[test]660    fn the_three_modes_read_the_same_layers_apart() {661        let shares: Vec<String> = ["cells", "edges", "corners"]662            .iter()663            .map(|mode| {664                let raster =665                    field(55, 256, "golden", 0.0, "odd", "plain", mode, "mean", 1).unwrap();666                let read =667                    stats(&raster, 256, 55, "golden", 0.0, "odd", "plain", "mean", 1).unwrap();668                format!("{:.6}", read.mean)669            })670            .collect();671        assert_eq!(shares.join(" "), "0.231722 0.111924 0.141211");672    }673674    #[test]675    fn the_quarter_turn_lattice_lists_its_angles_smallest_first() {676        let list = eyes(12);677        assert_eq!(list.len(), 185);678        let first: Vec<String> = list679            .iter()680            .take(8)681            .map(|eye| format!("{:.6}", eye.angle))682            .collect();683        assert_eq!(684            first.join(" "),685            "0.000000 7.500000 8.181818 9.000000 10.000000 11.250000 12.857143 15.000000"686        );687        assert_eq!((list[1].numer, list[1].denom), (90, 12));688        assert!(list.windows(2).all(|pair| pair[0].angle <= pair[1].angle));689    }690691    #[test]692    fn a_ninetieth_increment_folds_the_layers_into_angle_classes() {693        let read = |increment: f64, schedule: &str| {694            let raster = field(695                55, 64, schedule, increment, "odd", "plain", "cells", "mean", 1,696            )697            .unwrap();698            stats(699                &raster, 64, 55, schedule, increment, "odd", "plain", "mean", 1,700            )701            .unwrap()702        };703        let coarse = read(18.0, "degrees");704        assert_eq!(coarse.period, Some(5));705        assert_eq!((coarse.classes, coarse.pairs), (5, 65));706        let fine = read(7.5, "degrees");707        assert_eq!(fine.period, Some(12));708        assert_eq!(fine.classes, 12);709        let loose = read(0.0, "golden");710        assert_eq!((loose.classes, loose.pairs), (28, 0));711    }712}