ferrolearn_linear/

ransac.rs

//! RANSAC (RANdom SAmple Consensus) robust regression.
//!
//! This module provides [`RANSACRegressor`], a meta-estimator that fits a
//! base regressor to inlier data, automatically detecting and excluding
//! outliers.
//!
//! # Algorithm
//!
//! Mirrors scikit-learn 1.5.2's `RANSACRegressor.fit` decision rule
//! (`sklearn/linear_model/_ransac.py:451-606`). Initialize `n_inliers_best = 1`,
//! `score_best = -inf`, `inlier_mask_best = None`. Then for each trial:
//!
//! 1. Draw a `min_samples`-sized subset of indices (seedable; RNG-sequence
//!    parity with numpy is out of scope — see `## REQ status`).
//! 2. Fit the base estimator on the SUBSET.
//! 3. Predict on ALL of `X`; classify a sample as an inlier iff
//!    `|y − y_pred| <= residual_threshold` (boundary inclusive).
//! 4. If `n_inliers_subset < n_inliers_best`, skip.
//! 5. Compute `score_subset` = R² of the SUBSET model on its inlier set
//!    (`1 − SS_res/SS_tot`; if `SS_tot == 0` then `1.0` when `SS_res == 0` else
//!    `0.0`). No refit inside the loop.
//! 6. If `n_inliers_subset == n_inliers_best` and `score_subset < score_best`,
//!    skip — so higher R² wins ties on inlier count.
//! 7. Otherwise record the new best (`n_inliers_best`, `score_best`, and the
//!    SUBSET model's mask — never recomputed from a refit).
//!
//! After the loop, refit the base estimator ONCE on the best inlier set; that is
//! the stored model. The reported `inlier_mask` is the winning subset model's
//! mask. The default `residual_threshold` is the MAD of `y`
//! (`median(|y − median(y)|)`), which may be exactly `0` for a constant target.
//!
//! ## REQ status (per `.design/linear/ransac.md`, mirrors `sklearn/linear_model/_ransac.py` @ 1.5.2)
//!
//! Binary classification (R-DEFER-2): SHIPPED means impl plus non-test consumer
//! plus tests, all green; NOT-STARTED means an open blocker referenced by number.
//! The boundary estimator types [`RANSACRegressor`] and
//! [`FittedRANSACRegressor`] are re-exported at the crate root
//! (`pub use ransac::{...} in lib.rs`); under S5/R-DEFER-1 the public estimator
//! type IS the consumer surface (no `ferrolearn-python` RANSAC binding yet).
//!
//! **RNG non-parity caveat:** subset draws use `rand::rngs::StdRng` (Fisher-Yates),
//! NOT numpy's Mersenne-Twister `sample_without_replacement`. Same-seed
//! cross-implementation subset-sequence parity is infeasible and out of scope;
//! parity is asserted only on the deterministic decision rules below.
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (sampling loop) | SHIPPED | `fn sample_indices` draws `k` distinct indices via Fisher-Yates, called per trial in `fn fit`, seeded deterministically. Test: `test_ransac_reproducible_with_seed`. Structural only (RNG caveat above). |
//! | REQ-2 (MAD threshold default) | SHIPPED | `fn fit` sets the auto threshold to `fn mad` (`median(|y − median(y)|)`) when `residual_threshold` is `None`, mirroring `_ransac.py:401`. Test: `test_ransac_auto_threshold`. |
//! | REQ-3 (inlier classification) | SHIPPED | `fn fit`: `if (preds[i] - y[i]).abs() <= threshold { inlier_mask_subset[i] = true }`, boundary-inclusive `<=` per `_ransac.py:511`. Tests: `test_ransac_with_outlier`, `test_ransac_multiple_outliers`. |
//! | REQ-4 (selection: n_inliers then R²) | SHIPPED | `fn fit` ranks by `(n_inliers_subset, score_subset)` with `score_subset = fn r2_score` of the subset model on its inliers; ties skip when `score_subset < score_best` (higher R² wins), mirroring `_ransac.py:530-543`. Test: `ransac_selection_criterion_r2_not_residual_sum` (tests/divergence_ransac_fit.rs) — oracle picks group B `[F,F,F,T,T,T]`, predict([[1.0]])≈10.05. Closed #512. |
//! | REQ-5 (refit-once; mask from subset model) | SHIPPED | `fn fit` records `inlier_mask_best` from the SUBSET model (no in-loop refit/recompute) and refits the base estimator ONCE after the loop on `(x_inlier_best, y_inlier_best)`, mirroring `_ransac.py:544,602,605`. The stored `inlier_mask` is never recomputed from the refit. Verified by the green divergence suite + module unit tests. Closed #513. |
//! | REQ-6 (n_inliers_best init / acceptance gate) | SHIPPED | `fn fit` initializes `n_inliers_best = 1` (`_ransac.py:451`), skips only when `n_inliers_subset < n_inliers_best` (`_ransac.py:515`), and no longer gates on `n_inliers >= min_samples`. The up-front `n_samples < min_samples` guard mirrors `_ransac.py:393-397`. Closed #514. |
//! | REQ-7 (dynamic max_trials + stop criteria) | NOT-STARTED | open blocker #515. `fn fit` runs a FIXED `for _ in 0..self.max_trials` loop; no `_dynamic_max_trials` shrink, no `stop_n_inliers`/`stop_score`/`stop_probability`/`max_skips`, no `n_trials_`/`n_skips_*` tracking. |
//! | REQ-8 (loss='squared_error') | SHIPPED | impl: `pub enum RansacLoss { #[default] AbsoluteError, SquaredError }` + field `loss: RansacLoss` + `with_loss` builder (sklearn default `'absolute_error'`, `_ransac.py:301`). `fn fit` branches the per-sample residual: `AbsoluteError → (preds[i]-y[i]).abs()`, `SquaredError → { let d = preds[i]-y[i]; d*d }`, mirroring `_ransac.py:407,414` and applied at the `residuals <= residual_threshold` classification (`_ransac.py:508,511`); the MAD-default threshold stays loss-independent (`_ransac.py:399-401`). Consumer: boundary types re-exported at crate root + `RansacLoss` re-exported (`pub use ransac::{...RansacLoss} in lib.rs`). Tests (live-oracle, RNG-independent): `ransac_loss_squared_error_recovers_line`, `ransac_loss_default_absolute_error_byte_identical` (tests/divergence_ransac_fit.rs). Closed #516. |
//! | REQ-9 (MAD-zero parity) | SHIPPED | `fn fit` uses the MAD value directly (`mad(&y.to_vec())`) with no `1e-6` substitution, so a constant target yields threshold `0.0` per `_ransac.py:399-401`. Test: `ransac_mad_zero_threshold_excludes_tiny_deviation` (tests/divergence_ransac_fit.rs) — idx 7 (residual 1e-7) is an OUTLIER. Closed #517. |
//! | REQ-10 (introspection attributes) | NOT-STARTED | open blocker #518. `FittedRANSACRegressor` exposes only `inlier_mask()`; no `estimator_`/`n_trials_`/`n_skips_*`/`n_features_in_`. |
//! | REQ-11 (is_data_valid / is_model_valid / max_skips) | NOT-STARTED | open blocker #519. `RANSACRegressor` has no such fields. |
//! | REQ-12 (min_samples float fraction) | SHIPPED | impl: `pub enum MinSamples<F> { Count(usize), Fraction(F) }` + field `min_samples: Option<MinSamples<F>>` + builders `with_min_samples(usize) → Count`, `with_min_samples_fraction(F) → Fraction`, getter `min_samples()`. `fn fit` resolves `None → n_features+1` (unchanged), `Count(k) → k`, `Fraction(f) → ceil(f·n_samples)` validating `0 < f < 1` (else `FerroError::InvalidParameter`), with the resolved-count `> n_samples → InvalidParameter` guard mirroring `_ransac.py:382-397` (sklearn `ValueError`). Consumer: boundary types re-exported at crate root + `MinSamples` re-exported (`pub use ransac::{...MinSamples} in lib.rs`). Tests (live-oracle): `ransac_min_samples_fraction_resolves_ceil`, `ransac_min_samples_fraction_out_of_range_errors`, `ransac_min_samples_count_unchanged` (tests/divergence_ransac_fit.rs). Closed #520. |
//! | REQ-13 (ferray substrate) | NOT-STARTED | open blocker #521. Still on `ndarray` + `rand::rngs::StdRng`, not `ferray-core`/`ferray::random`. |
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::ransac::RANSACRegressor;
//! use ferrolearn_linear::LinearRegression;
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! // Data with an outlier at index 4.
//! let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
//! let y = array![2.0, 4.0, 6.0, 8.0, 100.0]; // last point is outlier
//!
//! let base = LinearRegression::<f64>::new();
//! let model = RANSACRegressor::new(base);
//! let fitted = model.fit(&x, &y).unwrap();
//!
//! // The outlier should be detected.
//! let mask = fitted.inlier_mask();
//! assert!(!mask[4], "outlier at index 4 should be detected");
//! ```

use ferrolearn_core::error::FerroError;
use ferrolearn_core::traits::{Fit, Predict};
use ndarray::{Array1, Array2, ScalarOperand};
use num_traits::Float;
use rand::Rng;
use rand::SeedableRng;

// ---------------------------------------------------------------------------
// Loss family (REQ-8)
// ---------------------------------------------------------------------------

/// Per-sample residual function used to classify inliers.
///
/// Mirrors scikit-learn's `RANSACRegressor(loss=...)` string options
/// (`sklearn/linear_model/_ransac.py:284,405-418`), constrained by
/// `StrOptions({"absolute_error", "squared_error"})`. The residual produced
/// here is compared (boundary-inclusive) against `residual_threshold`; the
/// MAD-based default threshold itself is computed identically regardless of the
/// loss (`_ransac.py:399-401` is unconditional — only the per-sample residual at
/// `_ransac.py:508` changes).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum RansacLoss {
    /// Absolute error per sample: `|y_true − y_pred|`
    /// (`_ransac.py:407`, sklearn default).
    #[default]
    AbsoluteError,
    /// Squared error per sample: `(y_true − y_pred)²` (`_ransac.py:414`).
    SquaredError,
}

// ---------------------------------------------------------------------------
// min_samples specification (REQ-12)
// ---------------------------------------------------------------------------

/// How `min_samples` (the per-trial subset size) is specified.
///
/// Mirrors scikit-learn's `min_samples` parameter, which is `int (>= 1)` for an
/// absolute count or `float ([0, 1])` for a relative fraction
/// (`sklearn/linear_model/_ransac.py:115-125`; constraint
/// `Interval(Integral, 1, None) | Interval(RealNotInt, 0, 1, closed="both")`,
/// `_ransac.py:262-266`). The resolution at fit time is:
/// `0 < f < 1 → ceil(f · n_samples)` (`_ransac.py:389-390`); an integer count is
/// used directly (`_ransac.py:391-392`). A resolved count larger than
/// `n_samples` is a `ValueError` (`_ransac.py:393-397`).
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum MinSamples<F> {
    /// An absolute number of samples per subset (`min_samples >= 1`).
    Count(usize),
    /// A fraction of `n_samples` per subset, resolved to
    /// `ceil(fraction · n_samples)` (sklearn `0 < min_samples < 1`).
    Fraction(F),
}

// ---------------------------------------------------------------------------
// RANSACRegressor (unfitted)
// ---------------------------------------------------------------------------

/// RANSAC robust regression meta-estimator.
///
/// Wraps a base regressor (e.g., [`LinearRegression`](crate::LinearRegression))
/// and repeatedly fits it on random subsets to find a model robust to
/// outliers.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
/// - `E`: The base estimator type.
#[derive(Debug, Clone)]
pub struct RANSACRegressor<F, E> {
    /// The base estimator.
    pub estimator: E,
    /// Minimum number of samples per subset. `None` resolves to
    /// `n_features + 1` (the `LinearRegression` default,
    /// `sklearn/linear_model/_ransac.py:388`); a [`MinSamples::Count`] is an
    /// absolute count, a [`MinSamples::Fraction`] resolves to
    /// `ceil(fraction · n_samples)` (`_ransac.py:389-390`).
    pub min_samples: Option<MinSamples<F>>,
    /// Residual threshold: points whose per-sample residual (under [`loss`]) is
    /// `<= threshold` are considered inliers. If `None`, uses the MAD of the
    /// target.
    ///
    /// [`loss`]: RANSACRegressor::loss
    pub residual_threshold: Option<F>,
    /// Maximum number of random trials.
    pub max_trials: usize,
    /// Per-sample residual loss used for inlier classification. Defaults to
    /// [`RansacLoss::AbsoluteError`] (sklearn default,
    /// `sklearn/linear_model/_ransac.py:301`).
    pub loss: RansacLoss,
    /// Optional random seed for reproducibility.
    pub random_state: Option<u64>,
}

impl<F: Float, E> RANSACRegressor<F, E> {
    /// Create a new `RANSACRegressor` with the given base estimator.
    ///
    /// Defaults: `min_samples = None` (auto: n_features + 1),
    /// `residual_threshold = None` (auto: MAD), `max_trials = 100`,
    /// `loss = AbsoluteError` (sklearn default), `random_state = None`.
    #[must_use]
    pub fn new(estimator: E) -> Self {
        Self {
            estimator,
            min_samples: None,
            residual_threshold: None,
            max_trials: 100,
            loss: RansacLoss::AbsoluteError,
            random_state: None,
        }
    }

    /// Set the minimum number of samples per subset as an absolute count
    /// ([`MinSamples::Count`]). Mirrors sklearn `min_samples` as an `int >= 1`.
    #[must_use]
    pub fn with_min_samples(mut self, min_samples: usize) -> Self {
        self.min_samples = Some(MinSamples::Count(min_samples));
        self
    }

    /// Set the minimum number of samples per subset as a fraction of
    /// `n_samples` ([`MinSamples::Fraction`]), resolved at fit time to
    /// `ceil(fraction · n_samples)`. Mirrors sklearn `min_samples` as a
    /// `float` in `(0, 1)` (`sklearn/linear_model/_ransac.py:389-390`).
    ///
    /// The fraction is validated at fit time: a value outside `(0, 1)` (or one
    /// whose resolved count exceeds `n_samples`) yields
    /// [`FerroError::InvalidParameter`], mirroring sklearn's `ValueError`
    /// (`_ransac.py:393-397`).
    #[must_use]
    pub fn with_min_samples_fraction(mut self, fraction: F) -> Self {
        self.min_samples = Some(MinSamples::Fraction(fraction));
        self
    }

    /// Returns the configured `min_samples` specification (`None` = auto).
    #[must_use]
    pub fn min_samples(&self) -> Option<MinSamples<F>> {
        self.min_samples
    }

    /// Set the per-sample residual loss used for inlier classification.
    ///
    /// Mirrors sklearn `RANSACRegressor(loss=...)`
    /// (`sklearn/linear_model/_ransac.py:284`). Defaults to
    /// [`RansacLoss::AbsoluteError`].
    #[must_use]
    pub fn with_loss(mut self, loss: RansacLoss) -> Self {
        self.loss = loss;
        self
    }

    /// Set the residual threshold for inlier detection.
    #[must_use]
    pub fn with_residual_threshold(mut self, threshold: F) -> Self {
        self.residual_threshold = Some(threshold);
        self
    }

    /// Set the maximum number of random trials.
    #[must_use]
    pub fn with_max_trials(mut self, max_trials: usize) -> Self {
        self.max_trials = max_trials;
        self
    }

    /// Set the random seed for reproducibility.
    #[must_use]
    pub fn with_random_state(mut self, seed: u64) -> Self {
        self.random_state = Some(seed);
        self
    }
}

// ---------------------------------------------------------------------------
// FittedRANSACRegressor
// ---------------------------------------------------------------------------

/// Fitted RANSAC robust regression model.
///
/// Stores the best estimator fitted on inlier data, and the inlier mask.
#[derive(Debug, Clone)]
pub struct FittedRANSACRegressor<Fitted> {
    /// The fitted base estimator (fitted on inliers).
    fitted_estimator: Fitted,
    /// Boolean mask: true if the sample was classified as an inlier.
    inlier_mask: Vec<bool>,
}

impl<Fitted> FittedRANSACRegressor<Fitted> {
    /// Returns the inlier mask. `true` indicates the sample was an inlier.
    #[must_use]
    pub fn inlier_mask(&self) -> &[bool] {
        &self.inlier_mask
    }
}

// ---------------------------------------------------------------------------
// Helper: Median Absolute Deviation
// ---------------------------------------------------------------------------

/// Compute the median of a slice of floats.
fn median<F: Float>(values: &[F]) -> F {
    let mut sorted: Vec<F> = values.to_vec();
    // Total order without `.unwrap()`: NaNs (absent for valid targets) sort last.
    sorted.sort_by(|a, b| a.partial_cmp(b).unwrap_or(core::cmp::Ordering::Equal));
    let n = sorted.len();
    if n == 0 {
        return F::zero();
    }
    if n.is_multiple_of(2) {
        (sorted[n / 2 - 1] + sorted[n / 2]) / (F::one() + F::one())
    } else {
        sorted[n / 2]
    }
}

/// Compute the Median Absolute Deviation (MAD) of a slice.
fn mad<F: Float>(values: &[F]) -> F {
    let med = median(values);
    let abs_devs: Vec<F> = values.iter().map(|&v| (v - med).abs()).collect();
    median(&abs_devs)
}

/// Coefficient of determination R² of `y_pred` against `y_true`.
///
/// Mirrors sklearn's `r2_score` (used through `estimator.score`,
/// `_ransac.py:530`): `R² = 1 - SS_res / SS_tot` where
/// `SS_res = Σ(y_true − y_pred)²` and `SS_tot = Σ(y_true − mean(y_true))²`.
///
/// Matches sklearn's constant-target edge case (`metrics/_regression.py`):
/// when `SS_tot == 0`, R² is `1.0` if `SS_res == 0` (perfect prediction) and
/// `0.0` otherwise.
fn r2_score<F: Float>(y_true: &[F], y_pred: &[F]) -> F {
    let n = y_true.len();
    if n == 0 {
        return F::zero();
    }
    let mut sum = F::zero();
    for &v in y_true {
        sum = sum + v;
    }
    let mean = sum / F::from(n).unwrap_or_else(F::one);

    let mut ss_res = F::zero();
    let mut ss_tot = F::zero();
    for (&t, &p) in y_true.iter().zip(y_pred.iter()) {
        let res = t - p;
        ss_res = ss_res + res * res;
        let dev = t - mean;
        ss_tot = ss_tot + dev * dev;
    }

    if ss_tot == F::zero() {
        if ss_res == F::zero() {
            F::one()
        } else {
            F::zero()
        }
    } else {
        F::one() - ss_res / ss_tot
    }
}

// ---------------------------------------------------------------------------
// Random subset sampling
// ---------------------------------------------------------------------------

/// Sample `k` distinct indices from `0..n` using Fisher-Yates.
fn sample_indices<R: Rng>(rng: &mut R, n: usize, k: usize) -> Vec<usize> {
    let mut indices: Vec<usize> = (0..n).collect();
    for i in 0..k {
        let j = rng.random_range(i..n);
        indices.swap(i, j);
    }
    indices.truncate(k);
    indices
}

/// Extract a subset of rows from a 2D array and a 1D array.
fn subset<F: Float>(x: &Array2<F>, y: &Array1<F>, indices: &[usize]) -> (Array2<F>, Array1<F>) {
    let n_features = x.ncols();
    let n = indices.len();
    let mut x_sub = Array2::<F>::zeros((n, n_features));
    let mut y_sub = Array1::<F>::zeros(n);
    for (row, &idx) in indices.iter().enumerate() {
        for col in 0..n_features {
            x_sub[[row, col]] = x[[idx, col]];
        }
        y_sub[row] = y[idx];
    }
    (x_sub, y_sub)
}

// ---------------------------------------------------------------------------
// Fit and Predict
// ---------------------------------------------------------------------------

impl<F, E, Ef> Fit<Array2<F>, Array1<F>> for RANSACRegressor<F, E>
where
    F: Float + Send + Sync + ScalarOperand + num_traits::FromPrimitive + 'static,
    E: Fit<Array2<F>, Array1<F>, Fitted = Ef, Error = FerroError> + Clone,
    Ef: Predict<Array2<F>, Output = Array1<F>, Error = FerroError> + Clone,
{
    type Fitted = FittedRANSACRegressor<Ef>;
    type Error = FerroError;

    /// Fit the RANSAC model by repeatedly sampling and fitting.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have different
    /// sample counts.
    /// Returns [`FerroError::ConvergenceFailure`] if no valid model is found
    /// after `max_trials` iterations.
    fn fit(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
    ) -> Result<FittedRANSACRegressor<E::Fitted>, FerroError> {
        let (n_samples, n_features) = x.dim();

        if n_samples != y.len() {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples],
                actual: vec![y.len()],
                context: "y length must match number of samples in X".into(),
            });
        }

        // Resolve `min_samples` per sklearn (`_ransac.py:382-397`):
        //   None        -> n_features + 1 (the LinearRegression default).
        //   Count(k)     -> k (integer count, `>= 1` branch); `k < 1` is a
        //                   ValueError per the parameter constraint
        //                   `Interval(Integral, 1, None, closed="left")`
        //                   (`_ransac.py:263`) — never silently coerced.
        //   Fraction(f)  -> ceil(f * n_samples) for `0 < f < 1`; a fraction
        //                   outside (0, 1) is a ValueError.
        // A resolved count `> n_samples` is a ValueError (`_ransac.py:393-397`).
        let min_samples = match self.min_samples {
            None => (n_features + 1).max(1),
            Some(MinSamples::Count(k)) => {
                // sklearn rejects `min_samples < 1` at parameter validation
                // (`Interval(Integral, 1, None, closed="left")`,
                // `_ransac.py:263`); do not coerce `0` to `1`.
                if k < 1 {
                    return Err(FerroError::InvalidParameter {
                        name: "min_samples".into(),
                        reason: "min_samples must be >= 1".into(),
                    });
                }
                k
            }
            Some(MinSamples::Fraction(f)) => {
                // sklearn's float branch is `0 < min_samples < 1`
                // (`_ransac.py:389`); a float `>= 1` would take the integer
                // branch — for the explicit `Fraction` variant we require a
                // genuine fraction in (0, 1).
                if !(f > F::zero() && f < F::one()) {
                    return Err(FerroError::InvalidParameter {
                        name: "min_samples".into(),
                        reason: "min_samples fraction must be in the open \
                                 interval (0, 1); use with_min_samples for an \
                                 absolute count"
                            .into(),
                    });
                }
                // ceil(f * n_samples) (`_ransac.py:390`).
                let raw = f * F::from(n_samples).unwrap_or_else(F::one);
                let resolved = raw.ceil();
                // n_samples >= 1 and 0 < f < 1 keep `resolved` finite and >= 1.
                let resolved = resolved.to_usize().unwrap_or(1).max(1);
                if resolved > n_samples {
                    return Err(FerroError::InvalidParameter {
                        name: "min_samples".into(),
                        reason: format!(
                            "`min_samples` may not be larger than number of \
                             samples: n_samples = {n_samples}."
                        ),
                    });
                }
                resolved
            }
        };

        if n_samples < min_samples {
            return Err(FerroError::InsufficientSamples {
                required: min_samples,
                actual: n_samples,
                context: "RANSAC requires at least min_samples samples".into(),
            });
        }

        // Compute residual threshold if not provided.
        //
        // sklearn (`_ransac.py:399-401`) sets the default threshold to the MAD
        // (median absolute deviation) of `y` with NO special-casing of zero:
        // `residual_threshold = np.median(np.abs(y - np.median(y)))`. A constant
        // (or near-constant) target therefore yields a threshold of exactly 0.0,
        // under which only samples with a zero residual are inliers.
        let threshold = match self.residual_threshold {
            Some(t) => t,
            None => mad(&y.to_vec()),
        };

        let mut rng = match self.random_state {
            Some(seed) => rand::rngs::StdRng::seed_from_u64(seed),
            None => rand::rngs::StdRng::seed_from_u64(42),
        };

        // sklearn-faithful selection state (`_ransac.py:451-456`):
        //   n_inliers_best = 1, score_best = -inf, inlier_mask_best = None.
        // The best inlier index set is remembered for the single final refit.
        let mut n_inliers_best: usize = 1;
        let mut score_best = F::neg_infinity();
        let mut inlier_mask_best: Option<Vec<bool>> = None;
        let mut inlier_best_idxs: Option<Vec<usize>> = None;

        // `while n_trials < max_trials` (`_ransac.py:467`). We keep the fixed
        // `self.max_trials` loop; dynamic max-trials / stop criteria are #515.
        for _ in 0..self.max_trials {
            // Choose a random sample set (`_ransac.py:478-482`). RNG-sequence
            // parity with numpy's `sample_without_replacement` is infeasible and
            // explicitly out of scope (see module REQ status).
            let subset_idxs = sample_indices(&mut rng, n_samples, min_samples);
            let (x_subset, y_subset) = subset(x, y, &subset_idxs);

            // Fit the base estimator on the SUBSET (`_ransac.py:497`).
            let fitted_subset = match self.estimator.fit(&x_subset, &y_subset) {
                Ok(f) => f,
                Err(_) => continue, // Skip failed fits (degenerate subset).
            };

            // Residuals of ALL data under the subset model (`_ransac.py:507-508`).
            let preds = match fitted_subset.predict(x) {
                Ok(p) => p,
                Err(_) => continue,
            };

            // Per-sample residual under the configured loss (`_ransac.py:508`,
            // `residuals_subset = loss_function(y, y_pred)`):
            //   AbsoluteError -> |y − y_pred|  (`_ransac.py:407`)
            //   SquaredError  -> (y − y_pred)² (`_ransac.py:414`)
            // The residual is then compared, boundary-inclusive, to
            // `residual_threshold` (`_ransac.py:511-512`); the MAD-default
            // threshold is loss-independent (`_ransac.py:399-401`).
            let mut inlier_mask_subset = vec![false; n_samples];
            let mut inlier_idxs_subset: Vec<usize> = Vec::new();
            for i in 0..n_samples {
                let residual = match self.loss {
                    RansacLoss::AbsoluteError => (preds[i] - y[i]).abs(),
                    RansacLoss::SquaredError => {
                        let d = preds[i] - y[i];
                        d * d
                    }
                };
                if residual <= threshold {
                    inlier_mask_subset[i] = true;
                    inlier_idxs_subset.push(i);
                }
            }
            let n_inliers_subset = inlier_idxs_subset.len();

            // Fewer inliers than the best so far -> skip (`_ransac.py:514-517`).
            if n_inliers_subset < n_inliers_best {
                continue;
            }

            // Score the SUBSET model on the inlier set: R² of the subset-fitted
            // model on `(X_inlier_subset, y_inlier_subset)` (`_ransac.py:530-534`).
            // No refit inside the loop.
            let y_inlier_subset: Vec<F> = inlier_idxs_subset.iter().map(|&i| y[i]).collect();
            let pred_inlier_subset: Vec<F> = inlier_idxs_subset.iter().map(|&i| preds[i]).collect();
            let score_subset = r2_score(&y_inlier_subset, &pred_inlier_subset);

            // Same inlier count but worse score -> skip (`_ransac.py:538-539`).
            // Higher R² wins ties on inlier count.
            if n_inliers_subset == n_inliers_best && score_subset < score_best {
                continue;
            }

            // Record the new best (`_ransac.py:542-547`). The stored mask is the
            // SUBSET model's mask — NOT recomputed from a refit.
            n_inliers_best = n_inliers_subset;
            score_best = score_subset;
            inlier_mask_best = Some(inlier_mask_subset);
            inlier_best_idxs = Some(inlier_idxs_subset);
        }

        // No valid consensus set found (`_ransac.py:561-580`).
        let (mask_best, idxs_best) = match (inlier_mask_best, inlier_best_idxs) {
            (Some(mask), Some(idxs)) => (mask, idxs),
            _ => {
                return Err(FerroError::ConvergenceFailure {
                    iterations: self.max_trials,
                    message: "RANSAC could not find a valid model after max_trials iterations"
                        .into(),
                });
            }
        };

        // Estimate the final model using the best inlier set ONCE, after the
        // loop (`_ransac.py:597-602`). `inlier_mask_` stays the subset model's
        // mask; it is never recomputed from this refit.
        let (x_inlier_best, y_inlier_best) = subset(x, y, &idxs_best);
        let fitted_estimator = self.estimator.fit(&x_inlier_best, &y_inlier_best)?;

        Ok(FittedRANSACRegressor {
            fitted_estimator,
            inlier_mask: mask_best,
        })
    }
}

impl<F, Fitted> Predict<Array2<F>> for FittedRANSACRegressor<Fitted>
where
    F: Float + Send + Sync + 'static,
    Fitted: Predict<Array2<F>, Output = Array1<F>, Error = FerroError>,
{
    type Output = Array1<F>;
    type Error = FerroError;

    /// Predict target values using the base estimator fitted on inliers.
    ///
    /// # Errors
    ///
    /// Returns any error from the base estimator's predict method.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        self.fitted_estimator.predict(x)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::LinearRegression;
    use approx::assert_relative_eq;
    use ndarray::array;

    #[test]
    fn test_ransac_no_outliers() {
        // Perfect linear data, no outliers.
        let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0];

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base)
            .with_random_state(42)
            .with_residual_threshold(1.0);
        let fitted = model.fit(&x, &y).unwrap();

        // All should be inliers.
        let mask = fitted.inlier_mask();
        assert!(mask.iter().all(|&v| v), "All should be inliers");

        // Predictions should be accurate.
        let preds = fitted.predict(&x).unwrap();
        for (p, &actual) in preds.iter().zip(y.iter()) {
            assert_relative_eq!(*p, actual, epsilon = 0.5);
        }
    }

    #[test]
    fn test_ransac_with_outlier() {
        // y = 2x, but one outlier.
        let x = Array2::from_shape_vec((6, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 100.0]; // outlier at idx 5

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base)
            .with_random_state(42)
            .with_max_trials(200)
            .with_residual_threshold(2.0);
        let fitted = model.fit(&x, &y).unwrap();

        let mask = fitted.inlier_mask();
        // The outlier at index 5 should be detected.
        assert!(!mask[5], "Outlier at index 5 should not be an inlier");

        // Most other points should be inliers.
        let n_inliers: usize = mask.iter().filter(|&&v| v).count();
        assert!(
            n_inliers >= 4,
            "Expected at least 4 inliers, got {n_inliers}"
        );

        // The prediction at x=3 should be close to 6.
        let x_test = Array2::from_shape_vec((1, 1), vec![3.0]).unwrap();
        let pred = fitted.predict(&x_test).unwrap();
        assert!(
            (pred[0] - 6.0).abs() < 3.0,
            "Prediction at x=3 should be near 6.0, got {}",
            pred[0]
        );
    }

    #[test]
    fn test_ransac_multiple_outliers() {
        // y = x + 1, with two outliers.
        let x =
            Array2::from_shape_vec((8, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0]).unwrap();
        let y = array![2.0, 3.0, 50.0, 5.0, 6.0, -40.0, 8.0, 9.0]; // outliers at 2 and 5

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base)
            .with_random_state(123)
            .with_max_trials(500)
            .with_residual_threshold(2.0);
        let fitted = model.fit(&x, &y).unwrap();

        let mask = fitted.inlier_mask();
        // Outliers at index 2 and 5 should be detected.
        assert!(!mask[2], "Outlier at index 2 should not be an inlier");
        assert!(!mask[5], "Outlier at index 5 should not be an inlier");
    }

    #[test]
    fn test_ransac_shape_mismatch() {
        let x = Array2::from_shape_vec((3, 1), vec![1.0, 2.0, 3.0]).unwrap();
        let y = array![1.0, 2.0];

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base);
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_ransac_insufficient_samples() {
        let x = Array2::from_shape_vec((1, 1), vec![1.0]).unwrap();
        let y = array![1.0];

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base).with_min_samples(3);
        assert!(model.fit(&x, &y).is_err());
    }

    #[test]
    fn test_ransac_reproducible_with_seed() {
        let x = Array2::from_shape_vec((6, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 100.0];

        let base1 = LinearRegression::<f64>::new();
        let model1 = RANSACRegressor::new(base1)
            .with_random_state(42)
            .with_residual_threshold(2.0);
        let fitted1 = model1.fit(&x, &y).unwrap();

        let base2 = LinearRegression::<f64>::new();
        let model2 = RANSACRegressor::new(base2)
            .with_random_state(42)
            .with_residual_threshold(2.0);
        let fitted2 = model2.fit(&x, &y).unwrap();

        // Same seed should produce same inlier mask.
        assert_eq!(fitted1.inlier_mask(), fitted2.inlier_mask());
    }

    #[test]
    fn test_ransac_auto_threshold() {
        // No explicit threshold — should use MAD.
        let x = Array2::from_shape_vec((6, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![2.0, 4.0, 6.0, 8.0, 10.0, 100.0];

        let base = LinearRegression::<f64>::new();
        let model = RANSACRegressor::new(base)
            .with_random_state(42)
            .with_max_trials(200);
        let fitted = model.fit(&x, &y).unwrap();

        let mask = fitted.inlier_mask();
        // At least some points should be inliers.
        let n_inliers: usize = mask.iter().filter(|&&v| v).count();
        assert!(
            n_inliers >= 3,
            "Expected at least 3 inliers, got {n_inliers}"
        );
    }
}