sklears-multioutput 0.1.1

//! Multi-output and multi-label evaluation metrics

// Use SciRS2-Core for arrays and random number generation (SciRS2 Policy)
use scirs2_core::ndarray::ArrayView2;
use sklears_core::error::{Result as SklResult, SklearsError};

/// Hamming loss for multi-label classification
///
/// The Hamming loss is the fraction of the wrong labels to the total
/// number of labels. It is a multi-label generalization of the zero-one loss.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// The Hamming loss between y_true and y_pred
pub fn hamming_loss(y_true: &ArrayView2<'_, i32>, y_pred: &ArrayView2<'_, i32>) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_errors = 0;
    let total_elements = n_samples * n_labels;

    for sample_idx in 0..n_samples {
        for label_idx in 0..n_labels {
            if y_true[[sample_idx, label_idx]] != y_pred[[sample_idx, label_idx]] {
                total_errors += 1;
            }
        }
    }

    Ok(total_errors as f64 / total_elements as f64)
}

/// Subset accuracy for multi-label classification
///
/// Subset accuracy is the most strict metric. It requires for each sample
/// that each label set be correctly predicted.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// The subset accuracy between y_true and y_pred
pub fn subset_accuracy(
    y_true: &ArrayView2<'_, i32>,
    y_pred: &ArrayView2<'_, i32>,
) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample".to_string(),
        ));
    }

    let mut correct_subsets = 0;

    for sample_idx in 0..n_samples {
        let mut subset_correct = true;
        for label_idx in 0..n_labels {
            if y_true[[sample_idx, label_idx]] != y_pred[[sample_idx, label_idx]] {
                subset_correct = false;
                break;
            }
        }
        if subset_correct {
            correct_subsets += 1;
        }
    }

    Ok(correct_subsets as f64 / n_samples as f64)
}

/// Jaccard similarity coefficient for multi-label classification
///
/// The Jaccard similarity coefficient is defined as the size of the intersection
/// divided by the size of the union of the sample sets.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// The Jaccard similarity coefficient
pub fn jaccard_score(y_true: &ArrayView2<'_, i32>, y_pred: &ArrayView2<'_, i32>) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample".to_string(),
        ));
    }

    let mut total_jaccard = 0.0;

    for sample_idx in 0..n_samples {
        let mut intersection = 0;
        let mut union = 0;

        for label_idx in 0..n_labels {
            let true_label = y_true[[sample_idx, label_idx]];
            let pred_label = y_pred[[sample_idx, label_idx]];

            if true_label == 1 && pred_label == 1 {
                intersection += 1;
            }
            if true_label == 1 || pred_label == 1 {
                union += 1;
            }
        }

        // Jaccard = intersection / union, handle division by zero
        let sample_jaccard = if union > 0 {
            intersection as f64 / union as f64
        } else {
            1.0 // If both sets are empty, Jaccard = 1
        };

        total_jaccard += sample_jaccard;
    }

    Ok(total_jaccard / n_samples as f64)
}

/// F1 score for multi-label classification
///
/// Compute the F1 score for each label and return the specified average.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
/// * `average` - The averaging strategy ('micro', 'macro', 'samples')
///
/// # Returns
///
/// The F1 score according to the specified averaging strategy
pub fn f1_score(
    y_true: &ArrayView2<'_, i32>,
    y_pred: &ArrayView2<'_, i32>,
    average: &str,
) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    match average {
        "micro" => {
            // Compute global precision and recall
            let mut total_tp = 0;
            let mut total_fp = 0;
            let mut total_false_negatives = 0;

            for sample_idx in 0..n_samples {
                for label_idx in 0..n_labels {
                    let true_label = y_true[[sample_idx, label_idx]];
                    let pred_label = y_pred[[sample_idx, label_idx]];

                    if true_label == 1 && pred_label == 1 {
                        total_tp += 1;
                    } else if true_label == 0 && pred_label == 1 {
                        total_fp += 1;
                    } else if true_label == 1 && pred_label == 0 {
                        total_false_negatives += 1;
                    }
                }
            }

            let precision = if total_tp + total_fp > 0 {
                total_tp as f64 / (total_tp + total_fp) as f64
            } else {
                0.0
            };

            let recall = if total_tp + total_false_negatives > 0 {
                total_tp as f64 / (total_tp + total_false_negatives) as f64
            } else {
                0.0
            };

            let f1 = if precision + recall > 0.0 {
                2.0 * precision * recall / (precision + recall)
            } else {
                0.0
            };

            Ok(f1)
        }
        "macro" => {
            // Compute F1 for each label and average
            let mut label_f1_scores = Vec::new();

            for label_idx in 0..n_labels {
                let mut tp = 0;
                let mut fp = 0;
                let mut false_negatives = 0;

                for sample_idx in 0..n_samples {
                    let true_label = y_true[[sample_idx, label_idx]];
                    let pred_label = y_pred[[sample_idx, label_idx]];

                    if true_label == 1 && pred_label == 1 {
                        tp += 1;
                    } else if true_label == 0 && pred_label == 1 {
                        fp += 1;
                    } else if true_label == 1 && pred_label == 0 {
                        false_negatives += 1;
                    }
                }

                let precision = if tp + fp > 0 {
                    tp as f64 / (tp + fp) as f64
                } else {
                    0.0
                };

                let recall = if tp + false_negatives > 0 {
                    tp as f64 / (tp + false_negatives) as f64
                } else {
                    0.0
                };

                let f1 = if precision + recall > 0.0 {
                    2.0 * precision * recall / (precision + recall)
                } else {
                    0.0
                };

                label_f1_scores.push(f1);
            }

            Ok(label_f1_scores.iter().sum::<f64>() / n_labels as f64)
        }
        "samples" => {
            // Compute F1 for each sample and average
            let mut sample_f1_scores = Vec::new();

            for sample_idx in 0..n_samples {
                let mut tp = 0;
                let mut fp = 0;
                let mut false_negatives = 0;

                for label_idx in 0..n_labels {
                    let true_label = y_true[[sample_idx, label_idx]];
                    let pred_label = y_pred[[sample_idx, label_idx]];

                    if true_label == 1 && pred_label == 1 {
                        tp += 1;
                    } else if true_label == 0 && pred_label == 1 {
                        fp += 1;
                    } else if true_label == 1 && pred_label == 0 {
                        false_negatives += 1;
                    }
                }

                let precision = if tp + fp > 0 {
                    tp as f64 / (tp + fp) as f64
                } else {
                    0.0
                };

                let recall = if tp + false_negatives > 0 {
                    tp as f64 / (tp + false_negatives) as f64
                } else {
                    0.0
                };

                let f1 = if precision + recall > 0.0 {
                    2.0 * precision * recall / (precision + recall)
                } else {
                    0.0
                };

                sample_f1_scores.push(f1);
            }

            Ok(sample_f1_scores.iter().sum::<f64>() / n_samples as f64)
        }
        _ => Err(SklearsError::InvalidInput(format!(
            "Unknown average type: {}. Valid options are 'micro', 'macro', 'samples'",
            average
        ))),
    }
}

/// Coverage error for multi-label ranking
///
/// Coverage error measures how far we need to go through the ranked scores
/// to cover all true labels. The best value is equal to the average number
/// of labels in y_true per sample.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_scores` - Target scores (predicted probabilities)
///
/// # Returns
///
/// The coverage error
pub fn coverage_error(
    y_true: &ArrayView2<'_, i32>,
    y_scores: &ArrayView2<'_, f64>,
) -> SklResult<f64> {
    if y_true.dim() != y_scores.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_scores must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_coverage = 0.0;

    for sample_idx in 0..n_samples {
        // Get the indices sorted by scores in descending order
        let mut score_label_pairs: Vec<(f64, usize)> = (0..n_labels)
            .map(|label_idx| (y_scores[[sample_idx, label_idx]], label_idx))
            .collect();
        score_label_pairs
            .sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));

        // Find the position of the last true label in the ranked list
        let mut last_true_position = 0;
        for (position, &(_, label_idx)) in score_label_pairs.iter().enumerate() {
            if y_true[[sample_idx, label_idx]] == 1 {
                last_true_position = position + 1; // Convert to 1-based indexing
            }
        }

        total_coverage += last_true_position as f64;
    }

    Ok(total_coverage / n_samples as f64)
}

/// Label ranking average precision for multi-label ranking
///
/// The label ranking average precision (LRAP) averages over the samples
/// the answer to the following question: for each ground truth label,
/// what fraction of higher-ranked labels were true labels?
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_scores` - Target scores (predicted probabilities)
///
/// # Returns
///
/// The label ranking average precision
pub fn label_ranking_average_precision(
    y_true: &ArrayView2<'_, i32>,
    y_scores: &ArrayView2<'_, f64>,
) -> SklResult<f64> {
    if y_true.dim() != y_scores.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_scores must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_lrap = 0.0;

    for sample_idx in 0..n_samples {
        // Get the indices sorted by scores in descending order
        let mut score_label_pairs: Vec<(f64, usize)> = (0..n_labels)
            .map(|label_idx| (y_scores[[sample_idx, label_idx]], label_idx))
            .collect();
        score_label_pairs
            .sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));

        // Count true labels for this sample
        let n_true_labels: i32 = (0..n_labels)
            .map(|label_idx| y_true[[sample_idx, label_idx]])
            .sum();

        if n_true_labels == 0 {
            continue; // Skip samples with no true labels
        }

        let mut precision_sum = 0.0;
        let mut true_labels_seen = 0;

        for (position, &(_, label_idx)) in score_label_pairs.iter().enumerate() {
            if y_true[[sample_idx, label_idx]] == 1 {
                true_labels_seen += 1;
                let precision_at_position = true_labels_seen as f64 / (position + 1) as f64;
                precision_sum += precision_at_position;
            }
        }

        let sample_lrap = precision_sum / n_true_labels as f64;
        total_lrap += sample_lrap;
    }

    Ok(total_lrap / n_samples as f64)
}

/// One-error for multi-label ranking
///
/// The one-error evaluates how many times the top-ranked label is not
/// in the set of true labels. The best performance is achieved when
/// one-error is 0, which means the top-ranked label is always correct.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_scores` - Target scores (predicted probabilities)
///
/// # Returns
///
/// The one-error (fraction of samples where top-ranked label is incorrect)
pub fn one_error(y_true: &ArrayView2<'_, i32>, y_scores: &ArrayView2<'_, f64>) -> SklResult<f64> {
    if y_true.dim() != y_scores.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_scores must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut errors = 0;

    for sample_idx in 0..n_samples {
        // Find the label with the highest score
        let mut max_score = f64::NEG_INFINITY;
        let mut top_label_idx = 0;

        for label_idx in 0..n_labels {
            let score = y_scores[[sample_idx, label_idx]];
            if score > max_score {
                max_score = score;
                top_label_idx = label_idx;
            }
        }

        // Check if the top-ranked label is correct
        if y_true[[sample_idx, top_label_idx]] != 1 {
            errors += 1;
        }
    }

    Ok(errors as f64 / n_samples as f64)
}

/// Ranking loss for multi-label ranking
///
/// The ranking loss evaluates the average fraction of label pairs that are
/// incorrectly ordered, given the predictions. The best performance is achieved
/// when ranking loss is 0.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_scores` - Target scores (predicted probabilities)
///
/// # Returns
///
/// The ranking loss
pub fn ranking_loss(
    y_true: &ArrayView2<'_, i32>,
    y_scores: &ArrayView2<'_, f64>,
) -> SklResult<f64> {
    if y_true.dim() != y_scores.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_scores must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_ranking_loss = 0.0;

    for sample_idx in 0..n_samples {
        let mut incorrect_pairs = 0;
        let mut total_pairs = 0;

        // Compare all pairs of labels
        for i in 0..n_labels {
            for j in 0..n_labels {
                if i != j {
                    let true_i = y_true[[sample_idx, i]];
                    let true_j = y_true[[sample_idx, j]];
                    let score_i = y_scores[[sample_idx, i]];
                    let score_j = y_scores[[sample_idx, j]];

                    // Check if this is a relevant pair (one positive, one negative)
                    if (true_i == 1 && true_j == 0) || (true_i == 0 && true_j == 1) {
                        total_pairs += 1;

                        // Check if the ordering is incorrect
                        if (true_i == 1 && true_j == 0 && score_i < score_j)
                            || (true_i == 0 && true_j == 1 && score_i > score_j)
                        {
                            incorrect_pairs += 1;
                        }
                    }
                }
            }
        }

        // Add to total ranking loss
        if total_pairs > 0 {
            total_ranking_loss += incorrect_pairs as f64 / total_pairs as f64;
        }
    }

    Ok(total_ranking_loss / n_samples as f64)
}

/// Average precision score for multi-label ranking
///
/// Computes the average precision for each sample and then averages
/// over all samples. This is different from label ranking average precision
/// as it focuses on precision-recall curves for each sample.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_scores` - Target scores (predicted probabilities)
///
/// # Returns
///
/// The average precision score
pub fn average_precision_score(
    y_true: &ArrayView2<'_, i32>,
    y_scores: &ArrayView2<'_, f64>,
) -> SklResult<f64> {
    if y_true.dim() != y_scores.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_scores must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_ap = 0.0;

    for sample_idx in 0..n_samples {
        // Get the indices sorted by scores in descending order
        let mut score_label_pairs: Vec<(f64, usize)> = (0..n_labels)
            .map(|label_idx| (y_scores[[sample_idx, label_idx]], label_idx))
            .collect();
        score_label_pairs
            .sort_by(|a, b| b.0.partial_cmp(&a.0).unwrap_or(std::cmp::Ordering::Equal));

        // Count true labels for this sample
        let n_true_labels: i32 = (0..n_labels)
            .map(|label_idx| y_true[[sample_idx, label_idx]])
            .sum();

        if n_true_labels == 0 {
            continue; // Skip samples with no true labels
        }

        let mut precision_sum = 0.0;
        let mut true_labels_seen = 0;

        for (position, &(_, label_idx)) in score_label_pairs.iter().enumerate() {
            if y_true[[sample_idx, label_idx]] == 1 {
                true_labels_seen += 1;
                let precision_at_position = true_labels_seen as f64 / (position + 1) as f64;
                precision_sum += precision_at_position;
            }
        }

        let sample_ap = precision_sum / n_true_labels as f64;
        total_ap += sample_ap;
    }

    Ok(total_ap / n_samples as f64)
}

/// Micro-averaged precision for multi-label classification
///
/// Calculate precision by counting the total true positives and false positives
/// across all labels and samples.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// The micro-averaged precision
pub fn precision_score_micro(
    y_true: &ArrayView2<'_, i32>,
    y_pred: &ArrayView2<'_, i32>,
) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_tp = 0;
    let mut total_fp = 0;

    for sample_idx in 0..n_samples {
        for label_idx in 0..n_labels {
            let true_label = y_true[[sample_idx, label_idx]];
            let pred_label = y_pred[[sample_idx, label_idx]];

            if true_label == 1 && pred_label == 1 {
                total_tp += 1;
            } else if true_label == 0 && pred_label == 1 {
                total_fp += 1;
            }
        }
    }

    let precision = if total_tp + total_fp > 0 {
        total_tp as f64 / (total_tp + total_fp) as f64
    } else {
        0.0
    };

    Ok(precision)
}

/// Micro-averaged recall for multi-label classification
///
/// Calculate recall by counting the total true positives and false negatives
/// across all labels and samples.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// The micro-averaged recall
pub fn recall_score_micro(
    y_true: &ArrayView2<'_, i32>,
    y_pred: &ArrayView2<'_, i32>,
) -> SklResult<f64> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut total_tp = 0;
    let mut total_fn = 0;

    for sample_idx in 0..n_samples {
        for label_idx in 0..n_labels {
            let true_label = y_true[[sample_idx, label_idx]];
            let pred_label = y_pred[[sample_idx, label_idx]];

            if true_label == 1 && pred_label == 1 {
                total_tp += 1;
            } else if true_label == 1 && pred_label == 0 {
                total_fn += 1;
            }
        }
    }

    let recall = if total_tp + total_fn > 0 {
        total_tp as f64 / (total_tp + total_fn) as f64
    } else {
        0.0
    };

    Ok(recall)
}

// Additional imports for statistical tests
use std::collections::HashMap;

/// Per-label performance metrics
///
/// Container for detailed performance metrics for each label, including
/// precision, recall, F1-score, support, and accuracy per label.
#[derive(Debug, Clone)]
pub struct PerLabelMetrics {
    /// Precision score for each label
    pub precision: Vec<f64>,
    /// Recall score for each label
    pub recall: Vec<f64>,
    /// F1 score for each label
    pub f1_score: Vec<f64>,
    /// Support (number of true instances) for each label
    pub support: Vec<usize>,
    /// Accuracy for each label (considering label as binary classification)
    pub accuracy: Vec<f64>,
    /// Number of labels
    pub n_labels: usize,
}

impl PerLabelMetrics {
    /// Get the macro average of a metric
    pub fn macro_average(&self, metric: &str) -> SklResult<f64> {
        let values = match metric {
            "precision" => &self.precision,
            "recall" => &self.recall,
            "f1_score" => &self.f1_score,
            "accuracy" => &self.accuracy,
            _ => return Err(SklearsError::InvalidInput(format!(
                "Unknown metric: {}. Valid options are 'precision', 'recall', 'f1_score', 'accuracy'",
                metric
            )))
        };

        Ok(values.iter().sum::<f64>() / values.len() as f64)
    }

    /// Get the weighted average of a metric (weighted by support)
    pub fn weighted_average(&self, metric: &str) -> SklResult<f64> {
        let values = match metric {
            "precision" => &self.precision,
            "recall" => &self.recall,
            "f1_score" => &self.f1_score,
            "accuracy" => &self.accuracy,
            _ => return Err(SklearsError::InvalidInput(format!(
                "Unknown metric: {}. Valid options are 'precision', 'recall', 'f1_score', 'accuracy'",
                metric
            )))
        };

        let total_support: usize = self.support.iter().sum();
        if total_support == 0 {
            return Ok(0.0);
        }

        let weighted_sum: f64 = values
            .iter()
            .zip(self.support.iter())
            .map(|(value, support)| value * (*support as f64))
            .sum();

        Ok(weighted_sum / total_support as f64)
    }
}

/// Compute detailed per-label performance metrics
///
/// Calculates precision, recall, F1-score, support, and accuracy for each label
/// individually, providing comprehensive per-label analysis.
///
/// # Arguments
///
/// * `y_true` - Ground truth (correct) labels
/// * `y_pred` - Predicted labels
///
/// # Returns
///
/// PerLabelMetrics containing detailed metrics for each label
pub fn per_label_metrics(
    y_true: &ArrayView2<'_, i32>,
    y_pred: &ArrayView2<'_, i32>,
) -> SklResult<PerLabelMetrics> {
    if y_true.dim() != y_pred.dim() {
        return Err(SklearsError::InvalidInput(
            "y_true and y_pred must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    let mut precision = Vec::with_capacity(n_labels);
    let mut recall = Vec::with_capacity(n_labels);
    let mut f1_score = Vec::with_capacity(n_labels);
    let mut support = Vec::with_capacity(n_labels);
    let mut accuracy = Vec::with_capacity(n_labels);

    // Calculate metrics for each label
    for label_idx in 0..n_labels {
        let mut tp = 0;
        let mut fp = 0;
        let mut fn_count = 0;
        let mut tn = 0;

        for sample_idx in 0..n_samples {
            let true_label = y_true[[sample_idx, label_idx]];
            let pred_label = y_pred[[sample_idx, label_idx]];

            match (true_label, pred_label) {
                (1, 1) => tp += 1,
                (0, 1) => fp += 1,
                (1, 0) => fn_count += 1,
                (0, 0) => tn += 1,
                _ => {} // Should not happen with binary labels
            }
        }

        // Calculate precision
        let label_precision = if tp + fp > 0 {
            tp as f64 / (tp + fp) as f64
        } else {
            0.0
        };

        // Calculate recall
        let label_recall = if tp + fn_count > 0 {
            tp as f64 / (tp + fn_count) as f64
        } else {
            0.0
        };

        // Calculate F1 score
        let label_f1 = if label_precision + label_recall > 0.0 {
            2.0 * label_precision * label_recall / (label_precision + label_recall)
        } else {
            0.0
        };

        // Calculate accuracy (for this label as binary classification)
        let label_accuracy = (tp + tn) as f64 / n_samples as f64;

        // Support is number of true instances for this label
        let label_support = (tp + fn_count) as usize;

        precision.push(label_precision);
        recall.push(label_recall);
        f1_score.push(label_f1);
        support.push(label_support);
        accuracy.push(label_accuracy);
    }

    Ok(PerLabelMetrics {
        precision,
        recall,
        f1_score,
        support,
        accuracy,
        n_labels,
    })
}

/// Statistical significance test result
#[derive(Debug, Clone)]
pub struct StatisticalTestResult {
    /// Test statistic value
    pub statistic: f64,
    /// P-value of the test
    pub p_value: f64,
    /// Whether the result is statistically significant (p < 0.05)
    pub is_significant: bool,
    /// Test name
    pub test_name: String,
    /// Additional information about the test
    pub additional_info: HashMap<String, f64>,
}

impl StatisticalTestResult {
    /// Create a new statistical test result
    pub fn new(
        statistic: f64,
        p_value: f64,
        test_name: String,
        additional_info: Option<HashMap<String, f64>>,
    ) -> Self {
        Self {
            statistic,
            p_value,
            is_significant: p_value < 0.05,
            test_name,
            additional_info: additional_info.unwrap_or_default(),
        }
    }
}

/// McNemar's test for comparing two classifiers
///
/// Tests whether two classifiers have significantly different error rates.
/// Appropriate for comparing paired predictions on the same test set.
///
/// # Arguments
///
/// * `y_true` - Ground truth labels
/// * `y_pred1` - Predictions from first classifier
/// * `y_pred2` - Predictions from second classifier
///
/// # Returns
///
/// Statistical test result with McNemar's test statistic and p-value
pub fn mcnemar_test(
    y_true: &ArrayView2<'_, i32>,
    y_pred1: &ArrayView2<'_, i32>,
    y_pred2: &ArrayView2<'_, i32>,
) -> SklResult<StatisticalTestResult> {
    if y_true.dim() != y_pred1.dim() || y_true.dim() != y_pred2.dim() {
        return Err(SklearsError::InvalidInput(
            "All input arrays must have the same shape".to_string(),
        ));
    }

    let (n_samples, n_labels) = y_true.dim();
    if n_samples == 0 || n_labels == 0 {
        return Err(SklearsError::InvalidInput(
            "Input arrays must have at least one sample and one label".to_string(),
        ));
    }

    // Count disagreements across all samples and labels
    let mut n01 = 0; // Classifier 1 correct, classifier 2 incorrect
    let mut n10 = 0; // Classifier 1 incorrect, classifier 2 correct

    for sample_idx in 0..n_samples {
        for label_idx in 0..n_labels {
            let true_label = y_true[[sample_idx, label_idx]];
            let pred1 = y_pred1[[sample_idx, label_idx]];
            let pred2 = y_pred2[[sample_idx, label_idx]];

            let correct1 = pred1 == true_label;
            let correct2 = pred2 == true_label;

            match (correct1, correct2) {
                (true, false) => n01 += 1,
                (false, true) => n10 += 1,
                _ => {} // Both correct or both incorrect
            }
        }
    }

    // McNemar's test statistic
    let total_disagreements = n01 + n10;
    if total_disagreements == 0 {
        return Ok(StatisticalTestResult::new(
            0.0,
            1.0, // No disagreements means no significant difference
            "McNemar".to_string(),
            Some({
                let mut info = HashMap::new();
                info.insert("n01".to_string(), n01 as f64);
                info.insert("n10".to_string(), n10 as f64);
                info.insert(
                    "total_disagreements".to_string(),
                    total_disagreements as f64,
                );
                info
            }),
        ));
    }

    // Use continuity correction for McNemar's test
    let statistic = ((n01 as f64 - n10 as f64).abs() - 1.0).max(0.0).powi(2) / (n01 + n10) as f64;

    // Chi-square distribution approximation for p-value (1 degree of freedom)
    let p_value = chi_square_p_value(statistic, 1);

    let mut info = HashMap::new();
    info.insert("n01".to_string(), n01 as f64);
    info.insert("n10".to_string(), n10 as f64);
    info.insert(
        "total_disagreements".to_string(),
        total_disagreements as f64,
    );

    Ok(StatisticalTestResult::new(
        statistic,
        p_value,
        "McNemar".to_string(),
        Some(info),
    ))
}

/// Paired t-test for metric comparisons
///
/// Tests whether the mean difference between paired metric values is
/// significantly different from zero.
///
/// # Arguments
///
/// * `metric_values1` - Metric values from first method
/// * `metric_values2` - Metric values from second method
///
/// # Returns
///
/// Statistical test result with t-statistic and p-value
pub fn paired_t_test(
    metric_values1: &[f64],
    metric_values2: &[f64],
) -> SklResult<StatisticalTestResult> {
    if metric_values1.len() != metric_values2.len() {
        return Err(SklearsError::InvalidInput(
            "Metric value arrays must have the same length".to_string(),
        ));
    }

    let n = metric_values1.len();
    if n < 2 {
        return Err(SklearsError::InvalidInput(
            "Need at least 2 paired observations for t-test".to_string(),
        ));
    }

    // Calculate differences
    let differences: Vec<f64> = metric_values1
        .iter()
        .zip(metric_values2.iter())
        .map(|(v1, v2)| v1 - v2)
        .collect();

    // Calculate mean difference
    let mean_diff = differences.iter().sum::<f64>() / n as f64;

    // Calculate standard deviation of differences
    let variance = differences
        .iter()
        .map(|d| (d - mean_diff).powi(2))
        .sum::<f64>()
        / (n - 1) as f64;
    let std_dev = variance.sqrt();

    // t-statistic
    let t_statistic = mean_diff / (std_dev / (n as f64).sqrt());

    // Degrees of freedom
    let df = n - 1;

    // Two-tailed p-value using t-distribution approximation
    let p_value = 2.0 * (1.0 - t_distribution_cdf(t_statistic.abs(), df as f64));

    let mut info = HashMap::new();
    info.insert("mean_difference".to_string(), mean_diff);
    info.insert("std_dev_diff".to_string(), std_dev);
    info.insert("degrees_of_freedom".to_string(), df as f64);
    info.insert("n_observations".to_string(), n as f64);

    Ok(StatisticalTestResult::new(
        t_statistic,
        p_value,
        "Paired t-test".to_string(),
        Some(info),
    ))
}

/// Wilcoxon signed-rank test for non-parametric metric comparison
///
/// Non-parametric alternative to paired t-test that doesn't assume
/// normal distribution of differences.
///
/// # Arguments
///
/// * `metric_values1` - Metric values from first method
/// * `metric_values2` - Metric values from second method
///
/// # Returns
///
/// Statistical test result with Wilcoxon statistic and approximate p-value
pub fn wilcoxon_signed_rank_test(
    metric_values1: &[f64],
    metric_values2: &[f64],
) -> SklResult<StatisticalTestResult> {
    if metric_values1.len() != metric_values2.len() {
        return Err(SklearsError::InvalidInput(
            "Metric value arrays must have the same length".to_string(),
        ));
    }

    let n = metric_values1.len();
    if n < 3 {
        return Err(SklearsError::InvalidInput(
            "Need at least 3 paired observations for Wilcoxon signed-rank test".to_string(),
        ));
    }

    // Calculate differences and their absolute values
    let mut differences_with_abs: Vec<(f64, f64, bool)> = metric_values1
        .iter()
        .zip(metric_values2.iter())
        .map(|(v1, v2)| {
            let diff = v1 - v2;
            (diff, diff.abs(), diff > 0.0)
        })
        .filter(|(_, abs_diff, _)| *abs_diff > 1e-10) // Remove ties (zero differences)
        .collect();

    let n_nonzero = differences_with_abs.len();
    if n_nonzero < 3 {
        return Err(SklearsError::InvalidInput(
            "Too many zero differences for Wilcoxon test".to_string(),
        ));
    }

    // Sort by absolute difference to assign ranks
    differences_with_abs.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("operation should succeed"));

    // Assign ranks (handling ties by averaging)
    let mut ranks = vec![0.0; n_nonzero];
    let mut i = 0;
    while i < n_nonzero {
        let current_abs_diff = differences_with_abs[i].1;
        let mut j = i;

        // Find end of tie group
        while j < n_nonzero && (differences_with_abs[j].1 - current_abs_diff).abs() < 1e-10 {
            j += 1;
        }

        // Assign average rank to tied values
        let avg_rank = (i + j + 1) as f64 / 2.0;
        for rank in ranks.iter_mut().take(j).skip(i) {
            *rank = avg_rank;
        }

        i = j;
    }

    // Calculate positive and negative rank sums
    let mut w_plus = 0.0;
    let mut w_minus = 0.0;

    for i in 0..n_nonzero {
        if differences_with_abs[i].2 {
            // Positive difference
            w_plus += ranks[i];
        } else {
            // Negative difference
            w_minus += ranks[i];
        }
    }

    // Test statistic is the smaller of the two rank sums
    let w_statistic = w_plus.min(w_minus);

    // Normal approximation for p-value (valid for n >= 10)
    let expected_w = (n_nonzero * (n_nonzero + 1)) as f64 / 4.0;
    let variance_w = (n_nonzero * (n_nonzero + 1) * (2 * n_nonzero + 1)) as f64 / 24.0;
    let std_w = variance_w.sqrt();

    // Continuity correction
    let z_score = ((w_statistic - expected_w).abs() - 0.5) / std_w;

    // Two-tailed p-value
    let p_value = 2.0 * (1.0 - standard_normal_cdf(z_score));

    let mut info = HashMap::new();
    info.insert("w_plus".to_string(), w_plus);
    info.insert("w_minus".to_string(), w_minus);
    info.insert("n_nonzero_differences".to_string(), n_nonzero as f64);
    info.insert("z_score".to_string(), z_score);

    Ok(StatisticalTestResult::new(
        w_statistic,
        p_value,
        "Wilcoxon signed-rank".to_string(),
        Some(info),
    ))
}

/// Confidence interval for a metric
#[derive(Debug, Clone)]
pub struct ConfidenceInterval {
    /// Lower bound of confidence interval
    pub lower: f64,
    /// Upper bound of confidence interval
    pub upper: f64,
    /// Point estimate (mean)
    pub point_estimate: f64,
    /// Confidence level (e.g., 0.95 for 95%)
    pub confidence_level: f64,
}

/// Calculate confidence interval for metric values
///
/// Computes confidence interval assuming normal distribution of metric values.
///
/// # Arguments
///
/// * `metric_values` - Array of metric values
/// * `confidence_level` - Confidence level (e.g., 0.95 for 95% CI)
///
/// # Returns
///
/// Confidence interval with lower and upper bounds
pub fn confidence_interval(
    metric_values: &[f64],
    confidence_level: f64,
) -> SklResult<ConfidenceInterval> {
    if metric_values.is_empty() {
        return Err(SklearsError::InvalidInput(
            "Metric values array cannot be empty".to_string(),
        ));
    }

    if confidence_level <= 0.0 || confidence_level >= 1.0 {
        return Err(SklearsError::InvalidInput(
            "Confidence level must be between 0 and 1".to_string(),
        ));
    }

    let n = metric_values.len();
    let mean = metric_values.iter().sum::<f64>() / n as f64;

    if n == 1 {
        return Ok(ConfidenceInterval {
            lower: mean,
            upper: mean,
            point_estimate: mean,
            confidence_level,
        });
    }

    // Calculate standard error
    let variance = metric_values
        .iter()
        .map(|v| (v - mean).powi(2))
        .sum::<f64>()
        / (n - 1) as f64;
    let std_error = (variance / n as f64).sqrt();

    // Critical value for t-distribution
    let alpha = 1.0 - confidence_level;
    let df = (n - 1) as f64;
    let t_critical = t_distribution_quantile(1.0 - alpha / 2.0, df);

    // Margin of error
    let margin_error = t_critical * std_error;

    Ok(ConfidenceInterval {
        lower: mean - margin_error,
        upper: mean + margin_error,
        point_estimate: mean,
        confidence_level,
    })
}

// Statistical distribution helper functions

/// Chi-square p-value approximation (1 degree of freedom)
fn chi_square_p_value(x: f64, df: usize) -> f64 {
    if df == 1 {
        // For 1 df, chi-square is the square of standard normal
        2.0 * (1.0 - standard_normal_cdf(x.sqrt()))
    } else {
        // Simplified approximation for other degrees of freedom
        let normalized = (x - df as f64) / (2.0 * df as f64).sqrt();
        2.0 * (1.0 - standard_normal_cdf(normalized.abs()))
    }
}

/// Standard normal CDF approximation
fn standard_normal_cdf(z: f64) -> f64 {
    0.5 * (1.0 + erf(z / 2.0_f64.sqrt()))
}

/// Error function approximation
fn erf(x: f64) -> f64 {
    // Abramowitz and Stegun approximation
    let a1 = 0.254829592;
    let a2 = -0.284496736;
    let a3 = 1.421413741;
    let a4 = -1.453152027;
    let a5 = 1.061405429;
    let p = 0.3275911;

    let sign = if x >= 0.0 { 1.0 } else { -1.0 };
    let x = x.abs();

    let t = 1.0 / (1.0 + p * x);
    let y = 1.0 - ((((a5 * t + a4) * t + a3) * t + a2) * t + a1) * t * (-x * x).exp();

    sign * y
}

/// t-distribution CDF approximation
fn t_distribution_cdf(t: f64, df: f64) -> f64 {
    if df > 30.0 {
        // For large df, t-distribution approaches standard normal
        standard_normal_cdf(t)
    } else {
        // Simplified approximation
        let normalized = t / (df + t * t).sqrt();
        0.5 + 0.5 * erf(normalized)
    }
}

/// t-distribution quantile approximation
fn t_distribution_quantile(p: f64, df: f64) -> f64 {
    if df > 100.0 {
        // For very large df, use normal quantile approximation
        normal_quantile(p)
    } else if df >= 2.0 {
        // Use Wilson-Hilferty approximation for better accuracy
        let z = normal_quantile(p);
        let h = 2.0 / (9.0 * df);
        let correction = z.powi(2) * h / 6.0;
        z * (1.0 + correction).max(0.1) // Ensure positive correction
    } else {
        // For very small df, use simpler approximation
        let z = normal_quantile(p);
        z * (1.0 + (z.powi(2) + 1.0) / (4.0 * df))
    }
}

/// Standard normal quantile approximation (inverse CDF) - Simple Box-Muller inspired approach
fn normal_quantile(p: f64) -> f64 {
    if p <= 0.0 {
        return f64::NEG_INFINITY;
    }
    if p >= 1.0 {
        return f64::INFINITY;
    }
    if (p - 0.5).abs() < f64::EPSILON {
        return 0.0;
    }

    // Use a lookup table approach for key values and interpolation for others
    let known_values = [
        (0.001, -3.090232),
        (0.005, -2.575829),
        (0.01, -2.326348),
        (0.025, -1.959964),
        (0.05, -1.644854),
        (0.1, -1.281552),
        (0.15, -1.036433),
        (0.2, -0.841621),
        (0.25, -0.674490),
        (0.3, -0.524401),
        (0.35, -0.385320),
        (0.4, -0.253347),
        (0.45, -0.125661),
        (0.5, 0.0),
        (0.55, 0.125661),
        (0.6, 0.253347),
        (0.65, 0.385320),
        (0.7, 0.524401),
        (0.75, 0.674490),
        (0.8, 0.841621),
        (0.85, 1.036433),
        (0.9, 1.281552),
        (0.95, 1.644854),
        (0.975, 1.959964),
        (0.99, 2.326348),
        (0.995, 2.575829),
        (0.999, 3.090232),
    ];

    // Find the closest tabulated values and interpolate
    if let Some(idx) = known_values.iter().position(|(prob, _)| *prob >= p) {
        if idx == 0 {
            return known_values[0].1;
        }

        let (p1, z1) = known_values[idx - 1];
        let (p2, z2) = known_values[idx];

        // Linear interpolation
        let weight = (p - p1) / (p2 - p1);
        z1 + weight * (z2 - z1)
    } else {
        // For very high probabilities, extrapolate
        3.5 // Conservative upper bound
    }
}

#[allow(non_snake_case)]
#[cfg(test)]
mod tests {
    use super::*;
    use scirs2_core::ndarray::{array, Array2};

    // Helper function to create test data
    fn create_test_data() -> (Array2<i32>, Array2<i32>) {
        let y_true = array![[1, 0, 1], [0, 1, 0], [1, 1, 1], [0, 0, 0], [1, 0, 1]];
        let y_pred = array![[1, 0, 0], [0, 1, 1], [1, 0, 1], [1, 0, 0], [1, 1, 1]];
        (y_true, y_pred)
    }

    #[test]
    fn test_per_label_metrics_basic() {
        let (y_true, y_pred) = create_test_data();
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        let metrics =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");

        assert_eq!(metrics.n_labels, 3);
        assert_eq!(metrics.precision.len(), 3);
        assert_eq!(metrics.recall.len(), 3);
        assert_eq!(metrics.f1_score.len(), 3);
        assert_eq!(metrics.support.len(), 3);
        assert_eq!(metrics.accuracy.len(), 3);

        // Check support (number of true instances per label)
        assert_eq!(metrics.support[0], 3); // Label 0 has 3 true instances
        assert_eq!(metrics.support[1], 2); // Label 1 has 2 true instances
        assert_eq!(metrics.support[2], 3); // Label 2 has 3 true instances

        // Verify all metrics are within valid range [0, 1]
        for i in 0..3 {
            assert!(metrics.precision[i] >= 0.0 && metrics.precision[i] <= 1.0);
            assert!(metrics.recall[i] >= 0.0 && metrics.recall[i] <= 1.0);
            assert!(metrics.f1_score[i] >= 0.0 && metrics.f1_score[i] <= 1.0);
            assert!(metrics.accuracy[i] >= 0.0 && metrics.accuracy[i] <= 1.0);
        }
    }

    #[test]
    fn test_per_label_metrics_perfect_prediction() {
        let y_perfect = array![[1, 0, 1], [0, 1, 0], [1, 1, 1]];
        let y_true_view = y_perfect.view();
        let y_pred_view = y_perfect.view(); // Same as true

        let metrics =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");

        // All metrics should be 1.0 for perfect prediction
        for i in 0..3 {
            assert!((metrics.precision[i] - 1.0).abs() < 1e-10);
            assert!((metrics.recall[i] - 1.0).abs() < 1e-10);
            assert!((metrics.f1_score[i] - 1.0).abs() < 1e-10);
            assert!((metrics.accuracy[i] - 1.0).abs() < 1e-10);
        }
    }

    #[test]
    fn test_per_label_metrics_all_zeros() {
        let y_true = array![[0, 0, 0], [0, 0, 0], [0, 0, 0]];
        let y_pred = array![[0, 0, 0], [0, 0, 0], [0, 0, 0]];
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        let metrics =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");

        // All support should be 0
        for i in 0..3 {
            assert_eq!(metrics.support[i], 0);
            assert!((metrics.accuracy[i] - 1.0).abs() < 1e-10); // All correct (all negative)
            assert!((metrics.precision[i] - 0.0).abs() < 1e-10); // No positive predictions
            assert!((metrics.recall[i] - 0.0).abs() < 1e-10); // No positive instances
        }
    }

    #[test]
    fn test_per_label_metrics_macro_average() {
        let (y_true, y_pred) = create_test_data();
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        let metrics =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");

        let macro_precision = metrics
            .macro_average("precision")
            .expect("operation should succeed");
        let macro_recall = metrics
            .macro_average("recall")
            .expect("operation should succeed");
        let macro_f1 = metrics
            .macro_average("f1_score")
            .expect("operation should succeed");
        let macro_accuracy = metrics
            .macro_average("accuracy")
            .expect("operation should succeed");

        // Verify macro averages are calculated correctly
        let expected_precision = metrics.precision.iter().sum::<f64>() / 3.0;
        let expected_recall = metrics.recall.iter().sum::<f64>() / 3.0;
        let expected_f1 = metrics.f1_score.iter().sum::<f64>() / 3.0;
        let expected_accuracy = metrics.accuracy.iter().sum::<f64>() / 3.0;

        assert!((macro_precision - expected_precision).abs() < 1e-10);
        assert!((macro_recall - expected_recall).abs() < 1e-10);
        assert!((macro_f1 - expected_f1).abs() < 1e-10);
        assert!((macro_accuracy - expected_accuracy).abs() < 1e-10);
    }

    #[test]
    fn test_per_label_metrics_weighted_average() {
        let (y_true, y_pred) = create_test_data();
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        let metrics =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");

        let weighted_precision = metrics
            .weighted_average("precision")
            .expect("operation should succeed");
        let weighted_recall = metrics
            .weighted_average("recall")
            .expect("operation should succeed");
        let weighted_f1 = metrics
            .weighted_average("f1_score")
            .expect("operation should succeed");
        let weighted_accuracy = metrics
            .weighted_average("accuracy")
            .expect("operation should succeed");

        // All should be valid values
        assert!((0.0..=1.0).contains(&weighted_precision));
        assert!((0.0..=1.0).contains(&weighted_recall));
        assert!((0.0..=1.0).contains(&weighted_f1));
        assert!((0.0..=1.0).contains(&weighted_accuracy));

        // Test invalid metric name
        assert!(metrics.weighted_average("invalid").is_err());
        assert!(metrics.macro_average("invalid").is_err());
    }

    #[test]
    fn test_per_label_metrics_error_handling() {
        let y_true = array![[1, 0], [0, 1]];
        let y_pred = array![[1, 0, 1], [0, 1, 0]]; // Different shape

        let result = per_label_metrics(&y_true.view(), &y_pred.view());
        assert!(result.is_err());

        // Empty arrays
        let empty_true = Array2::<i32>::zeros((0, 0));
        let empty_pred = Array2::<i32>::zeros((0, 0));
        let result = per_label_metrics(&empty_true.view(), &empty_pred.view());
        assert!(result.is_err());
    }

    #[test]
    fn test_mcnemar_test_identical_classifiers() {
        let y_true = array![[1, 0, 1], [0, 1, 0], [1, 1, 1]];
        let y_pred = array![[1, 0, 0], [0, 1, 1], [1, 0, 1]];

        // Test identical classifiers (should have p-value = 1.0)
        let result = mcnemar_test(&y_true.view(), &y_pred.view(), &y_pred.view())
            .expect("operation should succeed");
        assert_eq!(result.test_name, "McNemar");
        assert!((result.p_value - 1.0).abs() < 1e-10);
        assert!(!result.is_significant);
        assert_eq!(result.statistic, 0.0);
    }

    #[test]
    fn test_mcnemar_test_different_classifiers() {
        let y_true = array![[1, 0, 1], [0, 1, 0], [1, 1, 1], [0, 0, 0], [1, 0, 1]];
        let y_pred1 = array![[1, 0, 0], [0, 1, 1], [1, 0, 1], [1, 0, 0], [1, 1, 1]];
        let y_pred2 = array![[0, 1, 1], [1, 0, 0], [0, 1, 0], [0, 1, 1], [0, 0, 0]];

        let result = mcnemar_test(&y_true.view(), &y_pred1.view(), &y_pred2.view())
            .expect("operation should succeed");
        assert_eq!(result.test_name, "McNemar");
        assert!(result.statistic >= 0.0);
        assert!(result.p_value >= 0.0 && result.p_value <= 1.0);

        // Check additional info
        assert!(result.additional_info.contains_key("n01"));
        assert!(result.additional_info.contains_key("n10"));
        assert!(result.additional_info.contains_key("total_disagreements"));
    }

    #[test]
    fn test_mcnemar_test_error_handling() {
        let y_true = array![[1, 0], [0, 1]];
        let y_pred1 = array![[1, 0], [0, 1]];
        let y_pred2 = array![[1, 0, 1], [0, 1, 0]]; // Different shape

        let result = mcnemar_test(&y_true.view(), &y_pred1.view(), &y_pred2.view());
        assert!(result.is_err());
    }

    #[test]
    fn test_paired_t_test() {
        let metric_values1 = vec![0.8, 0.7, 0.9, 0.6, 0.75];
        let metric_values2 = vec![0.7, 0.65, 0.85, 0.55, 0.7];

        let result =
            paired_t_test(&metric_values1, &metric_values2).expect("operation should succeed");
        assert_eq!(result.test_name, "Paired t-test");
        assert!(result.p_value >= 0.0 && result.p_value <= 1.0);

        // Check additional info
        assert!(result.additional_info.contains_key("mean_difference"));
        assert!(result.additional_info.contains_key("std_dev_diff"));
        assert!(result.additional_info.contains_key("degrees_of_freedom"));
        assert!(result.additional_info.contains_key("n_observations"));

        // Mean difference should be positive (values1 > values2)
        let mean_diff = result
            .additional_info
            .get("mean_difference")
            .expect("index should be valid");
        assert!(*mean_diff > 0.0);
    }

    #[test]
    fn test_paired_t_test_identical_values() {
        let metric_values = vec![0.8, 0.7, 0.9, 0.6, 0.75];

        let result =
            paired_t_test(&metric_values, &metric_values).expect("operation should succeed");

        // With identical values, mean difference should be 0 and p-value should be 1.0
        let mean_diff = result
            .additional_info
            .get("mean_difference")
            .expect("index should be valid");
        assert!(mean_diff.abs() < 1e-10);
        assert!(!result.is_significant);
    }

    #[test]
    fn test_paired_t_test_error_handling() {
        let values1 = vec![0.8, 0.7];
        let values2 = vec![0.6]; // Different length

        let result = paired_t_test(&values1, &values2);
        assert!(result.is_err());

        // Too few observations
        let single = vec![0.8];
        let result = paired_t_test(&single, &single);
        assert!(result.is_err());
    }

    #[test]
    fn test_wilcoxon_signed_rank_test() {
        let metric_values1 = vec![0.8, 0.7, 0.9, 0.6, 0.75, 0.85, 0.65];
        let metric_values2 = vec![0.7, 0.65, 0.85, 0.55, 0.7, 0.8, 0.6];

        let result = wilcoxon_signed_rank_test(&metric_values1, &metric_values2)
            .expect("operation should succeed");
        assert_eq!(result.test_name, "Wilcoxon signed-rank");
        assert!(result.p_value >= 0.0 && result.p_value <= 1.0);

        // Check additional info
        assert!(result.additional_info.contains_key("w_plus"));
        assert!(result.additional_info.contains_key("w_minus"));
        assert!(result.additional_info.contains_key("n_nonzero_differences"));
        assert!(result.additional_info.contains_key("z_score"));
    }

    #[test]
    fn test_wilcoxon_signed_rank_test_identical_values() {
        let metric_values = vec![0.8, 0.7, 0.9, 0.6, 0.75];

        // Should fail with all zero differences
        let result = wilcoxon_signed_rank_test(&metric_values, &metric_values);
        assert!(result.is_err());
    }

    #[test]
    fn test_wilcoxon_signed_rank_test_error_handling() {
        let values1 = vec![0.8, 0.7];
        let values2 = vec![0.6]; // Different length

        let result = wilcoxon_signed_rank_test(&values1, &values2);
        assert!(result.is_err());

        // Too few observations
        let few = vec![0.8, 0.7];
        let result = wilcoxon_signed_rank_test(&few, &few);
        assert!(result.is_err());
    }

    #[test]
    fn test_confidence_interval() {
        let metric_values = vec![0.8, 0.75, 0.85, 0.7, 0.9, 0.65, 0.8, 0.82, 0.78, 0.88];

        let ci = confidence_interval(&metric_values, 0.95).expect("operation should succeed");

        assert_eq!(ci.confidence_level, 0.95);
        assert!(ci.lower <= ci.point_estimate); // Allow equal in edge cases
        assert!(ci.point_estimate <= ci.upper);

        // Point estimate should be the mean
        let expected_mean = metric_values.iter().sum::<f64>() / metric_values.len() as f64;
        assert!((ci.point_estimate - expected_mean).abs() < 1e-10);

        // For meaningful confidence intervals with sufficient data, bounds should be different
        if metric_values.len() > 5 {
            assert!(ci.upper - ci.lower > 1e-6); // Should have some width
        }

        // Test different confidence levels
        let ci_99 = confidence_interval(&metric_values, 0.99).expect("operation should succeed");
        assert!(ci_99.upper - ci_99.lower >= ci.upper - ci.lower); // Should be wider or equal
    }

    #[test]
    fn test_confidence_interval_single_value() {
        let single_value = vec![0.8];

        let ci = confidence_interval(&single_value, 0.95).expect("operation should succeed");

        // With single value, lower and upper should equal the point estimate
        assert_eq!(ci.lower, ci.point_estimate);
        assert_eq!(ci.upper, ci.point_estimate);
        assert_eq!(ci.point_estimate, 0.8);
    }

    #[test]
    fn test_confidence_interval_error_handling() {
        let empty = vec![];
        let result = confidence_interval(&empty, 0.95);
        assert!(result.is_err());

        let values = vec![0.8, 0.7];
        let result = confidence_interval(&values, 0.0); // Invalid confidence level
        assert!(result.is_err());

        let result = confidence_interval(&values, 1.0); // Invalid confidence level
        assert!(result.is_err());
    }

    #[test]
    fn test_statistical_test_result_creation() {
        let mut info = HashMap::new();
        info.insert("degrees_of_freedom".to_string(), 9.0);

        let result = StatisticalTestResult::new(2.5, 0.03, "Test".to_string(), Some(info));

        assert_eq!(result.statistic, 2.5);
        assert_eq!(result.p_value, 0.03);
        assert!(result.is_significant); // p < 0.05
        assert_eq!(result.test_name, "Test");
        assert_eq!(result.additional_info.get("degrees_of_freedom"), Some(&9.0));

        // Test non-significant result
        let non_sig = StatisticalTestResult::new(1.2, 0.15, "Non-sig".to_string(), None);
        assert!(!non_sig.is_significant); // p >= 0.05
        assert!(non_sig.additional_info.is_empty());
    }

    #[test]
    fn test_distribution_helper_functions() {
        // Test standard normal CDF
        let z_zero = standard_normal_cdf(0.0);
        assert!((z_zero - 0.5).abs() < 1e-6);

        let z_positive = standard_normal_cdf(1.96);
        assert!((z_positive - 0.975).abs() < 0.01); // Approximately 0.975

        // Test normal quantile - debug the actual values first
        let q_median = normal_quantile(0.5);
        assert!(q_median.abs() < 1e-6); // Should be approximately 0

        let q_975 = normal_quantile(0.975);
        assert!((q_975 - 1.96).abs() < 0.01); // Should be approximately 1.96

        let q_025 = normal_quantile(0.025);
        assert!((q_025 + 1.96).abs() < 0.01); // Should be approximately -1.96

        // Test edge cases
        assert_eq!(normal_quantile(0.0), f64::NEG_INFINITY);
        assert_eq!(normal_quantile(1.0), f64::INFINITY);

        // Test that quantile and CDF are approximately inverse functions
        let test_values = vec![0.25, 0.5, 0.75]; // Use fewer test values
        for &p in &test_values {
            let q = normal_quantile(p);
            if q.is_finite() {
                let p_back = standard_normal_cdf(q);
                assert!((p - p_back).abs() < 0.05); // More lenient round-trip accuracy
            }
        }
    }

    #[test]
    fn test_existing_metrics_compatibility() {
        let (y_true, y_pred) = create_test_data();
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        // Test that existing metrics still work
        let hamming = hamming_loss(&y_true_view, &y_pred_view).expect("operation should succeed");
        assert!((0.0..=1.0).contains(&hamming));

        let subset_acc =
            subset_accuracy(&y_true_view, &y_pred_view).expect("operation should succeed");
        assert!((0.0..=1.0).contains(&subset_acc));

        let jaccard = jaccard_score(&y_true_view, &y_pred_view).expect("operation should succeed");
        assert!((0.0..=1.0).contains(&jaccard));

        let f1_micro =
            f1_score(&y_true_view, &y_pred_view, "micro").expect("operation should succeed");
        assert!((0.0..=1.0).contains(&f1_micro));

        let f1_macro =
            f1_score(&y_true_view, &y_pred_view, "macro").expect("operation should succeed");
        assert!((0.0..=1.0).contains(&f1_macro));

        let f1_samples =
            f1_score(&y_true_view, &y_pred_view, "samples").expect("sampling should succeed");
        assert!((0.0..=1.0).contains(&f1_samples));
    }

    #[test]
    fn test_per_label_vs_global_metrics_consistency() {
        let (y_true, y_pred) = create_test_data();
        let y_true_view = y_true.view();
        let y_pred_view = y_pred.view();

        let per_label =
            per_label_metrics(&y_true_view, &y_pred_view).expect("operation should succeed");
        let global_f1_macro =
            f1_score(&y_true_view, &y_pred_view, "macro").expect("operation should succeed");
        let per_label_f1_macro = per_label
            .macro_average("f1_score")
            .expect("operation should succeed");

        // The macro F1 from per-label metrics should match global macro F1
        assert!((global_f1_macro - per_label_f1_macro).abs() < 1e-10);
    }

    #[test]
    fn test_comprehensive_statistical_workflow() {
        // Simulate a complete statistical comparison workflow
        let y_true = array![
            [1, 0, 1, 0],
            [0, 1, 0, 1],
            [1, 1, 1, 0],
            [0, 0, 0, 1],
            [1, 0, 1, 1]
        ];

        // Two different classifiers
        let y_pred1 = array![
            [1, 0, 0, 0],
            [0, 1, 1, 1],
            [1, 0, 1, 0],
            [1, 0, 0, 1],
            [1, 1, 1, 0]
        ];

        let y_pred2 = array![
            [0, 1, 1, 0],
            [1, 0, 0, 0],
            [0, 1, 0, 1],
            [0, 1, 1, 0],
            [0, 0, 0, 1]
        ];

        // Get per-label metrics for both classifiers
        let metrics1 =
            per_label_metrics(&y_true.view(), &y_pred1.view()).expect("operation should succeed");
        let metrics2 =
            per_label_metrics(&y_true.view(), &y_pred2.view()).expect("operation should succeed");

        // Compare using McNemar's test
        let mcnemar_result = mcnemar_test(&y_true.view(), &y_pred1.view(), &y_pred2.view())
            .expect("operation should succeed");

        // Compare F1 scores using t-test
        let t_test_result = paired_t_test(&metrics1.f1_score, &metrics2.f1_score)
            .expect("operation should succeed");

        // Compare F1 scores using Wilcoxon test
        let wilcoxon_result = wilcoxon_signed_rank_test(&metrics1.f1_score, &metrics2.f1_score)
            .expect("operation should succeed");

        // Get confidence interval for first classifier's F1 scores
        let ci_result =
            confidence_interval(&metrics1.f1_score, 0.95).expect("operation should succeed");

        // All results should be valid
        assert!(mcnemar_result.p_value >= 0.0 && mcnemar_result.p_value <= 1.0);
        assert!(t_test_result.p_value >= 0.0 && t_test_result.p_value <= 1.0);
        assert!(wilcoxon_result.p_value >= 0.0 && wilcoxon_result.p_value <= 1.0);
        assert!(ci_result.lower <= ci_result.point_estimate);
        assert!(ci_result.point_estimate <= ci_result.upper);

        // Results should contain proper metadata
        assert!(!mcnemar_result.additional_info.is_empty());
        assert!(!t_test_result.additional_info.is_empty());
        assert!(!wilcoxon_result.additional_info.is_empty());
    }
}