weights.rs

1use crate::{checked, rgba, theme, Fault, Pixels};2use mrlyrs::math::bang::Code;3use mrlyrs::math::two;4use wasm_bindgen::prelude::*;56const REACH: usize = 1024;7const SPREAD: f64 = 30.0;89// OBJECT1011struct Weighted {12    number: usize,13    corners: Vec<(usize, usize)>,14    weights: Vec<f64>,15    logs: Vec<f64>,16    rung: f64,17}1819fn corners_of(code: &str, number: usize, base: usize) -> Result<Vec<(usize, usize)>, Fault> {20    if number < 2 {21        return Err(Fault::new("the side must be at least two."));22    }23    let tile = two::create(Code::from(checked(code, 2, base)?), number, 1, 0, base)?;24    let types = tile.types().bytes()?.to_vec();25    let corners: Vec<(usize, usize)> = (0..number * number)26        .filter(|&at| types[at] != 0)27        .map(|at| (at / number, at % number))28        .collect();29    if corners.is_empty() {30        return Err(Fault::new("the empty design carries no mass."));31    }32    Ok(corners)33}3435fn ramp(t: f64) -> [u8; 4] {36    let ink = theme();37    let stops = [rgba(ink.ground), rgba(ink.blue), rgba(ink.yellow)];38    let reach = t.clamp(0.0, 1.0) * (stops.len() - 1) as f64;39    let low = (reach.floor() as usize).min(stops.len() - 2);40    let fade = reach - low as f64;41    let mix = |a: u8, b: u8| (f64::from(a) + (f64::from(b) - f64::from(a)) * fade).round() as u8;42    [43        mix(stops[low][0], stops[low + 1][0]),44        mix(stops[low][1], stops[low + 1][1]),45        mix(stops[low][2], stops[low + 1][2]),46        255,47    ]48}4950impl Weighted {51    fn read(code: &str, number: usize, base: usize, weights: &[f64]) -> Result<Weighted, Fault> {52        let corners = corners_of(code, number, base)?;53        if weights.len() != corners.len() {54            return Err(Fault::new(format!(55                "the design fills {} corners and {} weights arrived.",56                corners.len(),57                weights.len()58            )));59        }60        if weights.iter().any(|w| !w.is_finite() || *w <= 0.0) {61            return Err(Fault::new("every weight must be finite and above zero."));62        }63        let total: f64 = weights.iter().sum();64        let weights: Vec<f64> = weights.iter().map(|w| w / total).collect();65        let logs = weights.iter().map(|w| w.ln()).collect();66        Ok(Weighted {67            number,68            corners,69            weights,70            logs,71            rung: (number as f64).ln(),72        })73    }7475    fn tilt(&self, s: f64) -> (f64, f64) {76        let peak = self77            .logs78            .iter()79            .map(|l| s * l)80            .fold(f64::NEG_INFINITY, f64::max);81        let (mut sum, mut drift) = (0.0, 0.0);82        for l in &self.logs {83            let share = (s * l - peak).exp();84            sum += share;85            drift += share * l;86        }87        ((peak + sum.ln()) / self.rung, -drift / sum / self.rung)88    }8990    fn span(&self, level: usize) -> Result<usize, Fault> {91        if level < 1 {92            return Err(Fault::new("the level must be at least one."));93        }94        self.number95            .checked_pow(level as u32)96            .filter(|span| *span <= REACH)97            .ok_or_else(|| {98                Fault::new(format!(99                    "side {} at level {level} passes the {REACH} cells a side the grid allows.",100                    self.number101                ))102            })103    }104105    fn mass(&self, level: usize) -> Result<(usize, Vec<f64>), Fault> {106        let span = self.span(level)?;107        let mut wide = 1usize;108        let mut field = vec![1.0f64];109        for _ in 0..level {110            let next = wide * self.number;111            let mut grown = vec![0.0f64; next * next];112            for row in 0..wide {113                for col in 0..wide {114                    let held = field[row * wide + col];115                    if held == 0.0 {116                        continue;117                    }118                    for (at, &(a, b)) in self.corners.iter().enumerate() {119                        grown[(row * self.number + a) * next + col * self.number + b] =120                            held * self.weights[at];121                    }122                }123            }124            field = grown;125            wide = next;126        }127        Ok((span, field))128    }129}130131// EXPORTS132133/// 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.134///135/// The corner count `k` is half the length: a weighted design puts a probability vector on these136/// `k` cells, and the level-`L` cell of the digit word `f_1 ... f_L` carries the product of their137/// weights.138#[wasm_bindgen]139pub fn weights_corners(code: &str, number: usize, base: usize) -> Result<Vec<usize>, Fault> {140    Ok(corners_of(code, number, base)?141        .into_iter()142        .flat_map(|(row, col)| [row, col])143        .collect())144}145146/// Reads the four local dimensions of a weighted design: `alpha_min`, `alpha_max`, the information exponent `alpha(1)` and `tau(0)`.147///148/// `alpha_min = -log_q max_f w_f` and `alpha_max = -log_q min_f w_f` are the ends of the spectrum's149/// support, `alpha(1) = -sum_f w_f log_q w_f` is the exponent almost every point carries, and150/// `tau(0) = log_q k` is the box dimension of the support, which no weighting moves. The weights151/// are given per filled corner in the corner order and normalised to sum to one inside.152#[wasm_bindgen]153pub fn weights_dims(154    code: &str,155    number: usize,156    base: usize,157    weights: &[f64],158) -> Result<Vec<f64>, Fault> {159    let read = Weighted::read(code, number, base, weights)?;160    let peak = read.logs.iter().copied().fold(f64::MIN, f64::max);161    let least = read.logs.iter().copied().fold(f64::MAX, f64::min);162    let held: f64 = read163        .weights164        .iter()165        .zip(&read.logs)166        .map(|(w, l)| w * l)167        .sum();168    Ok(vec![169        -peak / read.rung,170        -least / read.rung,171        -held / read.rung,172        (read.corners.len() as f64).ln() / read.rung,173    ])174}175176/// 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.177///178/// The support is the unweighted design's own: only the contraction ratios enter it, and they stay179/// `1/q` at every weighting. The side is refused past 1024 cells.180#[wasm_bindgen]181pub fn weights_mass(182    code: &str,183    number: usize,184    level: usize,185    base: usize,186    weights: &[f64],187) -> Result<Vec<f32>, Fault> {188    let read = Weighted::read(code, number, base, weights)?;189    Ok(read190        .mass(level)?191        .1192        .into_iter()193        .map(|mass| mass as f32)194        .collect())195}196197/// 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.198///199/// The gamma is display alone: it lifts the light cells against the heavy ones without touching a200/// number the page prints.201#[wasm_bindgen]202pub fn weights_pixels(203    code: &str,204    number: usize,205    level: usize,206    base: usize,207    weights: &[f64],208    gamma: f64,209) -> Result<Pixels, Fault> {210    if !(gamma.is_finite() && gamma > 0.0) {211        return Err(Fault::new("the gamma must be finite and above zero."));212    }213    let read = Weighted::read(code, number, base, weights)?;214    let (span, field) = read.mass(level)?;215    let peak = field216        .iter()217        .copied()218        .fold(f64::MIN, f64::max)219        .max(f64::MIN_POSITIVE);220    let ground = rgba(theme().ground);221    let colors = field222        .iter()223        .map(|mass| {224            if *mass <= 0.0 {225                ground226            } else {227                ramp((mass / peak).powf(gamma))228            }229        })230        .collect();231    Ok(Pixels::of(span, span, colors))232}233234/// 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.235///236/// The closed form holds under the open set condition, which every design satisfies with the unit237/// cell. `tau(0) = log_q k` and `tau(1) = 0` at every probability vector.238#[wasm_bindgen]239pub fn weights_pressure(240    code: &str,241    number: usize,242    base: usize,243    weights: &[f64],244    s_lo: f64,245    s_hi: f64,246    samples: usize,247) -> Result<Vec<f32>, Fault> {248    if samples < 2 {249        return Err(Fault::new("the pressure needs at least two samples."));250    }251    if !(s_lo.is_finite() && s_hi.is_finite() && s_lo < s_hi) {252        return Err(Fault::new(253            "the range must run from a low s to a higher one.",254        ));255    }256    let read = Weighted::read(code, number, base, weights)?;257    let mut out = Vec::with_capacity(2 * samples);258    for step in 0..samples {259        let s = s_lo + (s_hi - s_lo) * step as f64 / (samples - 1) as f64;260        out.push(s as f32);261        out.push(read.tilt(s).0 as f32);262    }263    Ok(out)264}265266/// Reads the pressure at one `s` as the four numbers `s`, `tau(s)`, `alpha(s)` and `f(alpha(s))`.267///268/// `alpha(s) = -tau'(s) = -(sum_f w_f^s log_q w_f)/(sum_f w_f^s)` in closed form, and the Legendre269/// transform is attained there, `f(alpha(s)) = alpha(s) s + tau(s)`. At `s = 1` this is the270/// information exponent and its own spectrum value, since `tau(1) = 0`.271#[wasm_bindgen]272pub fn weights_point(273    code: &str,274    number: usize,275    base: usize,276    weights: &[f64],277    s: f64,278) -> Result<Vec<f64>, Fault> {279    if !s.is_finite() {280        return Err(Fault::new("the moment s must be finite."));281    }282    let (tau, alpha) = Weighted::read(code, number, base, weights)?.tilt(s);283    Ok(vec![s, tau, alpha, alpha * s + tau])284}285286/// Walks the multifractal spectrum, the pairs `alpha(s)` and `f(alpha(s))` over the whole range `[alpha_min, alpha_max]`, two numbers a sample.287///288/// The moment is swept on a cubic warp of `s` in minus thirty to thirty, dense where the spectrum289/// turns and sparse on its two flat tails, so the samples land evenly along the curve. Equal290/// weights collapse the range to the single point `(log_q k, log_q k)`.291#[wasm_bindgen]292pub fn weights_spectrum(293    code: &str,294    number: usize,295    base: usize,296    weights: &[f64],297    samples: usize,298) -> Result<Vec<f32>, Fault> {299    if samples < 2 {300        return Err(Fault::new("the spectrum needs at least two samples."));301    }302    let read = Weighted::read(code, number, base, weights)?;303    let mut out = Vec::with_capacity(2 * samples);304    for step in 0..samples {305        let walk = 2.0 * step as f64 / (samples - 1) as f64 - 1.0;306        let s = SPREAD * walk * walk * walk;307        let (tau, alpha) = read.tilt(s);308        out.push(alpha as f32);309        out.push((alpha * s + tau) as f32);310    }311    Ok(out)312}