plotters-statistical 0.1.0

//! Precision–recall curve: threshold sweep producing (recall, precision) points
//! and average precision (AP).

use super::StatsError;

/// A single operating point on a precision–recall curve.
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct PrPoint {
    /// Recall (true positive rate), `TP / P`.
    pub recall: f64,
    /// Precision, `TP / (TP + FP)`.
    pub precision: f64,
    /// The score threshold at which this point is reached.
    pub threshold: f64,
}

/// A precision–recall curve: operating points, average precision, and the
/// positive-class prevalence (the curve's chance baseline).
#[derive(Debug, Clone, PartialEq)]
pub struct PrCurve {
    /// Operating points, ascending in recall.
    pub points: Vec<PrPoint>,
    /// Average precision — the step-function area under the PR curve, using the
    /// `sum_n (R_n - R_{n-1}) * P_n` definition (matching scikit-learn's
    /// `average_precision_score`), **not** a trapezoidal interpolation.
    pub average_precision: f64,
    /// Positive-class prevalence, `P / (P + N)`. A no-skill classifier's PR
    /// curve is a horizontal line at this height — the correct PR baseline (a
    /// diagonal, which is the ROC baseline, would be wrong here).
    pub baseline: f64,
}

/// Compute a precision–recall curve from predicted `scores` and true binary
/// `labels` (`true` = positive class). Higher scores are "more positive".
///
/// # Errors
///
/// * [`StatsError::LengthMismatch`] if the inputs differ in length.
/// * [`StatsError::EmptyInput`] if empty.
/// * [`StatsError::NoPositiveLabels`] if there are no positives, since recall is
///   `0/0`. (Unlike ROC, PR does not require any negatives.)
pub fn precision_recall_curve(scores: &[f64], labels: &[bool]) -> Result<PrCurve, StatsError> {
    if scores.len() != labels.len() {
        return Err(StatsError::LengthMismatch {
            scores: scores.len(),
            labels: labels.len(),
        });
    }
    if scores.is_empty() {
        return Err(StatsError::EmptyInput);
    }
    let total_p = labels.iter().filter(|&&l| l).count();
    if total_p == 0 {
        return Err(StatsError::NoPositiveLabels);
    }
    let baseline = total_p as f64 / labels.len() as f64;

    let mut order: Vec<usize> = (0..scores.len()).collect();
    order.sort_by(|&a, &b| {
        scores[b]
            .partial_cmp(&scores[a])
            .unwrap_or(std::cmp::Ordering::Equal)
    });

    let p = total_p as f64;
    let mut points = Vec::new();
    let mut tp = 0usize;
    let mut fp = 0usize;
    let mut i = 0usize;
    let mut prev_recall = 0.0;
    let mut average_precision = 0.0;

    while i < order.len() {
        let score = scores[order[i]];
        while i < order.len() && scores[order[i]] == score {
            if labels[order[i]] {
                tp += 1;
            } else {
                fp += 1;
            }
            i += 1;
        }
        let recall = tp as f64 / p;
        let precision = if tp + fp == 0 {
            1.0
        } else {
            tp as f64 / (tp + fp) as f64
        };
        average_precision += (recall - prev_recall) * precision;
        prev_recall = recall;
        points.push(PrPoint {
            recall,
            precision,
            threshold: score,
        });
    }

    Ok(PrCurve {
        points,
        average_precision,
        baseline,
    })
}