Clustering

use clustering.qweruse cluster_kmeans, cluster_choose_k from clustering      // Quillo 0.6: only these

// A shop's customers grouped by what they spend and how often they come (k-means, with k chosen// by the elbow and silhouette), delivery drop-offs grouped into neighbourhoods with DBSCAN in// metres, and a few products grouped by hierarchical clustering.////   quillo run clustering.qwe      then open http://localhost:8080/customers and http://localhost:8080/dropsuse "clustering.qwer"let CUSTOMERS = [    { spend: 12.0, visits: 1 }, { spend: 15.0, visits: 2 }, { spend: 9.5, visits: 1 }, { spend: 14.0, visits: 1 },    { spend: 220.0, visits: 3 }, { spend: 260.0, visits: 2 }, { spend: 240.0, visits: 4 },    { spend: 55.0, visits: 18 }, { spend: 61.0, visits: 22 }, { spend: 48.0, visits: 20 }, { spend: 58.0, visits: 25 },]route GET /customers = {    groups: cluster_kmeans(CUSTOMERS, 3, { fields: ["spend", "visits"], seed: 7 })?,    choose: cluster_choose_k(CUSTOMERS, 5, { fields: ["spend", "visits"] })? | map((r) => `k=${r.k} silhouette=${r.silhouette}`),}// Drop-offs as [lat, lon]: two neighbourhoods of Lyon and one stray address.let DROPS = [[45.7640, 4.8357], [45.7642, 4.8360], [45.7638, 4.8352], [45.7645, 4.8363], [45.7484, 4.8467], [45.7487, 4.8470], [45.7481, 4.8462], [45.7790, 4.8020]]route GET /drops = cluster_dbscan(DROPS, 150.0, 3, { metric: "haversine" })?route GET /products = cluster_hierarchical([[1.0, 1.0], [1.2, 0.9], [5.0, 5.2], [5.1, 4.9], [9.0, 1.0]], 3, 0.0, "average", {})?

//! Data put into groups: k-means (k-means++ starts, seeded so runs repeat) with the elbow and silhouette to choose k, DBSCAN for groups of any shape and the outliers between them (in plain units or as latitude/longitude in metres), and hierarchical clustering cut at a distance or a number of groups; on rows of numbers or maps, scaled or not.//! version: 1.0.0//!//! [dependencies]use serde_json::{json, Value};/// Rows as vectors: lists of numbers, or maps whose `fields` are taken.fn rows(points: &Value, fields: &Value) -> Result<(Vec<Vec<f64>>, Vec<String>), String> {    let a = points.as_array().ok_or("points are a list")?;    if a.is_empty() || a.len() > 200_000 {        return Err("1 to 200,000 points".into());    }    let names: Vec<String> = match fields.as_array() {        Some(f) => f.iter().filter_map(|x| x.as_str().map(str::to_string)).collect(),        None => match &a[0] {            Value::Object(m) => m.iter().filter(|(_, v)| v.is_number()).map(|(k, _)| k.clone()).collect(),            _ => vec![],        },    };    let out: Result<Vec<Vec<f64>>, String> = a        .iter()        .enumerate()        .map(|(i, p)| match p {            Value::Array(xs) => xs.iter().map(|x| x.as_f64().ok_or_else(|| format!("point {}: {} is not a number", i, x))).collect(),            Value::Object(m) => names.iter().map(|k| m.get(k).and_then(Value::as_f64).ok_or_else(|| format!("point {}: no number {:?}", i, k))).collect(),            other => Err(format!("point {}: {} is not a list or a map", i, other)),        })        .collect();    let out = out?;    let d = out[0].len();    if d == 0 || out.iter().any(|r| r.len() != d) {        return Err("every point needs the same number of values".into());    }    Ok((out, names))}/// Each column scaled to mean 0 and standard deviation 1 (so metres and kilograms weigh alike).fn standardize(x: &mut [Vec<f64>]) {    let d = x[0].len();    let n = x.len() as f64;    for j in 0..d {        let m = x.iter().map(|r| r[j]).sum::<f64>() / n;        let s = (x.iter().map(|r| (r[j] - m).powi(2)).sum::<f64>() / n).sqrt();        for r in x.iter_mut() {            r[j] = if s > 0.0 { (r[j] - m) / s } else { 0.0 };        }    }}fn d2(a: &[f64], b: &[f64]) -> f64 {    a.iter().zip(b).map(|(x, y)| (x - y) * (x - y)).sum()}/// Metres between two [lat, lon] points.fn haversine(a: &[f64], b: &[f64]) -> f64 {    let (la1, la2) = (a[0].to_radians(), b[0].to_radians());    let dla = la2 - la1;    let dlo = (b[1] - a[1]).to_radians();    let h = (dla / 2.0).sin().powi(2) + la1.cos() * la2.cos() * (dlo / 2.0).sin().powi(2);    2.0 * 6_371_008.8 * h.sqrt().asin()}struct Rng(u64);impl Rng {    fn next(&mut self) -> f64 {        self.0 ^= self.0 << 13;        self.0 ^= self.0 >> 7;        self.0 ^= self.0 << 17;        (self.0 >> 11) as f64 / (1u64 << 53) as f64    }}/// One k-means run: centres, labels, inertia (the sum of squared distances to the centres).fn kmeans_run(x: &[Vec<f64>], k: usize, seed: u64, iters: usize) -> (Vec<Vec<f64>>, Vec<usize>, f64) {    let mut rng = Rng(seed.wrapping_mul(0x9E3779B97F4A7C15) | 1);    // k-means++: each next centre far from those chosen.    let mut c: Vec<Vec<f64>> = vec![x[(rng.next() * x.len() as f64) as usize % x.len()].clone()];    let mut dist: Vec<f64> = x.iter().map(|p| d2(p, &c[0])).collect();    while c.len() < k {        let total: f64 = dist.iter().sum();        let mut r = rng.next() * total;        let mut pick = x.len() - 1;        for (i, w) in dist.iter().enumerate() {            r -= w;            if r <= 0.0 {                pick = i;                break;            }        }        c.push(x[pick].clone());        for (i, p) in x.iter().enumerate() {            dist[i] = dist[i].min(d2(p, c.last().unwrap()));        }    }    let mut labels = vec![0usize; x.len()];    for _ in 0..iters {        let mut changed = false;        for (i, p) in x.iter().enumerate() {            let best = (0..k).min_by(|a, b| d2(p, &c[*a]).partial_cmp(&d2(p, &c[*b])).unwrap()).unwrap();            if labels[i] != best {                labels[i] = best;                changed = true;            }        }        let d = x[0].len();        let mut sums = vec![vec![0.0; d]; k];        let mut counts = vec![0usize; k];        for (i, p) in x.iter().enumerate() {            counts[labels[i]] += 1;            for j in 0..d {                sums[labels[i]][j] += p[j];            }        }        for g in 0..k {            if counts[g] > 0 {                c[g] = sums[g].iter().map(|s| s / counts[g] as f64).collect();            }        }        if !changed {            break;        }    }    let inertia = x.iter().enumerate().map(|(i, p)| d2(p, &c[labels[i]])).sum();    (c, labels, inertia)}fn round(v: f64) -> f64 {    (v * 1e6).round() / 1e6}/// The mean silhouette (from -1 to 1; above 0.5 the groups are well apart), on at most 2,000 points.fn silhouette(x: &[Vec<f64>], labels: &[i64]) -> Option<f64> {    let n = x.len().min(2000);    let groups: Vec<i64> = {        let mut g: Vec<i64> = labels.iter().copied().filter(|l| *l >= 0).collect();        g.sort();        g.dedup();        g    };    if groups.len() < 2 {        return None;    }    let mut total = 0.0;    let mut m = 0;    for i in 0..n {        if labels[i] < 0 {            continue;        }        let mut mean = |g: i64| -> Option<f64> {            let ds: Vec<f64> = (0..n).filter(|j| *j != i && labels[*j] == g).map(|j| d2(&x[i], &x[j]).sqrt()).collect();            if ds.is_empty() { None } else { Some(ds.iter().sum::<f64>() / ds.len() as f64) }        };        let Some(a) = mean(labels[i]) else { continue };        let b = groups.iter().filter(|g| **g != labels[i]).filter_map(|g| mean(*g)).fold(f64::INFINITY, f64::min);        if b.is_finite() {            total += (b - a) / a.max(b);            m += 1;        }    }    if m == 0 { None } else { Some(round(total / m as f64)) }}/// K-means: `k` groups of `points` (lists of numbers, or maps with `options.fields`). `options`: `{ fields, scale: true (standardize the columns), seed: 1, runs: 10 (the best of several starts), max_iterations: 300 }`./// Gives `{ labels (a group for each point, from 0), centres (in the points' own units), sizes, inertia, silhouette }`.pub fn cluster_kmeans(points: serde_json::Value, k: i64, options: serde_json::Value) -> Result<serde_json::Value, String> {    let (raw, names) = rows(&points, &options["fields"])?;    let k = k.clamp(1, raw.len() as i64) as usize;    let mut x = raw.clone();    let scale = options["scale"].as_bool().unwrap_or(true);    if scale {        standardize(&mut x);    }    let seed = options["seed"].as_u64().unwrap_or(1);    let runs = options["runs"].as_u64().unwrap_or(10).clamp(1, 100);    let iters = options["max_iterations"].as_u64().unwrap_or(300).clamp(1, 10_000) as usize;    let best = (0..runs).map(|r| kmeans_run(&x, k, seed + r, iters)).min_by(|a, b| a.2.partial_cmp(&b.2).unwrap()).unwrap();    let (_, labels, inertia) = best;    // Centres in the points' own units: the mean of each group's raw values.    let d = raw[0].len();    let mut centres = vec![vec![0.0; d]; k];    let mut sizes = vec![0usize; k];    for (i, p) in raw.iter().enumerate() {        sizes[labels[i]] += 1;        for j in 0..d {            centres[labels[i]][j] += p[j];        }    }    let centres: Vec<Value> = centres        .iter()        .zip(&sizes)        .map(|(c, n)| {            let v: Vec<f64> = c.iter().map(|s| round(s / (*n).max(1) as f64)).collect();            if names.is_empty() { json!(v) } else { Value::Object(names.iter().cloned().zip(v.into_iter().map(|f| json!(f))).collect()) }        })        .collect();    let li: Vec<i64> = labels.iter().map(|l| *l as i64).collect();    Ok(json!({ "labels": li, "centres": centres, "sizes": sizes, "inertia": round(inertia), "silhouette": silhouette(&x, &li) }))}/// K-means for each k from 1 to `max_k`, to choose k: `[{ k, inertia, silhouette }]` (look for the elbow, or the highest silhouette).pub fn cluster_choose_k(points: serde_json::Value, max_k: i64, options: serde_json::Value) -> Result<Vec<serde_json::Value>, String> {    let (raw, _) = rows(&points, &options["fields"])?;    let mut x = raw;    if options["scale"].as_bool().unwrap_or(true) {        standardize(&mut x);    }    let seed = options["seed"].as_u64().unwrap_or(1);    Ok((1..=max_k.clamp(1, 30).min(x.len() as i64) as usize)        .map(|k| {            let (_, labels, inertia) = (0..5).map(|r| kmeans_run(&x, k, seed + r, 300)).min_by(|a, b| a.2.partial_cmp(&b.2).unwrap()).unwrap();            let li: Vec<i64> = labels.iter().map(|l| *l as i64).collect();            json!({ "k": k, "inertia": round(inertia), "silhouette": silhouette(&x, &li) })        })        .collect())}/// DBSCAN: groups of points with at least `min_points` within `eps` of each other, any shape; points in none are noise (label -1)./// `options`: `{ fields, metric: "euclidean" | "haversine" ([lat, lon] points, `eps` in metres), scale: false }`. Gives `{ labels, groups, noise, sizes }`.pub fn cluster_dbscan(points: serde_json::Value, eps: f64, min_points: i64, options: serde_json::Value) -> Result<serde_json::Value, String> {    let (mut x, _) = rows(&points, &options["fields"])?;    if x.len() > 50_000 {        return Err("DBSCAN here takes at most 50,000 points".into());    }    let geo = options["metric"].as_str() == Some("haversine");    if geo && x[0].len() != 2 {        return Err("haversine takes [lat, lon] points".into());    }    if !geo && options["scale"] == true {        standardize(&mut x);    }    let dist = |a: &[f64], b: &[f64]| if geo { haversine(a, b) } else { d2(a, b).sqrt() };    let n = x.len();    let neighbours = |i: usize| -> Vec<usize> { (0..n).filter(|j| dist(&x[i], &x[*j]) <= eps).collect() };    let mut labels = vec![-2i64; n]; // -2 unvisited, -1 noise    let mut g = 0i64;    for i in 0..n {        if labels[i] != -2 {            continue;        }        let nb = neighbours(i);        if (nb.len() as i64) < min_points {            labels[i] = -1;            continue;        }        labels[i] = g;        let mut queue: Vec<usize> = nb;        let mut q = 0;        while q < queue.len() {            let j = queue[q];            q += 1;            if labels[j] == -1 {                labels[j] = g;            }            if labels[j] != -2 {                continue;            }            labels[j] = g;            let nb2 = neighbours(j);            if (nb2.len() as i64) >= min_points {                queue.extend(nb2);            }        }        g += 1;    }    let sizes: Vec<usize> = (0..g).map(|k| labels.iter().filter(|l| **l == k).count()).collect();    Ok(json!({ "labels": labels, "groups": g, "noise": labels.iter().filter(|l| **l == -1).count(), "sizes": sizes }))}/// Hierarchical (agglomerative) clustering, merging the nearest groups (`linkage` "single", "complete" or "average") until `groups` remain, or until the nearest are farther than `max_distance` (give 0 for none)./// Gives `{ labels, groups, merges: [{ a, b, distance, size }] }`; at most 3,000 points.pub fn cluster_hierarchical(points: serde_json::Value, groups: i64, max_distance: f64, linkage: &str, options: serde_json::Value) -> Result<serde_json::Value, String> {    let (mut x, _) = rows(&points, &options["fields"])?;    if x.len() > 3000 {        return Err("hierarchical clustering here takes at most 3,000 points".into());    }    if options["scale"] == true {        standardize(&mut x);    }    let n = x.len();    let mut members: Vec<Option<Vec<usize>>> = (0..n).map(|i| Some(vec![i])).collect();    let pd: Vec<Vec<f64>> = (0..n).map(|i| (0..n).map(|j| d2(&x[i], &x[j]).sqrt()).collect()).collect();    let link = |a: &[usize], b: &[usize]| -> f64 {        let ds = a.iter().flat_map(|i| b.iter().map(move |j| (i, j))).map(|(i, j)| pd[*i][*j]);        match linkage {            "single" => ds.fold(f64::INFINITY, f64::min),            "complete" => ds.fold(0.0, f64::max),            _ => {                let v: Vec<f64> = ds.collect();                v.iter().sum::<f64>() / v.len() as f64            }        }    };    if !["single", "complete", "average"].contains(&linkage) {        return Err(format!("linkage {:?}: single, complete or average", linkage));    }    let target = groups.max(1) as usize;    let mut merges = vec![];    let mut alive = n;    while alive > target {        let ids: Vec<usize> = (0..members.len()).filter(|i| members[*i].is_some()).collect();        let mut best = (0, 0, f64::INFINITY);        for (ai, &a) in ids.iter().enumerate() {            for &b in &ids[ai + 1..] {                let d = link(members[a].as_ref().unwrap(), members[b].as_ref().unwrap());                if d < best.2 {                    best = (a, b, d);                }            }        }        if max_distance > 0.0 && best.2 > max_distance {            break;        }        let mut m = members[best.0].take().unwrap();        m.extend(members[best.1].take().unwrap());        merges.push(json!({ "a": best.0, "b": best.1, "distance": round(best.2), "size": m.len() }));        members.push(Some(m));        alive -= 1;    }    let mut labels = vec![0i64; n];    let mut g = 0;    for m in members.iter().flatten() {        for i in m {            labels[*i] = g;        }        g += 1;    }    Ok(json!({ "labels": labels, "groups": g, "merges": merges }))}