tube.rs

1use crate::{checked, Fault};2use mrlyrs::math::bang::Code;3use mrlyrs::math::two;4use wasm_bindgen::prelude::*;56const REACH: usize = 729;7const FLOOR: f64 = 3.0;89// TRANSFORM1011fn lower(f: &[f64], d: &mut [f64], hull: &mut [usize], edge: &mut [f64]) {12    let n = f.len();13    let mut k = 0;14    hull[0] = 0;15    edge[0] = f64::NEG_INFINITY;16    edge[1] = f64::INFINITY;17    for q in 1..n {18        let mut cut;19        loop {20            let p = hull[k];21            cut = ((f[q] + (q * q) as f64) - (f[p] + (p * p) as f64)) / (2 * (q - p)) as f64;22            if k > 0 && cut <= edge[k] {23                k -= 1;24            } else {25                break;26            }27        }28        k += 1;29        hull[k] = q;30        edge[k] = cut;31        edge[k + 1] = f64::INFINITY;32    }33    k = 0;34    for (q, slot) in d.iter_mut().enumerate() {35        while edge[k + 1] < q as f64 {36            k += 1;37        }38        let gap = q as f64 - hull[k] as f64;39        *slot = gap * gap + f[hull[k]];40    }41}4243fn transform(types: &[u8], side: usize) -> Vec<f64> {44    let far = (4 * side * side) as f64;45    let mut square: Vec<f64> = types46        .iter()47        .map(|&b| if b != 0 { 0.0 } else { far })48        .collect();49    let mut lane = vec![0.0; side];50    let mut out = vec![0.0; side];51    let mut hull = vec![0usize; side];52    let mut edge = vec![0.0; side + 1];53    for col in 0..side {54        for row in 0..side {55            lane[row] = square[row * side + col];56        }57        lower(&lane, &mut out, &mut hull, &mut edge);58        for row in 0..side {59            square[row * side + col] = out[row];60        }61    }62    for row in 0..side {63        lane.copy_from_slice(&square[row * side..(row + 1) * side]);64        lower(&lane, &mut out, &mut hull, &mut edge);65        square[row * side..(row + 1) * side].copy_from_slice(&out);66    }67    square.iter().map(|v| v.sqrt()).collect()68}6970// FIELD7172struct Field {73    side: usize,74    number: usize,75    dimension: f64,76    dist: Vec<f64>,77}7879impl Field {80    fn read(code: &str, number: usize, level: usize, base: usize) -> Result<Field, Fault> {81        if number < 2 {82            return Err(Fault::new("the side must be at least two."));83        }84        if level < 1 {85            return Err(Fault::new("the level must be at least one."));86        }87        let side = number88            .checked_pow(level as u32)89            .filter(|&side| side <= REACH)90            .ok_or_else(|| {91                Fault::new(format!(92                    "side {number} at level {level} passes the {REACH} cells a side the tube allows."93                ))94            })?;95        let code = checked(code, 2, base)?;96        let tile = two::create(Code::from(code), number, 1, 0, base)?;97        let digits = tile.types().bytes()?.iter().filter(|&&b| b != 0).count();98        if digits == 0 {99            return Err(Fault::new("the empty design has no tube."));100        }101        let cell = two::create(Code::from(code), number, level, 0, base)?;102        if cell.width() != side || cell.height() != side {103            return Err(Fault::new("the design is not a square grid."));104        }105        let dist = transform(cell.types().bytes()?, side);106        Ok(Field {107            side,108            number,109            dimension: (digits as f64).ln() / (number as f64).ln(),110            dist,111        })112    }113114    fn volume(&self, eps: f64) -> f64 {115        self.dist116            .iter()117            .map(|&d| (eps - d + 1.0).clamp(0.0, 1.0))118            .sum()119    }120}121122// EXPORTS123124/// Reads the distance field of the design's level-L grid: the Euclidean distance in cell widths from every cell centre to the nearest filled cell centre, row-major, zero on the filled cells.125///126/// The transform is exact, the two-pass lower envelope of parabolas, so a cell inside a square127/// hole of side `s` reads the whole number of cells to the hole's wall. The grid is refused past128/// 729 cells a side.129#[wasm_bindgen]130pub fn tube_distance(131    code: &str,132    number: usize,133    level: usize,134    base: usize,135) -> Result<Vec<f32>, Fault> {136    let field = Field::read(code, number, level, base)?;137    Ok(field.dist.iter().map(|&v| v as f32).collect())138}139140/// Measures the inner tube at radius `eps_cells`: the area within that distance of the design, in cell widths squared, the filled cells counted whole.141///142/// A cell carries the share of itself the radius reaches, `eps - dist + 1` clamped to the unit143/// interval, which is exact at whole radii on a design whose holes are squares: the carpet at144/// level 6 and radius 21 cells reads `5912/6561` of the unit square.145#[wasm_bindgen]146pub fn tube_volume(147    code: &str,148    number: usize,149    level: usize,150    base: usize,151    eps_cells: f64,152) -> Result<f64, Fault> {153    Ok(Field::read(code, number, level, base)?.volume(eps_cells))154}155156/// Walks the Minkowski profile across the radii the grid resolves: the pairs `u = ln(1/eps)` and `M = eps^(d-2) V(eps)` in unit-square units, flat, two numbers a sample.157///158/// The radius runs from `1/q` down to three cells, evenly in `u`, and the pairs come back empty159/// when the grid is too coarse to hold that range. `d` is `log_q k`, the design's fractal160/// dimension, so `M` is the Minkowski content reading whose limit the closed profile is.161#[wasm_bindgen]162pub fn tube_profile(163    code: &str,164    number: usize,165    level: usize,166    base: usize,167    samples: usize,168) -> Result<Vec<f32>, Fault> {169    let field = Field::read(code, number, level, base)?;170    let side = field.side as f64;171    let top = side / field.number as f64;172    if top <= FLOOR || samples < 2 {173        return Ok(Vec::new());174    }175    let (low, high) = ((side / top).ln(), (side / FLOOR).ln());176    let mut pairs = Vec::with_capacity(2 * samples);177    for step in 0..samples {178        let u = low + (high - low) * step as f64 / (samples - 1) as f64;179        let cells = side * (-u).exp();180        let eps = cells / side;181        let volume = field.volume(cells) / (side * side);182        pairs.push(u as f32);183        pairs.push((eps.powf(field.dimension - 2.0) * volume) as f32);184    }185    Ok(pairs)186}187188/// Reads the closed limit profile `G(t)` of the interior-hole class at base `q` with `k` filled cells, the hole sum in closed form.189///190/// Every hole of such a design is an isolated open square whose boundary lies in the design, so191/// the inner tube is the exact sum over the levels, `V(eps) = sum k^(m-1) h(q^-m, eps)`, and192/// `M(eps) = eps^(d-2) V(eps)` runs onto `G`, log-periodic with period `ln q`. The carpet is193/// `q = 3`, `k = 8`, where `G(1/3) = G(1) = 379/280` and the swing above zero is the whole proof194/// that the design is not Minkowski measurable.195#[wasm_bindgen]196pub fn tube_closed(q: usize, k: usize, t: f64) -> Result<f64, Fault> {197    if q < 2 || k <= q || k >= q * q {198        return Err(Fault::new(format!(199            "the class needs {q} < k < {}, not {k}.",200            q * q201        )));202    }203    if !(t > 0.0 && t <= 1.0) {204        return Err(Fault::new("the phase must lie in zero to one."));205    }206    let (base, mass) = (q as f64, k as f64);207    let mut step = 0i32;208    while t * base.powi(step) < 0.5 && step < 512 {209        step += 1;210    }211    while step > -512 && t * base.powi(step - 1) >= 0.5 {212        step -= 1;213    }214    let tail = (mass / (base * base)).powi(step - 1);215    let phase = t * base.powi(step) / base;216    let dimension = mass.ln() / base.ln();217    let sum = 1.0 / (base * base - mass) + 4.0 * phase / (mass - base)218        - 4.0 * phase * phase / (mass - 1.0);219    Ok(t.powf(dimension - 2.0) * tail * sum)220}221222/// Decides whether the design's holes are the isolated interior squares the closed profile needs: every empty cell of the level-one tile off the border, and no two of them touching, even at a corner.223#[wasm_bindgen]224pub fn tube_class(code: &str, number: usize, base: usize) -> Result<bool, Fault> {225    if number < 3 {226        return Ok(false);227    }228    let tile = two::create(Code::from(checked(code, 2, base)?), number, 1, 0, base)?;229    let types = tile.types().bytes()?.to_vec();230    if types.len() != number * number {231        return Ok(false);232    }233    let holes: Vec<(usize, usize)> = (0..types.len())234        .filter(|&at| types[at] == 0)235        .map(|at| (at / number, at % number))236        .collect();237    if holes.is_empty() {238        return Ok(false);239    }240    if holes241        .iter()242        .any(|&(row, col)| row == 0 || col == 0 || row + 1 == number || col + 1 == number)243    {244        return Ok(false);245    }246    for (at, &(row, col)) in holes.iter().enumerate() {247        for &(other, next) in &holes[at + 1..] {248            if row.abs_diff(other) <= 1 && col.abs_diff(next) <= 1 {249                return Ok(false);250            }251        }252    }253    Ok(true)254}