face-core 0.1.0

//! §5.2 Natural breaks (Fisher–Jenks) clustering on the axis field's
//! numeric value.
//!
//! - Same record-extraction policy as Bands and Quantiles: missing
//!   path drops silently, present-but-non-numeric drops with a
//!   [`SkipReport`].
//! - Sorts the resolved values ascending and runs the classic
//!   Fisher–Jenks one-dimensional optimal-break algorithm with `k =
//!   count`. Naive O(N²·k) implementation is fine for face's expected
//!   sizes (≤10K items per the brief). No external dependency.
//! - Output: one cluster per Jenks class, sorted numeric descending by
//!   class lower bound so highest-scoring ranges appear first. Class
//!   bounds use the actual record min/max
//!   inside the class (mirroring Quantiles — for Natural the "edges"
//!   are determined by the data, not a synthetic interval).
//! - `count == 0` is treated as `1` (degenerate single class). When
//!   `count >= N`, every distinct sorted position becomes its own
//!   class (the algorithm degenerates to one record per class).
//! - Labels: en-dash, integer-vs-fractional heuristic from Bands.

use serde_json::json;

use crate::{
    Cluster, ClusterId, Diagnostics, FaceError, Record, Strategy,
    cluster_tree::{
        AxisPlan, build_tree_impl, extend_id,
        util::{resolve_numeric, score_range},
    },
};

/// Cluster `items` by Jenks natural breaks on the resolved numeric
/// value at `plan.axis.field`, recursing per cluster into `plan.within`.
pub(super) fn cluster<D: Diagnostics + ?Sized>(
    plan: &AxisPlan,
    items: Vec<Record>,
    parent: &ClusterId,
    diag: &mut D,
) -> Result<Vec<Cluster>, FaceError> {
    let Strategy::Natural { count } = &plan.axis.strategy else {
        unreachable!("natural clusterer is only called for Strategy::Natural")
    };
    let count = (*count).max(1) as usize;

    let mut resolved = resolve_numeric(items, &plan.axis.field, diag);
    if resolved.is_empty() {
        return Ok(Vec::new());
    }

    // Sort ascending — Jenks is defined on a sorted vector.
    resolved.sort_by(|a, b| a.0.total_cmp(&b.0));
    let n = resolved.len();

    if count == 1 || n == 1 {
        let lo = resolved[0].0;
        let hi = resolved[n - 1].0;
        let group: Vec<Record> = resolved.into_iter().map(|(_, r)| r).collect();
        return Ok(vec![build_cluster(
            plan,
            parent,
            format_label(lo, hi, &[lo, hi]),
            lo,
            hi,
            group,
            diag,
        )?]);
    }

    // Cap k at n: with more requested classes than records, each
    // record becomes its own class.
    let k = count.min(n);

    // Compute Jenks breaks on the sorted values. `breaks[i]` is the
    // sorted index of the *first* record in class `i`. Always
    // `breaks[0] == 0` and `breaks[k] == n`.
    let values: Vec<f64> = resolved.iter().map(|(v, _)| *v).collect();
    let breaks = jenks_breaks(&values, k);

    let mut out = Vec::with_capacity(k);
    let mut iter = resolved.into_iter();
    for class in 0..k {
        let start = breaks[class];
        let end = breaks[class + 1];
        if start >= end {
            // Defensive: a degenerate breakpoint that produces an
            // empty class. Shouldn't happen for the standard
            // algorithm but if it did we'd skip rather than emit
            // a zero-total cluster.
            continue;
        }
        let len = end - start;
        let group: Vec<Record> = (&mut iter).take(len).map(|(_, r)| r).collect();
        debug_assert_eq!(group.len(), len, "iterator drained mid-class");
        let lo = values[start];
        let hi = values[end - 1];
        let label = format_label(lo, hi, &[lo, hi]);
        out.push(build_cluster(plan, parent, label, lo, hi, group, diag)?);
    }
    out.reverse();
    Ok(out)
}

/// Compute Fisher–Jenks natural breaks on a sorted slice `values`,
/// partitioning into `k` classes. Returns `k+1` indices: `breaks[0] == 0`,
/// `breaks[k] == n`, and `breaks[i]` is the (sorted) index of the first
/// record in class `i`. The classes are `[breaks[i], breaks[i+1])`.
///
/// The algorithm minimizes the sum of within-class variance (equiv.
/// sum-of-squared-deviations from the class mean) over all partitions.
/// We build:
///
/// - `lower_class_limits[i][m]` (1-indexed in classic refs; here
///   we use 1-based in the matrices for parity with the published
///   pseudocode and convert to 0-based at the boundary).
/// - `variance_combinations[i][m]` — the optimal SSD when the first
///   `i` records are split into `m` classes.
///
/// Cost: O(N²·k) time, O(N·k) memory. Fine for ≤10K records (the
/// face brief's expected upper bound).
fn jenks_breaks(values: &[f64], k: usize) -> Vec<usize> {
    let n = values.len();
    debug_assert!(n > 0 && k > 0 && k <= n);

    if k == 1 {
        return vec![0, n];
    }
    if k == n {
        // One record per class.
        return (0..=n).collect();
    }

    // 1-based dimensions (size n+1 by k+1) so the indexing matches
    // the published pseudocode and the recurrence. The 0th row /
    // column carry the base cases.
    let big = f64::MAX / 2.0;
    let mut lower_class_limits = vec![vec![0usize; k + 1]; n + 1];
    let mut variance_combinations = vec![vec![big; k + 1]; n + 1];

    for j in 1..=k {
        lower_class_limits[1][j] = 1;
        variance_combinations[1][j] = 0.0;
    }

    for l in 2..=n {
        // Running statistics for the inner window [m..=l] (1-based).
        let mut sum = 0.0f64;
        let mut sum_squares = 0.0f64;
        let mut w = 0.0f64;
        for m in 1..=l {
            // Lower class limit walks from l down to 1.
            let lower_class_limit = l + 1 - m;
            let val = values[lower_class_limit - 1];
            w += 1.0;
            sum += val;
            sum_squares += val * val;
            // SSD of [lower_class_limit..=l].
            let variance = sum_squares - (sum * sum) / w;
            let i4 = lower_class_limit.saturating_sub(1);
            if i4 != 0 {
                for j in 2..=k {
                    let candidate = variance + variance_combinations[i4][j - 1];
                    if variance_combinations[l][j] >= candidate {
                        lower_class_limits[l][j] = lower_class_limit;
                        variance_combinations[l][j] = candidate;
                    }
                }
            }
        }
        lower_class_limits[l][1] = 1;
        variance_combinations[l][1] = sum_squares - (sum * sum) / w;
    }

    // Backtrack: `lower_class_limits[k][n]` (1-based) is the start
    // of the last class. Walk down through the `j` dimension.
    let mut breaks = vec![0usize; k + 1];
    breaks[k] = n;
    let mut k_idx = n;
    for j in (2..=k).rev() {
        // 1-based start of class `j` is `lower_class_limits[k_idx][j]`.
        let start_one_based = lower_class_limits[k_idx][j];
        // 0-based start of class `j-1`'s end / class `j`'s start is
        // `start_one_based - 1`.
        breaks[j - 1] = start_one_based - 1;
        // For the next iteration up, the upper boundary of class
        // `j-1` is its start - 1 (still 1-based).
        k_idx = start_one_based - 1;
    }
    breaks[0] = 0;
    breaks
}

/// Build one cluster from a (label, bounds, records) triple, recursing
/// into `plan.within` if any.
fn build_cluster<D: Diagnostics + ?Sized>(
    plan: &AxisPlan,
    parent: &ClusterId,
    label: String,
    lo: f64,
    hi: f64,
    group: Vec<Record>,
    diag: &mut D,
) -> Result<Cluster, FaceError> {
    let id = extend_id(parent, &plan.axis.field, &label);
    let total = group.len() as u64;
    let (score_min, score_max) = score_range(&group);

    let children = if plan.within.is_empty() {
        Vec::new()
    } else {
        let mut acc: Vec<Cluster> = Vec::new();
        for child_plan in &plan.within {
            let child_clusters = build_tree_impl(child_plan, group.clone(), &id, diag)?;
            acc.extend(child_clusters);
        }
        acc
    };

    Ok(Cluster {
        id,
        label,
        axis: plan.axis.field.clone(),
        value: json!({"min": lo, "max": hi}),
        total,
        score_min,
        score_max,
        clusters: children,
    })
}

/// Render `lo–hi` with the en-dash (U+2013). Same heuristic as Bands /
/// Quantiles: integer-valued bounds → no decimals; otherwise two-decimal
/// precision.
fn format_label(lo: f64, hi: f64, bounds: &[f64]) -> String {
    if bounds.iter().all(|b| b.fract() == 0.0) {
        format!("{}\u{2013}{}", lo as i64, hi as i64)
    } else {
        format!("{lo:.2}\u{2013}{hi:.2}")
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{Axis, NullDiagnostics, Strategy};
    use serde_json::json;

    fn axis(field: &str, count: u8) -> Axis {
        Axis {
            field: field.into(),
            strategy: Strategy::Natural { count },
            auto: false,
        }
    }

    #[test]
    fn breaks_clustered_input() {
        // Two well-separated clusters; Jenks with k=2 should split
        // exactly between them.
        let plan = AxisPlan::leaf(axis("v", 2));
        let items = vec![
            json!({"v": 1.0}),
            json!({"v": 1.1}),
            json!({"v": 1.2}),
            json!({"v": 9.0}),
            json!({"v": 9.1}),
            json!({"v": 9.2}),
        ];
        let mut diag = NullDiagnostics;
        let clusters = super::cluster(
            &plan,
            Record::from_items(items, None),
            &ClusterId::default(),
            &mut diag,
        )
        .unwrap();
        assert_eq!(clusters.len(), 2);
        // Each class carries 3 records (the two natural groups).
        assert_eq!(clusters[0].total, 3);
        assert_eq!(clusters[1].total, 3);
        // Descending by lower bound.
        assert_eq!(clusters[0].value, json!({"min": 9.0, "max": 9.2}));
        assert_eq!(clusters[1].value, json!({"min": 1.0, "max": 1.2}));
    }

    #[test]
    fn single_record() {
        let plan = AxisPlan::leaf(axis("v", 4));
        let items = vec![json!({"v": 0.42})];
        let mut diag = NullDiagnostics;
        let clusters = super::cluster(
            &plan,
            Record::from_items(items, None),
            &ClusterId::default(),
            &mut diag,
        )
        .unwrap();
        assert_eq!(clusters.len(), 1);
        assert_eq!(clusters[0].total, 1);
        assert_eq!(clusters[0].value, json!({"min": 0.42, "max": 0.42}));
    }

    #[test]
    fn empty_records_returns_empty() {
        let plan = AxisPlan::leaf(axis("v", 4));
        let items: Vec<serde_json::Value> = vec![];
        let mut diag = NullDiagnostics;
        let clusters = super::cluster(
            &plan,
            Record::from_items(items, None),
            &ClusterId::default(),
            &mut diag,
        )
        .unwrap();
        assert!(clusters.is_empty());
    }

    #[test]
    fn count_one_returns_single_class() {
        let plan = AxisPlan::leaf(axis("v", 1));
        let items = vec![json!({"v": 1.0}), json!({"v": 5.0}), json!({"v": 9.0})];
        let mut diag = NullDiagnostics;
        let clusters = super::cluster(
            &plan,
            Record::from_items(items, None),
            &ClusterId::default(),
            &mut diag,
        )
        .unwrap();
        assert_eq!(clusters.len(), 1);
        assert_eq!(clusters[0].total, 3);
        assert_eq!(clusters[0].value, json!({"min": 1.0, "max": 9.0}));
    }
}