weights.rs

1use crate::{checked, rgba, theme, Fault, Pixels};2use mrlymath::two;3use wasm_bindgen::prelude::*;45const REACH: usize = 1024;6const SPREAD: f64 = 30.0;78// OBJECT910struct Weighted {11    number: usize,12    corners: Vec<(usize, usize)>,13    weights: Vec<f64>,14    logs: Vec<f64>,15    rung: f64,16}1718fn corners_of(code: &str, number: usize, base: usize) -> Result<Vec<(usize, usize)>, Fault> {19    if number < 2 {20        return Err(Fault::new("the side must be at least two."));21    }22    let tile = two::create(checked(code, 2, base)?, number, 1, 0, base)?;23    let types = tile.types().bytes().to_vec();24    let corners: Vec<(usize, usize)> = (0..number * number)25        .filter(|&at| types[at] != 0)26        .map(|at| (at / number, at % number))27        .collect();28    if corners.is_empty() {29        return Err(Fault::new("the empty design carries no mass."));30    }31    Ok(corners)32}3334fn ramp(t: f64) -> [u8; 4] {35    let ink = theme();36    let stops = [rgba(ink.ground), rgba(ink.blue), rgba(ink.yellow)];37    let reach = t.clamp(0.0, 1.0) * (stops.len() - 1) as f64;38    let low = (reach.floor() as usize).min(stops.len() - 2);39    let fade = reach - low as f64;40    let mix = |a: u8, b: u8| (f64::from(a) + (f64::from(b) - f64::from(a)) * fade).round() as u8;41    [42        mix(stops[low][0], stops[low + 1][0]),43        mix(stops[low][1], stops[low + 1][1]),44        mix(stops[low][2], stops[low + 1][2]),45        255,46    ]47}4849impl Weighted {50    fn read(code: &str, number: usize, base: usize, weights: &[f64]) -> Result<Weighted, Fault> {51        let corners = corners_of(code, number, base)?;52        if weights.len() != corners.len() {53            return Err(Fault::new(format!(54                "the design fills {} corners and {} weights arrived.",55                corners.len(),56                weights.len()57            )));58        }59        if weights.iter().any(|w| !w.is_finite() || *w <= 0.0) {60            return Err(Fault::new("every weight must be finite and above zero."));61        }62        let total: f64 = weights.iter().sum();63        let weights: Vec<f64> = weights.iter().map(|w| w / total).collect();64        let logs = weights.iter().map(|w| w.ln()).collect();65        Ok(Weighted {66            number,67            corners,68            weights,69            logs,70            rung: (number as f64).ln(),71        })72    }7374    fn tilt(&self, s: f64) -> (f64, f64) {75        let peak = self76            .logs77            .iter()78            .map(|l| s * l)79            .fold(f64::NEG_INFINITY, f64::max);80        let (mut sum, mut drift) = (0.0, 0.0);81        for l in &self.logs {82            let share = (s * l - peak).exp();83            sum += share;84            drift += share * l;85        }86        ((peak + sum.ln()) / self.rung, -drift / sum / self.rung)87    }8889    fn span(&self, level: usize) -> Result<usize, Fault> {90        if level < 1 {91            return Err(Fault::new("the level must be at least one."));92        }93        self.number94            .checked_pow(level as u32)95            .filter(|span| *span <= REACH)96            .ok_or_else(|| {97                Fault::new(format!(98                    "side {} at level {level} passes the {REACH} cells a side the grid allows.",99                    self.number100                ))101            })102    }103104    fn mass(&self, level: usize) -> Result<(usize, Vec<f64>), Fault> {105        let span = self.span(level)?;106        let mut wide = 1usize;107        let mut field = vec![1.0f64];108        for _ in 0..level {109            let next = wide * self.number;110            let mut grown = vec![0.0f64; next * next];111            for row in 0..wide {112                for col in 0..wide {113                    let held = field[row * wide + col];114                    if held == 0.0 {115                        continue;116                    }117                    for (at, &(a, b)) in self.corners.iter().enumerate() {118                        grown[(row * self.number + a) * next + col * self.number + b] =119                            held * self.weights[at];120                    }121                }122            }123            field = grown;124            wide = next;125        }126        Ok((span, field))127    }128}129130// EXPORTS131132/// Lists the filled corners of the design's level-one tile as row and column pairs, two numbers a corner, the order every weight vector is read in.133///134/// The corner count `k` is half the length: a weighted design puts a probability vector on these135/// `k` cells, and the level-`L` cell of the digit word `f_1 ... f_L` carries the product of their136/// weights.137#[wasm_bindgen]138pub fn weights_corners(code: &str, number: usize, base: usize) -> Result<Vec<usize>, Fault> {139    Ok(corners_of(code, number, base)?140        .into_iter()141        .flat_map(|(row, col)| [row, col])142        .collect())143}144145/// Reads the four local dimensions of a weighted design: `alpha_min`, `alpha_max`, the information exponent `alpha(1)` and `tau(0)`.146///147/// `alpha_min = -log_q max_f w_f` and `alpha_max = -log_q min_f w_f` are the ends of the spectrum's148/// support, `alpha(1) = -sum_f w_f log_q w_f` is the exponent almost every point carries, and149/// `tau(0) = log_q k` is the box dimension of the support, which no weighting moves. The weights150/// are given per filled corner in the corner order and normalised to sum to one inside.151#[wasm_bindgen]152pub fn weights_dims(153    code: &str,154    number: usize,155    base: usize,156    weights: &[f64],157) -> Result<Vec<f64>, Fault> {158    let read = Weighted::read(code, number, base, weights)?;159    let peak = read.logs.iter().copied().fold(f64::MIN, f64::max);160    let least = read.logs.iter().copied().fold(f64::MAX, f64::min);161    let held: f64 = read162        .weights163        .iter()164        .zip(&read.logs)165        .map(|(w, l)| w * l)166        .sum();167    Ok(vec![168        -peak / read.rung,169        -least / read.rung,170        -held / read.rung,171        (read.corners.len() as f64).ln() / read.rung,172    ])173}174175/// Builds the level-`L` mass field of a weighted design row-major on the `q^L` by `q^L` grid, the masses summing to one and the empty cells zero.176///177/// The support is the unweighted design's own: only the contraction ratios enter it, and they stay178/// `1/q` at every weighting. The side is refused past 1024 cells.179#[wasm_bindgen]180pub fn weights_mass(181    code: &str,182    number: usize,183    level: usize,184    base: usize,185    weights: &[f64],186) -> Result<Vec<f32>, Fault> {187    let read = Weighted::read(code, number, base, weights)?;188    Ok(read189        .mass(level)?190        .1191        .into_iter()192        .map(|mass| mass as f32)193        .collect())194}195196/// Paints the level-`L` mass field through the ground-blue-gold ramp, every cell read at `(mass / peak)^gamma`, the empty cells left on the ground.197///198/// The gamma is display alone: it lifts the light cells against the heavy ones without touching a199/// number the page prints.200#[wasm_bindgen]201pub fn weights_pixels(202    code: &str,203    number: usize,204    level: usize,205    base: usize,206    weights: &[f64],207    gamma: f64,208) -> Result<Pixels, Fault> {209    if !(gamma.is_finite() && gamma > 0.0) {210        return Err(Fault::new("the gamma must be finite and above zero."));211    }212    let read = Weighted::read(code, number, base, weights)?;213    let (span, field) = read.mass(level)?;214    let peak = field215        .iter()216        .copied()217        .fold(f64::MIN, f64::max)218        .max(f64::MIN_POSITIVE);219    let ground = rgba(theme().ground);220    let colors = field221        .iter()222        .map(|mass| {223            if *mass <= 0.0 {224                ground225            } else {226                ramp((mass / peak).powf(gamma))227            }228        })229        .collect();230    Ok(Pixels::of(span, span, colors))231}232233/// Walks the pressure `tau(s) = log_q sum_f w_f^s` across the range, the pairs `s` and `tau(s)`, two numbers a sample.234///235/// The closed form holds under the open set condition, which every design satisfies with the unit236/// cell. `tau(0) = log_q k` and `tau(1) = 0` at every probability vector.237#[wasm_bindgen]238pub fn weights_pressure(239    code: &str,240    number: usize,241    base: usize,242    weights: &[f64],243    s_lo: f64,244    s_hi: f64,245    samples: usize,246) -> Result<Vec<f32>, Fault> {247    if samples < 2 {248        return Err(Fault::new("the pressure needs at least two samples."));249    }250    if !(s_lo.is_finite() && s_hi.is_finite() && s_lo < s_hi) {251        return Err(Fault::new(252            "the range must run from a low s to a higher one.",253        ));254    }255    let read = Weighted::read(code, number, base, weights)?;256    let mut out = Vec::with_capacity(2 * samples);257    for step in 0..samples {258        let s = s_lo + (s_hi - s_lo) * step as f64 / (samples - 1) as f64;259        out.push(s as f32);260        out.push(read.tilt(s).0 as f32);261    }262    Ok(out)263}264265/// Reads the pressure at one `s` as the four numbers `s`, `tau(s)`, `alpha(s)` and `f(alpha(s))`.266///267/// `alpha(s) = -tau'(s) = -(sum_f w_f^s log_q w_f)/(sum_f w_f^s)` in closed form, and the Legendre268/// transform is attained there, `f(alpha(s)) = alpha(s) s + tau(s)`. At `s = 1` this is the269/// information exponent and its own spectrum value, since `tau(1) = 0`.270#[wasm_bindgen]271pub fn weights_point(272    code: &str,273    number: usize,274    base: usize,275    weights: &[f64],276    s: f64,277) -> Result<Vec<f64>, Fault> {278    if !s.is_finite() {279        return Err(Fault::new("the moment s must be finite."));280    }281    let (tau, alpha) = Weighted::read(code, number, base, weights)?.tilt(s);282    Ok(vec![s, tau, alpha, alpha * s + tau])283}284285/// Walks the multifractal spectrum, the pairs `alpha(s)` and `f(alpha(s))` over the whole range `[alpha_min, alpha_max]`, two numbers a sample.286///287/// The moment is swept on a cubic warp of `s` in minus thirty to thirty, dense where the spectrum288/// turns and sparse on its two flat tails, so the samples land evenly along the curve. Equal289/// weights collapse the range to the single point `(log_q k, log_q k)`.290#[wasm_bindgen]291pub fn weights_spectrum(292    code: &str,293    number: usize,294    base: usize,295    weights: &[f64],296    samples: usize,297) -> Result<Vec<f32>, Fault> {298    if samples < 2 {299        return Err(Fault::new("the spectrum needs at least two samples."));300    }301    let read = Weighted::read(code, number, base, weights)?;302    let mut out = Vec::with_capacity(2 * samples);303    for step in 0..samples {304        let walk = 2.0 * step as f64 / (samples - 1) as f64 - 1.0;305        let s = SPREAD * walk * walk * walk;306        let (tau, alpha) = read.tilt(s);307        out.push(alpha as f32);308        out.push((alpha * s + tau) as f32);309    }310    Ok(out)311}