fdars_core/

clustering.rs

//! Clustering algorithms for functional data.
//!
//! This module provides k-means and fuzzy c-means clustering algorithms
//! for functional data.

use crate::helpers::{l2_distance, simpsons_weights, NUMERICAL_EPS};
use crate::matrix::FdMatrix;
use crate::{iter_maybe_parallel, slice_maybe_parallel};
use rand::prelude::*;
#[cfg(feature = "parallel")]
use rayon::iter::ParallelIterator;

/// Result of k-means clustering.
pub struct KmeansResult {
    /// Cluster assignments for each observation
    pub cluster: Vec<usize>,
    /// Cluster centers (k x m matrix)
    pub centers: FdMatrix,
    /// Within-cluster sum of squares for each cluster
    pub withinss: Vec<f64>,
    /// Total within-cluster sum of squares
    pub tot_withinss: f64,
    /// Number of iterations
    pub iter: usize,
    /// Whether the algorithm converged
    pub converged: bool,
}

/// K-means++ initialization: select initial centers with probability proportional to D^2.
///
/// Uses incremental distance tracking: maintains a `min_dist_sq` vector and only
/// computes distances to the newest center on each iteration, avoiding redundant
/// distance computations to all existing centers.
///
/// # Arguments
/// * `curves` - Flat row-major buffer of curves (n curves, each m values)
/// * `n` - Number of curves
/// * `m` - Number of evaluation points per curve
/// * `k` - Number of clusters
/// * `weights` - Integration weights for L2 distance
/// * `rng` - Random number generator
///
/// # Returns
/// Flat buffer of k initial cluster centers (k * m values)
/// Select an index with probability proportional to the given weights.
fn weighted_random_select(dist_sq: &[f64], rng: &mut StdRng) -> usize {
    let total: f64 = dist_sq.iter().sum();
    if total < NUMERICAL_EPS {
        return rng.gen_range(0..dist_sq.len());
    }
    let r = rng.gen::<f64>() * total;
    let mut cumsum = 0.0;
    for (i, &d) in dist_sq.iter().enumerate() {
        cumsum += d;
        if cumsum >= r {
            return i;
        }
    }
    dist_sq.len() - 1
}

fn kmeans_plusplus_init(
    curves: &[f64],
    n: usize,
    m: usize,
    k: usize,
    weights: &[f64],
    rng: &mut StdRng,
) -> Vec<f64> {
    let mut centers = vec![0.0; k * m];

    // First center: random
    let first_idx = rng.gen_range(0..n);
    centers[..m].copy_from_slice(&curves[first_idx * m..(first_idx + 1) * m]);

    // Initialize min_dist_sq with squared distances to first center
    let center0 = &centers[..m];
    let mut min_dist_sq: Vec<f64> = (0..n)
        .map(|i| {
            let d = l2_distance(&curves[i * m..(i + 1) * m], center0, weights);
            d * d
        })
        .collect();

    // Remaining centers: probability proportional to D^2
    for c_idx in 1..k {
        let chosen = weighted_random_select(&min_dist_sq, rng);
        centers[c_idx * m..(c_idx + 1) * m].copy_from_slice(&curves[chosen * m..(chosen + 1) * m]);

        // Update min_dist_sq: only compute distance to the newest center
        let new_center = &centers[c_idx * m..(c_idx + 1) * m];
        for i in 0..n {
            let d_sq = l2_distance(&curves[i * m..(i + 1) * m], new_center, weights).powi(2);
            if d_sq < min_dist_sq[i] {
                min_dist_sq[i] = d_sq;
            }
        }
    }

    centers
}

/// Compute fuzzy membership values for a single observation, writing into `out`.
///
/// # Arguments
/// * `distances` - Distances from the observation to each cluster center
/// * `k` - Number of clusters
/// * `exponent` - Exponent for fuzzy membership (2 / (fuzziness - 1))
/// * `out` - Output slice (length k) to write membership values into
fn compute_fuzzy_membership_into(distances: &[f64], k: usize, exponent: f64, out: &mut [f64]) {
    out[..k].fill(0.0);

    // Check if observation is very close to any center
    for (c, &dist) in distances[..k].iter().enumerate() {
        if dist < NUMERICAL_EPS {
            // Assign full membership to this cluster
            out[c] = 1.0;
            return;
        }
    }

    // Normal fuzzy membership computation
    for c in 0..k {
        let mut sum = 0.0;
        for c2 in 0..k {
            if distances[c2] > NUMERICAL_EPS {
                sum += (distances[c] / distances[c2]).powf(exponent);
            }
        }
        out[c] = if sum > NUMERICAL_EPS { 1.0 / sum } else { 1.0 };
    }
}

/// Build an FdMatrix (k x m) from flat row-major centers buffer.
fn centers_to_matrix(centers: &[f64], k: usize, m: usize) -> FdMatrix {
    let mut flat = vec![0.0; k * m];
    for c in 0..k {
        for j in 0..m {
            flat[c + j * k] = centers[c * m + j];
        }
    }
    FdMatrix::from_column_major(flat, k, m).unwrap()
}

/// Initialize a random membership matrix (n x k) with rows summing to 1.
fn init_random_membership(n: usize, k: usize, rng: &mut StdRng) -> FdMatrix {
    let mut membership = FdMatrix::zeros(n, k);
    for i in 0..n {
        let mut row_sum = 0.0;
        for c in 0..k {
            let val = rng.gen::<f64>();
            membership[(i, c)] = val;
            row_sum += val;
        }
        for c in 0..k {
            membership[(i, c)] /= row_sum;
        }
    }
    membership
}

/// Group sample indices by their cluster assignment.
fn cluster_member_indices(cluster: &[usize], k: usize) -> Vec<Vec<usize>> {
    let mut indices = vec![Vec::new(); k];
    for (i, &c) in cluster.iter().enumerate() {
        indices[c].push(i);
    }
    indices
}

/// Assign each curve to its nearest center, returning cluster indices.
fn assign_clusters(
    curves: &[f64],
    n: usize,
    m: usize,
    centers: &[f64],
    k: usize,
    weights: &[f64],
) -> Vec<usize> {
    // Build a slice of curve slices for parallel iteration
    let curve_indices: Vec<usize> = (0..n).collect();
    slice_maybe_parallel!(curve_indices)
        .map(|&i| {
            let curve = &curves[i * m..(i + 1) * m];
            let mut best_cluster = 0;
            let mut best_dist = f64::INFINITY;
            for c in 0..k {
                let center = &centers[c * m..(c + 1) * m];
                let dist = l2_distance(curve, center, weights);
                if dist < best_dist {
                    best_dist = dist;
                    best_cluster = c;
                }
            }
            best_cluster
        })
        .collect()
}

/// Compute new cluster centers from curve assignments.
fn update_kmeans_centers(
    curves: &[f64],
    n: usize,
    m: usize,
    assignments: &[usize],
    old_centers: &[f64],
    k: usize,
) -> Vec<f64> {
    let mut centers = vec![0.0; k * m];
    let mut counts = vec![0usize; k];

    for i in 0..n {
        let c = assignments[i];
        counts[c] += 1;
        let curve = &curves[i * m..(i + 1) * m];
        let center = &mut centers[c * m..(c + 1) * m];
        for j in 0..m {
            center[j] += curve[j];
        }
    }

    for c in 0..k {
        if counts[c] > 0 {
            let center = &mut centers[c * m..(c + 1) * m];
            let n_members = counts[c] as f64;
            for j in 0..m {
                center[j] /= n_members;
            }
        } else {
            // Keep old center for empty clusters
            centers[c * m..(c + 1) * m].copy_from_slice(&old_centers[c * m..(c + 1) * m]);
        }
    }

    centers
}

/// Compute within-cluster sum of squares for each cluster.
fn compute_within_ss(
    curves: &[f64],
    n: usize,
    m: usize,
    centers: &[f64],
    assignments: &[usize],
    k: usize,
    weights: &[f64],
) -> Vec<f64> {
    let mut withinss = vec![0.0; k];
    for i in 0..n {
        let c = assignments[i];
        let dist = l2_distance(
            &curves[i * m..(i + 1) * m],
            &centers[c * m..(c + 1) * m],
            weights,
        );
        withinss[c] += dist * dist;
    }
    withinss
}

/// Update fuzzy c-means cluster centers from membership values.
fn update_fuzzy_centers(
    curves: &[f64],
    n: usize,
    m: usize,
    membership: &FdMatrix,
    k: usize,
    fuzziness: f64,
) -> Vec<f64> {
    let mut centers = vec![0.0; k * m];
    for c in 0..k {
        let mut denominator = 0.0;
        let center = &mut centers[c * m..(c + 1) * m];

        for i in 0..n {
            let weight = membership[(i, c)].powf(fuzziness);
            let curve = &curves[i * m..(i + 1) * m];
            for j in 0..m {
                center[j] += weight * curve[j];
            }
            denominator += weight;
        }

        if denominator > NUMERICAL_EPS {
            for j in 0..m {
                center[j] /= denominator;
            }
        }
    }
    centers
}

/// Update fuzzy membership values and compute max change.
fn update_fuzzy_membership_step(
    curves: &[f64],
    n: usize,
    m: usize,
    centers: &[f64],
    k: usize,
    old_membership: &FdMatrix,
    exponent: f64,
    weights: &[f64],
) -> (FdMatrix, f64) {
    let mut new_membership = FdMatrix::zeros(n, k);
    let mut max_change = 0.0;
    let mut distances = vec![0.0; k];
    let mut memberships = vec![0.0; k];

    for i in 0..n {
        let curve = &curves[i * m..(i + 1) * m];
        for c in 0..k {
            distances[c] = l2_distance(curve, &centers[c * m..(c + 1) * m], weights);
        }

        compute_fuzzy_membership_into(&distances, k, exponent, &mut memberships);

        for c in 0..k {
            new_membership[(i, c)] = memberships[c];
            let change = (memberships[c] - old_membership[(i, c)]).abs();
            if change > max_change {
                max_change = change;
            }
        }
    }

    (new_membership, max_change)
}

/// Compute mean L2 distance from a curve to a set of curve indices.
fn mean_cluster_distance(
    curve: &[f64],
    curves: &[f64],
    m: usize,
    indices: &[usize],
    weights: &[f64],
) -> f64 {
    if indices.is_empty() {
        return 0.0;
    }
    let sum: f64 = indices
        .iter()
        .map(|&j| l2_distance(curve, &curves[j * m..(j + 1) * m], weights))
        .sum();
    sum / indices.len() as f64
}

/// Compute cluster centers, global mean, and counts from curves and assignments.
fn compute_centers_and_global_mean(
    curves: &[f64],
    n: usize,
    m: usize,
    assignments: &[usize],
    k: usize,
) -> (Vec<f64>, Vec<f64>, Vec<usize>) {
    let mut global_mean = vec![0.0; m];
    for i in 0..n {
        let curve = &curves[i * m..(i + 1) * m];
        for j in 0..m {
            global_mean[j] += curve[j];
        }
    }
    for j in 0..m {
        global_mean[j] /= n as f64;
    }

    let mut centers = vec![0.0; k * m];
    let mut counts = vec![0usize; k];
    for i in 0..n {
        let c = assignments[i];
        counts[c] += 1;
        let curve = &curves[i * m..(i + 1) * m];
        let center = &mut centers[c * m..(c + 1) * m];
        for j in 0..m {
            center[j] += curve[j];
        }
    }
    for c in 0..k {
        if counts[c] > 0 {
            let center = &mut centers[c * m..(c + 1) * m];
            for j in 0..m {
                center[j] /= counts[c] as f64;
            }
        }
    }

    (centers, global_mean, counts)
}

/// Run one k-means iteration: assign clusters, update centers, compute movement.
fn kmeans_step(
    curves: &[f64],
    n: usize,
    m: usize,
    centers: &[f64],
    k: usize,
    weights: &[f64],
) -> (Vec<usize>, Vec<f64>, f64) {
    let new_cluster = assign_clusters(curves, n, m, centers, k, weights);
    let new_centers = update_kmeans_centers(curves, n, m, &new_cluster, centers, k);
    let max_movement = (0..k)
        .map(|c| {
            l2_distance(
                &centers[c * m..(c + 1) * m],
                &new_centers[c * m..(c + 1) * m],
                weights,
            )
        })
        .fold(0.0, f64::max);
    (new_cluster, new_centers, max_movement)
}

/// Run the k-means iteration loop until convergence or max iterations.
fn kmeans_iterate(
    curves: &[f64],
    n: usize,
    m: usize,
    mut centers: Vec<f64>,
    k: usize,
    weights: &[f64],
    max_iter: usize,
    tol: f64,
) -> (Vec<usize>, Vec<f64>, usize, bool) {
    let mut cluster = vec![0usize; n];
    let mut converged = false;
    let mut iter = 0;

    for iteration in 0..max_iter {
        iter = iteration + 1;
        let (new_cluster, new_centers, max_movement) =
            kmeans_step(curves, n, m, &centers, k, weights);

        if new_cluster == cluster {
            converged = true;
            break;
        }
        cluster = new_cluster;
        centers = new_centers;

        if max_movement < tol {
            converged = true;
            break;
        }
    }

    (cluster, centers, iter, converged)
}

/// K-means clustering for functional data.
///
/// # Arguments
/// * `data` - Functional data matrix (n x m)
/// * `argvals` - Evaluation points
/// * `k` - Number of clusters
/// * `max_iter` - Maximum iterations
/// * `tol` - Convergence tolerance
/// * `seed` - Random seed
pub fn kmeans_fd(
    data: &FdMatrix,
    argvals: &[f64],
    k: usize,
    max_iter: usize,
    tol: f64,
    seed: u64,
) -> KmeansResult {
    let n = data.nrows();
    let m = data.ncols();

    if n == 0 || m == 0 || k == 0 || k > n || argvals.len() != m {
        return KmeansResult {
            cluster: Vec::new(),
            centers: FdMatrix::zeros(0, 0),
            withinss: Vec::new(),
            tot_withinss: 0.0,
            iter: 0,
            converged: false,
        };
    }

    let weights = simpsons_weights(argvals);
    let mut rng = StdRng::seed_from_u64(seed);

    // Extract curves as flat row-major buffer
    let curves = data.to_row_major();

    // K-means++ initialization
    let centers = kmeans_plusplus_init(&curves, n, m, k, &weights, &mut rng);

    let (cluster, centers, iter, converged) =
        kmeans_iterate(&curves, n, m, centers, k, &weights, max_iter, tol);

    let withinss = compute_within_ss(&curves, n, m, &centers, &cluster, k, &weights);
    let tot_withinss: f64 = withinss.iter().sum();
    let centers_mat = centers_to_matrix(&centers, k, m);

    KmeansResult {
        cluster,
        centers: centers_mat,
        withinss,
        tot_withinss,
        iter,
        converged,
    }
}

/// Result of fuzzy c-means clustering.
pub struct FuzzyCmeansResult {
    /// Membership matrix (n x k)
    pub membership: FdMatrix,
    /// Cluster centers (k x m)
    pub centers: FdMatrix,
    /// Number of iterations
    pub iter: usize,
    /// Whether the algorithm converged
    pub converged: bool,
}

/// Fuzzy c-means clustering for functional data.
///
/// # Arguments
/// * `data` - Functional data matrix (n x m)
/// * `argvals` - Evaluation points
/// * `k` - Number of clusters
/// * `fuzziness` - Fuzziness parameter (> 1)
/// * `max_iter` - Maximum iterations
/// * `tol` - Convergence tolerance
/// * `seed` - Random seed
pub fn fuzzy_cmeans_fd(
    data: &FdMatrix,
    argvals: &[f64],
    k: usize,
    fuzziness: f64,
    max_iter: usize,
    tol: f64,
    seed: u64,
) -> FuzzyCmeansResult {
    let n = data.nrows();
    let m = data.ncols();

    if n == 0 || m == 0 || k == 0 || k > n || argvals.len() != m || fuzziness <= 1.0 {
        return FuzzyCmeansResult {
            membership: FdMatrix::zeros(0, 0),
            centers: FdMatrix::zeros(0, 0),
            iter: 0,
            converged: false,
        };
    }

    let weights = simpsons_weights(argvals);
    let mut rng = StdRng::seed_from_u64(seed);

    // Extract curves as flat row-major buffer
    let curves = data.to_row_major();

    let mut membership = init_random_membership(n, k, &mut rng);

    let mut centers = vec![0.0; k * m];
    let mut converged = false;
    let mut iter = 0;
    let exponent = 2.0 / (fuzziness - 1.0);

    for iteration in 0..max_iter {
        iter = iteration + 1;

        centers = update_fuzzy_centers(&curves, n, m, &membership, k, fuzziness);

        let (new_membership, max_change) = update_fuzzy_membership_step(
            &curves,
            n,
            m,
            &centers,
            k,
            &membership,
            exponent,
            &weights,
        );

        membership = new_membership;

        if max_change < tol {
            converged = true;
            break;
        }
    }

    let centers_mat = centers_to_matrix(&centers, k, m);

    FuzzyCmeansResult {
        membership,
        centers: centers_mat,
        iter,
        converged,
    }
}

/// Compute silhouette score for clustering result.
pub fn silhouette_score(data: &FdMatrix, argvals: &[f64], cluster: &[usize]) -> Vec<f64> {
    let n = data.nrows();
    let m = data.ncols();

    if n == 0 || m == 0 || cluster.len() != n || argvals.len() != m {
        return Vec::new();
    }

    let weights = simpsons_weights(argvals);
    let curves = data.to_row_major();

    let k = cluster.iter().cloned().max().unwrap_or(0) + 1;
    let members = cluster_member_indices(cluster, k);

    iter_maybe_parallel!(0..n)
        .map(|i| {
            let my_cluster = cluster[i];
            let curve_i = &curves[i * m..(i + 1) * m];

            let same_indices: Vec<usize> = members[my_cluster]
                .iter()
                .copied()
                .filter(|&j| j != i)
                .collect();
            let a_i = mean_cluster_distance(curve_i, &curves, m, &same_indices, &weights);

            let mut b_i = f64::INFINITY;
            for c in 0..k {
                if c != my_cluster && !members[c].is_empty() {
                    b_i = b_i.min(mean_cluster_distance(
                        curve_i,
                        &curves,
                        m,
                        &members[c],
                        &weights,
                    ));
                }
            }

            if b_i.is_infinite() {
                0.0
            } else {
                let max_ab = a_i.max(b_i);
                if max_ab > NUMERICAL_EPS {
                    (b_i - a_i) / max_ab
                } else {
                    0.0
                }
            }
        })
        .collect()
}

/// Compute Calinski-Harabasz index for clustering result.
pub fn calinski_harabasz(data: &FdMatrix, argvals: &[f64], cluster: &[usize]) -> f64 {
    let n = data.nrows();
    let m = data.ncols();

    if n == 0 || m == 0 || cluster.len() != n || argvals.len() != m {
        return 0.0;
    }

    let weights = simpsons_weights(argvals);
    let curves = data.to_row_major();

    let k = cluster.iter().cloned().max().unwrap_or(0) + 1;
    if k < 2 {
        return 0.0;
    }

    let (centers, global_mean, counts) = compute_centers_and_global_mean(&curves, n, m, cluster, k);

    let mut bgss = 0.0;
    for c in 0..k {
        let dist = l2_distance(&centers[c * m..(c + 1) * m], &global_mean, &weights);
        bgss += counts[c] as f64 * dist * dist;
    }

    let wgss_vec = compute_within_ss(&curves, n, m, &centers, cluster, k, &weights);
    let wgss: f64 = wgss_vec.iter().sum();

    if wgss < NUMERICAL_EPS {
        return f64::INFINITY;
    }

    (bgss / (k - 1) as f64) / (wgss / (n - k) as f64)
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::f64::consts::PI;

    /// Generate a uniform grid of points
    fn uniform_grid(n: usize) -> Vec<f64> {
        (0..n).map(|i| i as f64 / (n - 1) as f64).collect()
    }

    /// Generate two clearly separated clusters of curves as an FdMatrix
    fn generate_two_clusters(n_per_cluster: usize, m: usize) -> (FdMatrix, Vec<f64>) {
        let t = uniform_grid(m);
        let n = 2 * n_per_cluster;
        let mut col_major = vec![0.0; n * m];

        // Cluster 0: sine waves with low amplitude
        for i in 0..n_per_cluster {
            for (j, &ti) in t.iter().enumerate() {
                col_major[i + j * n] =
                    (2.0 * PI * ti).sin() + 0.1 * (i as f64 / n_per_cluster as f64);
            }
        }

        // Cluster 1: sine waves shifted up by 5
        for i in 0..n_per_cluster {
            for (j, &ti) in t.iter().enumerate() {
                col_major[(i + n_per_cluster) + j * n] =
                    (2.0 * PI * ti).sin() + 5.0 + 0.1 * (i as f64 / n_per_cluster as f64);
            }
        }

        (FdMatrix::from_column_major(col_major, n, m).unwrap(), t)
    }

    // ============== K-means tests ==============

    #[test]
    fn test_kmeans_fd_basic() {
        let m = 50;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);

        assert_eq!(result.cluster.len(), n);
        assert!(result.converged);
        assert!(result.iter > 0 && result.iter <= 100);
    }

    #[test]
    fn test_kmeans_fd_finds_clusters() {
        let m = 50;
        let n_per = 10;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);

        // First half should be one cluster, second half the other
        let cluster_0 = result.cluster[0];
        let cluster_1 = result.cluster[n_per];

        assert_ne!(cluster_0, cluster_1, "Clusters should be different");

        // Check that first half is in same cluster
        for i in 0..n_per {
            assert_eq!(result.cluster[i], cluster_0);
        }

        // Check that second half is in same cluster
        for i in n_per..n {
            assert_eq!(result.cluster[i], cluster_1);
        }
    }

    #[test]
    fn test_kmeans_fd_deterministic() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);

        let result1 = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);
        let result2 = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);

        // Same seed should give same results
        assert_eq!(result1.cluster, result2.cluster);
    }

    #[test]
    fn test_kmeans_fd_withinss() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);

        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);

        // Within-cluster sum of squares should be non-negative
        for &wss in &result.withinss {
            assert!(wss >= 0.0);
        }

        // Total should equal sum
        let sum: f64 = result.withinss.iter().sum();
        assert!((sum - result.tot_withinss).abs() < 1e-10);
    }

    #[test]
    fn test_kmeans_fd_centers_shape() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let k = 3;

        let result = kmeans_fd(&data, &t, k, 100, 1e-6, 42);

        // Centers should be k x m matrix
        assert_eq!(result.centers.nrows(), k);
        assert_eq!(result.centers.ncols(), m);
    }

    #[test]
    fn test_kmeans_fd_invalid_input() {
        let t = uniform_grid(30);

        // Empty data
        let data = FdMatrix::zeros(0, 0);
        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);
        assert!(result.cluster.is_empty());
        assert!(!result.converged);

        // k > n
        let data = FdMatrix::zeros(5, 30);
        let result = kmeans_fd(&data, &t, 10, 100, 1e-6, 42);
        assert!(result.cluster.is_empty());
    }

    #[test]
    fn test_kmeans_fd_single_cluster() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 10;
        let data = FdMatrix::zeros(n, m);

        let result = kmeans_fd(&data, &t, 1, 100, 1e-6, 42);

        // All should be in cluster 0
        for &c in &result.cluster {
            assert_eq!(c, 0);
        }
    }

    // ============== Fuzzy C-means tests ==============

    #[test]
    fn test_fuzzy_cmeans_fd_basic() {
        let m = 50;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let result = fuzzy_cmeans_fd(&data, &t, 2, 2.0, 100, 1e-6, 42);

        assert_eq!(result.membership.nrows(), n);
        assert_eq!(result.membership.ncols(), 2);
        assert!(result.iter > 0);
    }

    #[test]
    fn test_fuzzy_cmeans_fd_membership_sums_to_one() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;
        let k = 2;

        let result = fuzzy_cmeans_fd(&data, &t, k, 2.0, 100, 1e-6, 42);

        // Each observation's membership should sum to 1
        for i in 0..n {
            let sum: f64 = (0..k).map(|c| result.membership[(i, c)]).sum();
            assert!(
                (sum - 1.0).abs() < 1e-6,
                "Membership should sum to 1, got {}",
                sum
            );
        }
    }

    #[test]
    fn test_fuzzy_cmeans_fd_membership_in_range() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);

        let result = fuzzy_cmeans_fd(&data, &t, 2, 2.0, 100, 1e-6, 42);

        // All memberships should be in [0, 1]
        for &mem in result.membership.as_slice() {
            assert!((0.0..=1.0 + 1e-10).contains(&mem));
        }
    }

    #[test]
    fn test_fuzzy_cmeans_fd_fuzziness_effect() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);

        let result_low = fuzzy_cmeans_fd(&data, &t, 2, 1.5, 100, 1e-6, 42);
        let result_high = fuzzy_cmeans_fd(&data, &t, 2, 3.0, 100, 1e-6, 42);

        // Higher fuzziness should give more diffuse memberships
        // Measure by entropy-like metric
        let entropy_low: f64 = result_low
            .membership
            .as_slice()
            .iter()
            .map(|&m| if m > 1e-10 { -m * m.ln() } else { 0.0 })
            .sum();

        let entropy_high: f64 = result_high
            .membership
            .as_slice()
            .iter()
            .map(|&m| if m > 1e-10 { -m * m.ln() } else { 0.0 })
            .sum();

        assert!(
            entropy_high >= entropy_low - 0.1,
            "Higher fuzziness should give higher entropy"
        );
    }

    #[test]
    fn test_fuzzy_cmeans_fd_invalid_fuzziness() {
        let t = uniform_grid(30);
        let data = FdMatrix::zeros(10, 30);

        // Fuzziness <= 1 should fail
        let result = fuzzy_cmeans_fd(&data, &t, 2, 1.0, 100, 1e-6, 42);
        assert!(result.membership.is_empty());

        let result = fuzzy_cmeans_fd(&data, &t, 2, 0.5, 100, 1e-6, 42);
        assert!(result.membership.is_empty());
    }

    #[test]
    fn test_fuzzy_cmeans_fd_centers_shape() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 10;
        let k = 3;
        let data = FdMatrix::zeros(n, m);

        let result = fuzzy_cmeans_fd(&data, &t, k, 2.0, 100, 1e-6, 42);

        assert_eq!(result.centers.nrows(), k);
        assert_eq!(result.centers.ncols(), m);
    }

    // ============== Silhouette score tests ==============

    #[test]
    fn test_silhouette_score_well_separated() {
        let m = 30;
        let n_per = 10;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        // Perfect clustering: first half in 0, second in 1
        let cluster: Vec<usize> = (0..n).map(|i| if i < n_per { 0 } else { 1 }).collect();

        let scores = silhouette_score(&data, &t, &cluster);

        assert_eq!(scores.len(), n);

        // Well-separated clusters should have high silhouette scores
        let mean_score: f64 = scores.iter().sum::<f64>() / n as f64;
        assert!(
            mean_score > 0.5,
            "Well-separated clusters should have high silhouette: {}",
            mean_score
        );
    }

    #[test]
    fn test_silhouette_score_range() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let cluster: Vec<usize> = (0..n).map(|i| if i < n_per { 0 } else { 1 }).collect();

        let scores = silhouette_score(&data, &t, &cluster);

        // Silhouette scores should be in [-1, 1]
        for &s in &scores {
            assert!((-1.0 - 1e-10..=1.0 + 1e-10).contains(&s));
        }
    }

    #[test]
    fn test_silhouette_score_single_cluster() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 10;
        let data = FdMatrix::zeros(n, m);

        // All in one cluster
        let cluster = vec![0usize; n];

        let scores = silhouette_score(&data, &t, &cluster);

        // Single cluster should give zeros
        for &s in &scores {
            assert!(s.abs() < 1e-10);
        }
    }

    #[test]
    fn test_silhouette_score_invalid_input() {
        let t = uniform_grid(30);

        // Empty data
        let data = FdMatrix::zeros(0, 0);
        let scores = silhouette_score(&data, &t, &[]);
        assert!(scores.is_empty());

        // Mismatched cluster length
        let data = FdMatrix::zeros(10, 30);
        let cluster = vec![0; 5]; // Wrong length
        let scores = silhouette_score(&data, &t, &cluster);
        assert!(scores.is_empty());
    }

    // ============== Calinski-Harabasz tests ==============

    #[test]
    fn test_calinski_harabasz_well_separated() {
        let m = 30;
        let n_per = 10;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let cluster: Vec<usize> = (0..n).map(|i| if i < n_per { 0 } else { 1 }).collect();

        let ch = calinski_harabasz(&data, &t, &cluster);

        // Well-separated clusters should have high CH index
        assert!(
            ch > 1.0,
            "Well-separated clusters should have high CH: {}",
            ch
        );
    }

    #[test]
    fn test_calinski_harabasz_positive() {
        let m = 30;
        let n_per = 5;
        let (data, t) = generate_two_clusters(n_per, m);
        let n = 2 * n_per;

        let cluster: Vec<usize> = (0..n).map(|i| if i < n_per { 0 } else { 1 }).collect();

        let ch = calinski_harabasz(&data, &t, &cluster);

        assert!(ch >= 0.0, "CH index should be non-negative");
    }

    #[test]
    fn test_calinski_harabasz_single_cluster() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 10;
        let data = FdMatrix::zeros(n, m);

        // All in one cluster
        let cluster = vec![0usize; n];

        let ch = calinski_harabasz(&data, &t, &cluster);

        // Single cluster should give 0
        assert!(ch.abs() < 1e-10);
    }

    #[test]
    fn test_calinski_harabasz_invalid_input() {
        let t = uniform_grid(30);

        // Empty data
        let data = FdMatrix::zeros(0, 0);
        let ch = calinski_harabasz(&data, &t, &[]);
        assert!(ch.abs() < 1e-10);
    }

    #[test]
    fn test_identical_curves_kmeans() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 10;
        // All curves identical
        let data_vec: Vec<f64> = (0..n * m)
            .map(|idx| (2.0 * PI * (idx % m) as f64 / (m - 1) as f64).sin())
            .collect();
        let data = FdMatrix::from_column_major(data_vec, n, m).unwrap();
        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);
        // Should not panic with identical data
        assert_eq!(result.cluster.len(), n);
    }

    #[test]
    fn test_k_equals_n() {
        let m = 30;
        let t = uniform_grid(m);
        let n = 5;
        let (data, _) = generate_two_clusters(n, m);
        let result = kmeans_fd(&data, &t, 2 * n, 100, 1e-6, 42);
        // k == n: each curve is its own cluster
        assert_eq!(result.cluster.len(), 2 * n);
    }

    #[test]
    fn test_n2_kmeans() {
        let m = 30;
        let t = uniform_grid(m);
        let mut data = FdMatrix::zeros(2, m);
        for j in 0..m {
            data[(0, j)] = 0.0;
            data[(1, j)] = 10.0;
        }
        let result = kmeans_fd(&data, &t, 2, 100, 1e-6, 42);
        assert_eq!(result.cluster.len(), 2);
        assert_ne!(result.cluster[0], result.cluster[1]);
    }
}