ferrolearn_linear/

isotonic.rs

//! Isotonic (monotonic) regression.
//!
//! This module provides [`IsotonicRegression`], a non-parametric regression
//! model that fits a piecewise-constant (step) function subject to a
//! monotonicity constraint. The fitted model uses linear interpolation
//! between breakpoints for prediction.
//!
//! # Algorithm
//!
//! Uses the **Pool Adjacent Violators (PAV)** algorithm, which runs in
//! `O(n)` time.
//!
//! # Examples
//!
//! ```
//! use ferrolearn_linear::isotonic::{IsotonicRegression, OutOfBounds};
//! use ferrolearn_core::{Fit, Predict};
//! use ndarray::{array, Array1, Array2};
//!
//! let model = IsotonicRegression::<f64>::new();
//! let x = Array2::from_shape_vec((5, 1), vec![1.0, 2.0, 3.0, 4.0, 5.0]).unwrap();
//! let y = array![1.0, 3.0, 2.0, 5.0, 4.0];
//!
//! let fitted = model.fit(&x, &y).unwrap();
//! let preds = fitted.predict(&x).unwrap();
//! // Predictions are monotonically non-decreasing.
//! for i in 1..preds.len() {
//!     assert!(preds[i] >= preds[i - 1]);
//! }
//! ```
//!
//! ## REQ status (per `.design/linear/isotonic.md`, mirrors `sklearn/isotonic.py` @ 1.5.2)
//!
//! | REQ | Status | Evidence |
//! |---|---|---|
//! | REQ-1 (increasing PAVA fit) | SHIPPED | `fn make_unique` → `pav_increasing_unique_weighted`; distinct-X fit matches the live oracle (`X=[1..6],y=[1,4,2,5,3,7]` → `[1,3,3,4,4,7]`). Consumer: `Fit for IsotonicRegression`. |
//! | REQ-2 (decreasing) | SHIPPED | negate-fit-negate path; decreasing dup-X `[4,2,3,1]` → `[3,3,1]` matches oracle. |
//! | REQ-3 (predict piecewise-LINEAR interpolation) | SHIPPED | `predict_single` does `y0 + t*(y1-y0)` (`scipy interp1d(kind='linear')`); test `test_interpolation`. |
//! | REQ-4 (out_of_bounds nan/clip/raise; default `nan`) | SHIPPED | `OutOfBounds::{Nan,Clip,Raise}`; `new()` defaults `Nan` (`isotonic.py:274`); test `isotonic_default_out_of_bounds_nan`. Closed #565. |
//! | REQ-8 (`_make_unique` weighted duplicate-X collapse) | SHIPPED | `fn make_unique` collapses equal-X runs to `(x, Σwy/Σw, Σw)` + weighted PAVA (`isotonic.py:308-325`); test `isotonic_make_unique_duplicate_x` (`[1,1,2,3]/[1,3,2,4]` → `[2,2,4]`). Closed #569. |
//! | REQ-5 (y_min/y_max clipping) | SHIPPED | `IsotonicRegression` gains `pub y_min`/`pub y_max: Option<F>` fields (default `None` in `new`, matching `isotonic.py:274`) + `#[must_use] with_y_min`/`with_y_max` builders; `fn fit_with_sample_weight` clips each pooled `y_threshold` to `[y_min.unwrap_or(-inf), y_max.unwrap_or(+inf)]` AFTER PAVA (and after the decreasing negate-fit-negate is undone), mirroring `np.clip(y, y_min, y_max, y)` (`isotonic.py:163-170`). Both-`None` is a no-op (byte-identical unclipped path). Consumer: `Fit::fit` → `FittedIsotonicRegression` (crate-root export). Test: `isotonic_y_min_y_max` (divergence suite, live oracle `y_min=2`→`[2,2,3,4,5]`, `y_max=4`→`[1,2,3,4,4]`, both→`[2,2,3,4,4]`). #566. |
//! | REQ-6 (increasing='auto' via Spearman) | SHIPPED | `enum Increasing::Auto` + `fn with_increasing_auto`/`fn with_increasing_mode`; `fn fit_with_sample_weight` resolves `Auto` via the free `fn check_increasing` (Spearman rho sign, `sklearn/isotonic.py:32-98,306-307`) and stores the bool in `FittedIsotonicRegression::increasing`. Consumer: `Fit::fit`. Test: `isotonic_increasing_auto` (divergence suite, live oracle `X=[1..4],y=[4,3,2,1]`→decreasing, `increasing_==false`). #567. |
//! | REQ-7 (sample_weight public API) | SHIPPED | `fn fit_with_sample_weight` threads per-sample weights into weighted `make_unique` (weighted-mean collapse) + `pav_increasing_unique_weighted` (weighted pool), mirroring `IsotonicRegression.fit(X,y,sample_weight)` → `_build_y` `_make_unique`/`isotonic_regression` (`isotonic.py:251`,`:300-328`). Consumer: `Fit::fit` delegates with an all-ones weight vector. Test: `isotonic_sample_weight` (divergence suite). Closed #568. |
//! | REQ-9 (X_min_/X_max_/X_thresholds_/y_thresholds_/increasing_) | SHIPPED | `FittedIsotonicRegression::{x_min,x_max,x_thresholds,y_thresholds,increasing}` accessors mirror `X_min_`/`X_max_`/`X_thresholds_`/`y_thresholds_`/`increasing_` (`sklearn/isotonic.py:331,393,307-309`); `fit_with_sample_weight` applies sklearn's `trim_duplicates` interior-plateau trim (`isotonic.py:333-341`) to the stored thresholds. Consumer: `Fit::fit` → these accessors are read by the predict path (`x_min`/`x_max` bound the interpolant). Test: `isotonic_fitted_attributes` (live oracle `X=[1..4],y=[1,3,2,4]`→`x_min=1,x_max=4,x_thr=[1,2,3,4],y_thr=[1,2.5,2.5,4],increasing=true`). #570. |
//! | REQ-10 (free `isotonic_regression` + `check_increasing`) | SHIPPED | `pub fn check_increasing` (Spearman rho sign, `isotonic.py:32-98`) + `pub fn isotonic_regression` (free PAVA with `sample_weight`/`y_min`/`y_max`/`increasing`, `isotonic.py:111-171`). Consumer: `check_increasing` consumed by `fit_with_sample_weight`'s `Auto` resolution; `isotonic_regression` reuses the internal weighted-PAVA machinery and is itself a production free function. Tests: `isotonic_free_check_increasing`, `isotonic_free_isotonic_regression` (live oracle). #571. |
//! | REQ-11 (ferray substrate) | NOT-STARTED | #572 (crate-wide-deferred, cf. ridge #391). |

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

// ---------------------------------------------------------------------------
// Out-of-bounds strategy
// ---------------------------------------------------------------------------

/// Strategy for handling predictions outside the training range.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum OutOfBounds {
    /// Clip predictions to the range of training values.
    Clip,
    /// Return NaN for out-of-range inputs.
    Nan,
    /// Return an error for out-of-range inputs.
    Raise,
}

// ---------------------------------------------------------------------------
// Increasing mode
// ---------------------------------------------------------------------------

/// Monotonicity direction for the fitted function.
///
/// Mirrors scikit-learn's `increasing` constructor parameter, whose
/// `_parameter_constraints` allows `["boolean", StrOptions({"auto"})]` with
/// default `True` (`sklearn/isotonic.py:271-274`):
///
/// - [`Increasing::True`] — force a non-decreasing fit (`increasing=True`).
/// - [`Increasing::False`] — force a non-increasing fit (`increasing=False`).
/// - [`Increasing::Auto`] — resolve the direction from the data at fit time via
///   a Spearman correlation test (`increasing='auto'`,
///   `sklearn/isotonic.py:306-307`: `self.increasing_ = check_increasing(X, y)`).
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum Increasing {
    /// Force a non-decreasing (increasing) fit. The default, matching
    /// `IsotonicRegression(increasing=True)` (`sklearn/isotonic.py:274`).
    #[default]
    True,
    /// Force a non-increasing (decreasing) fit.
    False,
    /// Resolve the direction from the data via a Spearman correlation test
    /// (`increasing='auto'`).
    Auto,
}

impl From<bool> for Increasing {
    /// `true` → [`Increasing::True`], `false` → [`Increasing::False`].
    ///
    /// This preserves the prior `with_increasing(bool)` API semantics.
    fn from(b: bool) -> Self {
        if b {
            Increasing::True
        } else {
            Increasing::False
        }
    }
}

// ---------------------------------------------------------------------------
// IsotonicRegression (unfitted)
// ---------------------------------------------------------------------------

/// Isotonic regression configuration.
///
/// Fits a piecewise-constant monotonic function using the Pool Adjacent
/// Violators (PAV) algorithm.
///
/// # Type Parameters
///
/// - `F`: The floating-point type (`f32` or `f64`).
#[derive(Debug, Clone)]
pub struct IsotonicRegression<F> {
    /// Monotonicity direction: increasing, decreasing, or auto-resolved from
    /// the data via a Spearman test. Mirrors scikit-learn's `increasing`
    /// constructor parameter (`sklearn/isotonic.py:271-274`, default `True`).
    pub increasing: Increasing,
    /// Strategy for predictions outside the training range.
    pub out_of_bounds: OutOfBounds,
    /// Lower bound on the lowest predicted value. `None` (the default) means
    /// `-inf` — no lower clip. Mirrors scikit-learn's `y_min=None`
    /// (`sklearn/isotonic.py:274`); the pooled `y_thresholds` are clipped to
    /// `[y_min, y_max]` after PAVA (`isotonic.py:163-170`).
    pub y_min: Option<F>,
    /// Upper bound on the highest predicted value. `None` (the default) means
    /// `+inf` — no upper clip. Mirrors scikit-learn's `y_max=None`
    /// (`sklearn/isotonic.py:274`); the pooled `y_thresholds` are clipped to
    /// `[y_min, y_max]` after PAVA (`isotonic.py:163-170`).
    pub y_max: Option<F>,
    _marker: std::marker::PhantomData<F>,
}

impl<F: Float + Send + Sync + 'static> IsotonicRegression<F> {
    /// Fit the isotonic regression model with per-sample weights.
    ///
    /// This is the weighted generalization of [`Fit::fit`]. Each sample
    /// `(x[i], y[i])` carries weight `sample_weight[i]`; the weights flow into
    /// the `_make_unique` duplicate-`X` collapse (each equal-`X` run collapses
    /// to its sample-weighted mean `Σ wᵢ yᵢ / Σ wᵢ` and summed weight) and into
    /// the weighted PAV pool, mirroring scikit-learn's
    /// `IsotonicRegression.fit(X, y, sample_weight)` → `_build_y`
    /// (`sklearn/isotonic.py:300-328`, the `_make_unique` + `isotonic_regression`
    /// weighted pipeline) at tag 1.5.2.
    ///
    /// Zero-weight samples are removed before fitting, matching
    /// `_build_y`'s `mask = sample_weight > 0` filter (`isotonic.py:314-315`).
    ///
    /// [`Fit::fit`] is exactly this method with an all-ones weight vector, so
    /// the default (unweighted) path is byte-identical.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` does not have exactly one
    /// column, or if `y`/`sample_weight` lengths do not match the sample count.
    /// Returns [`FerroError::InvalidParameter`] if any weight is negative
    /// (mirroring sklearn's `_check_sample_weight` non-negativity contract).
    /// Returns [`FerroError::InsufficientSamples`] if fewer than 2 positively
    /// weighted samples remain.
    pub fn fit_with_sample_weight(
        &self,
        x: &Array2<F>,
        y: &Array1<F>,
        sample_weight: &Array1<F>,
    ) -> Result<FittedIsotonicRegression<F>, FerroError> {
        let (n_samples, n_features) = x.dim();

        if n_features != 1 {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples, 1],
                actual: vec![n_samples, n_features],
                context: "IsotonicRegression requires exactly 1 feature".into(),
            });
        }

        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(),
            });
        }

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

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

        // Non-negativity: sklearn's `_check_sample_weight` rejects negative
        // weights (and `_build_y` then drops the zero-weight rows).
        if sample_weight.iter().any(|&w| w < F::zero()) {
            return Err(FerroError::InvalidParameter {
                name: "sample_weight".into(),
                reason: "sample weights must be non-negative".into(),
            });
        }

        // Extract the single feature column, dropping zero-weight rows
        // (`isotonic.py:314-315`: `mask = sample_weight > 0`).
        let mut xs: Vec<F> = Vec::with_capacity(n_samples);
        let mut ys: Vec<F> = Vec::with_capacity(n_samples);
        let mut ws: Vec<F> = Vec::with_capacity(n_samples);
        let col = x.column(0);
        for i in 0..n_samples {
            if sample_weight[i] > F::zero() {
                xs.push(col[i]);
                ys.push(y[i]);
                ws.push(sample_weight[i]);
            }
        }

        if xs.len() < 2 {
            return Err(FerroError::InsufficientSamples {
                required: 2,
                actual: xs.len(),
                context: "IsotonicRegression requires at least 2 positively weighted samples"
                    .into(),
            });
        }

        // Resolve the monotonicity direction. For `Increasing::Auto` this runs
        // a Spearman correlation test over the (positively weighted) sample
        // `(x, y)` pairs, mirroring scikit-learn's `_build_y`
        // (`sklearn/isotonic.py:306-307`: `if self.increasing == "auto":
        // self.increasing_ = check_increasing(X, y)`). NOTE: sklearn resolves
        // BEFORE the zero-weight `mask` filter (`isotonic.py:306` precedes
        // `:314-315`), so `check_increasing` sees all rows; we replicate that by
        // resolving on the full `x.column(0)` / `y` rather than the filtered
        // `xs`/`ys`.
        let increasing: bool = match self.increasing {
            Increasing::True => true,
            Increasing::False => false,
            Increasing::Auto => {
                let x_full: Vec<F> = col.to_vec();
                let y_full: Vec<F> = y.to_vec();
                check_increasing(&x_full, &y_full)
            }
        };

        let (mut result_x, mut result_y) = if increasing {
            let (ux, uy, uw) = make_unique(&xs, &ys, &ws);
            pav_increasing_unique_weighted(&ux, &uy, &uw)
        } else {
            // For decreasing: negate y, run weighted increasing PAV, negate
            // result — threading the same per-sample weights through.
            let neg_ys: Vec<F> = ys.iter().map(|&v| -v).collect();
            let (ux, uy, uw) = make_unique(&xs, &neg_ys, &ws);
            let (rx, ry) = pav_increasing_unique_weighted(&ux, &uy, &uw);
            let ry_neg: Vec<F> = ry.iter().map(|&v| -v).collect();
            (rx, ry_neg)
        };

        // Clip the pooled `y_thresholds` to `[y_min, y_max]` AFTER PAVA (and
        // after the decreasing negate-fit-negate is undone), mirroring
        // scikit-learn's `np.clip(y, y_min, y_max, y)` on the pooled values
        // (`sklearn/isotonic.py:163-170`). Unset bounds default to the open
        // `±inf` (`isotonic.py:165-168`), so when both `y_min`/`y_max` are
        // `None` this is `y.max(-inf).min(+inf)` — a no-op leaving every
        // threshold byte-identical to the unclipped path. The clip is applied
        // to the STORED thresholds so `predict` (linear interpolation between
        // them) stays within `[y_min, y_max]`.
        if self.y_min.is_some() || self.y_max.is_some() {
            let lo = self.y_min.unwrap_or_else(F::neg_infinity);
            let hi = self.y_max.unwrap_or_else(F::infinity);
            for y in &mut result_y {
                *y = y.max(lo).min(hi);
            }
        }

        // Trim interior plateau points so the stored thresholds mirror
        // scikit-learn's `X_thresholds_`/`y_thresholds_` exactly: aside from the
        // first and last point, drop any point whose `y` equals BOTH its
        // neighbors (`sklearn/isotonic.py:333-341`, the `trim_duplicates`
        // branch: `keep_data[1:-1] = not_equal(y[1:-1], y[:-2]) | not_equal(
        // y[1:-1], y[2:])`). This is purely a storage compaction — the
        // piecewise-linear interpolant is unchanged because the dropped points
        // lie on a flat segment between two retained breakpoints with the same
        // `y`.
        if result_y.len() > 2 {
            let n = result_y.len();
            let mut kept_x = Vec::with_capacity(n);
            let mut kept_y = Vec::with_capacity(n);
            kept_x.push(result_x[0]);
            kept_y.push(result_y[0]);
            for i in 1..n - 1 {
                if result_y[i] != result_y[i - 1] || result_y[i] != result_y[i + 1] {
                    kept_x.push(result_x[i]);
                    kept_y.push(result_y[i]);
                }
            }
            kept_x.push(result_x[n - 1]);
            kept_y.push(result_y[n - 1]);
            result_x = kept_x;
            result_y = kept_y;
        }

        // Ensure at least 2 breakpoints.
        if result_x.len() < 2 {
            // All same x value: duplicate.
            if result_x.len() == 1 {
                result_x.push(result_x[0]);
                result_y.push(result_y[0]);
            } else {
                return Err(FerroError::NumericalInstability {
                    message: "PAV produced no breakpoints".into(),
                });
            }
        }

        Ok(FittedIsotonicRegression {
            x_thresholds: result_x,
            y_thresholds: result_y,
            out_of_bounds: self.out_of_bounds,
            increasing,
        })
    }
}

impl<F: Float> IsotonicRegression<F> {
    /// Create a new `IsotonicRegression` with default settings.
    ///
    /// Defaults: `increasing = true`, `out_of_bounds = Nan`, `y_min = None`,
    /// `y_max = None`.
    ///
    /// The `out_of_bounds` default matches scikit-learn's
    /// `IsotonicRegression(out_of_bounds="nan")` (`sklearn/isotonic.py:274`):
    /// a default-constructed estimator returns `NaN` for predictions outside
    /// the training range `[X_min_, X_max_]`.
    ///
    /// The `y_min`/`y_max` defaults of `None` match scikit-learn's
    /// `IsotonicRegression(y_min=None, y_max=None)` (`sklearn/isotonic.py:274`):
    /// with both unset the pooled `y_thresholds` are clipped to
    /// `[-inf, +inf]`, i.e. not clipped at all.
    #[must_use]
    pub fn new() -> Self {
        Self {
            increasing: Increasing::True,
            out_of_bounds: OutOfBounds::Nan,
            y_min: None,
            y_max: None,
            _marker: std::marker::PhantomData,
        }
    }

    /// Set the lower bound on the lowest predicted value (`y_min`).
    ///
    /// The pooled `y_thresholds` produced by PAVA are clipped so none falls
    /// below `y_min`, mirroring scikit-learn's `np.clip(y, y_min, y_max, y)`
    /// after pooling (`sklearn/isotonic.py:163-170`; constructor `y_min`,
    /// `isotonic.py:274`).
    #[must_use]
    pub fn with_y_min(mut self, y_min: F) -> Self {
        self.y_min = Some(y_min);
        self
    }

    /// Set the upper bound on the highest predicted value (`y_max`).
    ///
    /// The pooled `y_thresholds` produced by PAVA are clipped so none rises
    /// above `y_max`, mirroring scikit-learn's `np.clip(y, y_min, y_max, y)`
    /// after pooling (`sklearn/isotonic.py:163-170`; constructor `y_max`,
    /// `isotonic.py:274`).
    #[must_use]
    pub fn with_y_max(mut self, y_max: F) -> Self {
        self.y_max = Some(y_max);
        self
    }

    /// Set whether the fitted function should be increasing.
    ///
    /// `true` → [`Increasing::True`] (non-decreasing), `false` →
    /// [`Increasing::False`] (non-increasing). This preserves the prior
    /// `with_increasing(bool)` API; for the data-resolved `'auto'` direction use
    /// [`with_increasing_auto`](Self::with_increasing_auto) or
    /// [`with_increasing_mode`](Self::with_increasing_mode).
    #[must_use]
    pub fn with_increasing(mut self, increasing: bool) -> Self {
        self.increasing = Increasing::from(increasing);
        self
    }

    /// Resolve the monotonicity direction from the data at fit time via a
    /// Spearman correlation test, mirroring scikit-learn's
    /// `IsotonicRegression(increasing='auto')` (`sklearn/isotonic.py:306-307`).
    ///
    /// The resolved direction is exposed by
    /// [`FittedIsotonicRegression::increasing`].
    #[must_use]
    pub fn with_increasing_auto(mut self) -> Self {
        self.increasing = Increasing::Auto;
        self
    }

    /// Set the monotonicity direction directly via the [`Increasing`] enum.
    ///
    /// Mirrors scikit-learn's `increasing` parameter
    /// (`sklearn/isotonic.py:271-274`), which accepts `True`/`False`/`'auto'`.
    #[must_use]
    pub fn with_increasing_mode(mut self, increasing: Increasing) -> Self {
        self.increasing = increasing;
        self
    }

    /// Set the out-of-bounds strategy.
    #[must_use]
    pub fn with_out_of_bounds(mut self, strategy: OutOfBounds) -> Self {
        self.out_of_bounds = strategy;
        self
    }
}

impl<F: Float> Default for IsotonicRegression<F> {
    fn default() -> Self {
        Self::new()
    }
}

// ---------------------------------------------------------------------------
// FittedIsotonicRegression
// ---------------------------------------------------------------------------

/// Fitted isotonic regression model.
///
/// Stores the breakpoints of the fitted step function and uses linear
/// interpolation between them for prediction.
#[derive(Debug, Clone)]
pub struct FittedIsotonicRegression<F> {
    /// Sorted x-values of breakpoints.
    x_thresholds: Vec<F>,
    /// Corresponding y-values (monotonic).
    y_thresholds: Vec<F>,
    /// Out-of-bounds strategy.
    out_of_bounds: OutOfBounds,
    /// Whether the function is increasing.
    increasing: bool,
}

impl<F: Float> FittedIsotonicRegression<F> {
    /// Returns whether the fitted function is increasing.
    #[must_use]
    pub fn is_increasing(&self) -> bool {
        self.increasing
    }

    /// The resolved monotonicity direction (`true` = increasing).
    ///
    /// Mirrors scikit-learn's fitted `increasing_` attribute
    /// (`sklearn/isotonic.py:307-309`). When the estimator was configured with
    /// [`Increasing::Auto`] this is the direction resolved from the data via the
    /// Spearman test; otherwise it equals the requested direction.
    #[must_use]
    pub fn increasing(&self) -> bool {
        self.increasing
    }

    /// The minimum training `X` value (`X_min_`).
    ///
    /// Mirrors scikit-learn's fitted `X_min_` attribute
    /// (`sklearn/isotonic.py:331`: `self.X_min_, self.X_max_ = np.min(X),
    /// np.max(X)`). The thresholds are sorted ascending, so this is the first
    /// stored threshold.
    #[must_use]
    pub fn x_min(&self) -> F {
        self.x_thresholds[0]
    }

    /// The maximum training `X` value (`X_max_`).
    ///
    /// Mirrors scikit-learn's fitted `X_max_` attribute
    /// (`sklearn/isotonic.py:331`). The thresholds are sorted ascending, so this
    /// is the last stored threshold.
    #[must_use]
    pub fn x_max(&self) -> F {
        self.x_thresholds[self.x_thresholds.len() - 1]
    }

    /// The breakpoint `X` values of the fitted step function (`X_thresholds_`).
    ///
    /// Mirrors scikit-learn's fitted `X_thresholds_` attribute
    /// (`sklearn/isotonic.py:393`), after the interior-plateau trim
    /// (`isotonic.py:333-341`).
    #[must_use]
    pub fn x_thresholds(&self) -> &[F] {
        &self.x_thresholds
    }

    /// The breakpoint `y` values of the fitted step function (`y_thresholds_`),
    /// monotonic in the resolved direction.
    ///
    /// Mirrors scikit-learn's fitted `y_thresholds_` attribute
    /// (`sklearn/isotonic.py:393`), after the interior-plateau trim
    /// (`isotonic.py:333-341`).
    #[must_use]
    pub fn y_thresholds(&self) -> &[F] {
        &self.y_thresholds
    }

    /// Predict a single value using linear interpolation.
    fn predict_single(&self, x: F) -> Result<F, FerroError> {
        if self.x_thresholds.is_empty() {
            return Err(FerroError::NumericalInstability {
                message: "isotonic model has no breakpoints".into(),
            });
        }

        let x_min = self.x_thresholds[0];
        let x_max = *self.x_thresholds.last().unwrap();

        if x < x_min {
            return match self.out_of_bounds {
                OutOfBounds::Clip => Ok(self.y_thresholds[0]),
                OutOfBounds::Nan => Ok(F::nan()),
                OutOfBounds::Raise => Err(FerroError::InvalidParameter {
                    name: "x".into(),
                    reason: "value is below training range".into(),
                }),
            };
        }

        if x > x_max {
            return match self.out_of_bounds {
                OutOfBounds::Clip => Ok(*self.y_thresholds.last().unwrap()),
                OutOfBounds::Nan => Ok(F::nan()),
                OutOfBounds::Raise => Err(FerroError::InvalidParameter {
                    name: "x".into(),
                    reason: "value is above training range".into(),
                }),
            };
        }

        // Binary search for the interval containing x.
        let n = self.x_thresholds.len();

        // Handle exact match at the last point.
        if x == x_max {
            return Ok(*self.y_thresholds.last().unwrap());
        }

        // Find the interval [x_thresholds[i], x_thresholds[i+1]) containing x.
        let mut lo = 0;
        let mut hi = n - 1;
        while lo < hi - 1 {
            let mid = usize::midpoint(lo, hi);
            if self.x_thresholds[mid] <= x {
                lo = mid;
            } else {
                hi = mid;
            }
        }

        let x0 = self.x_thresholds[lo];
        let x1 = self.x_thresholds[hi];
        let y0 = self.y_thresholds[lo];
        let y1 = self.y_thresholds[hi];

        if (x1 - x0).abs() < F::epsilon() {
            return Ok(y0);
        }

        // Linear interpolation.
        let t = (x - x0) / (x1 - x0);
        Ok(y0 + t * (y1 - y0))
    }
}

// ---------------------------------------------------------------------------
// Pool Adjacent Violators (PAV) algorithm
// ---------------------------------------------------------------------------

/// Collapse maximal runs of equal `X` into a single point, mirroring
/// scikit-learn's `_make_unique` (`sklearn/_isotonic.pyx`).
///
/// The inputs are first ordered by `X` (ties broken by `y`, matching
/// `np.lexsort((y, X))` at `sklearn/isotonic.py:317`). Each run of equal `X`
/// then collapses to one point whose `x` is the shared value, whose `y` is the
/// **sample-weight-weighted mean** of the run (`Σ wᵢ yᵢ / Σ wᵢ`), and whose
/// weight is the **summed** weight of the run.
///
/// For unit weights this reduces to the plain mean and a count, so the
/// returned weights double as the run multiplicities consumed by the weighted
/// PAVA. Returns `(x_unique, y_unique, w_unique)`.
fn make_unique<F: Float>(xs: &[F], ys: &[F], ws: &[F]) -> (Vec<F>, Vec<F>, Vec<F>) {
    let n = xs.len();
    if n == 0 {
        return (Vec::new(), Vec::new(), Vec::new());
    }

    // Order by X (primary), y (secondary) — np.lexsort((y, X)). total_cmp
    // gives a total order without panicking on NaN (goal.md R-APG-1).
    let mut indices: Vec<usize> = (0..n).collect();
    indices.sort_by(|&a, &b| {
        xs[a]
            .partial_cmp(&xs[b])
            .unwrap_or(std::cmp::Ordering::Equal)
            .then_with(|| {
                ys[a]
                    .partial_cmp(&ys[b])
                    .unwrap_or(std::cmp::Ordering::Equal)
            })
    });

    let mut x_out = Vec::new();
    let mut y_out = Vec::new();
    let mut w_out = Vec::new();

    let mut cur_x = xs[indices[0]];
    let mut cur_w = F::zero();
    let mut cur_wy = F::zero();

    for &idx in &indices {
        let x = xs[idx];
        let w = ws[idx];
        if x != cur_x {
            // Close the previous run.
            x_out.push(cur_x);
            w_out.push(cur_w);
            y_out.push(cur_wy / cur_w);

            cur_x = x;
            cur_w = w;
            cur_wy = ys[idx] * w;
        } else {
            cur_w = cur_w + w;
            cur_wy = cur_wy + ys[idx] * w;
        }
    }
    // Close the final run.
    x_out.push(cur_x);
    w_out.push(cur_w);
    y_out.push(cur_wy / cur_w);

    (x_out, y_out, w_out)
}

/// Run the **weighted** PAV algorithm on points pre-ordered and de-duplicated
/// by `X` (see [`make_unique`]), producing a monotonically non-decreasing set
/// of `(x, y)` breakpoints.
///
/// When two adjacent blocks violate monotonicity they are pooled: the merged
/// block's value is the weighted mean `(w₁·v₁ + w₂·v₂)/(w₁ + w₂)` and its
/// weight is `w₁ + w₂`, mirroring sklearn's
/// `_inplace_contiguous_isotonic_regression` (`sklearn/_isotonic.pyx`). The
/// `xs`/`ys`/`ws` slices must already be sorted by `x` with unique `x` values.
fn pav_increasing_unique_weighted<F: Float>(xs: &[F], ys: &[F], ws: &[F]) -> (Vec<F>, Vec<F>) {
    let n = xs.len();

    // PAV: merge adjacent blocks that violate monotonicity.
    // Each block carries the weighted sum, total weight, and x extent.
    struct Block<F> {
        wsum: F,
        weight: F,
        first_idx: usize,
        last_idx: usize,
    }

    let mut blocks: Vec<Block<F>> = Vec::with_capacity(n);

    for i in 0..n {
        blocks.push(Block {
            wsum: ys[i] * ws[i],
            weight: ws[i],
            first_idx: i,
            last_idx: i,
        });

        // Merge with previous blocks as needed.
        while blocks.len() > 1 {
            let len = blocks.len();
            let prev_mean = blocks[len - 2].wsum / blocks[len - 2].weight;
            let curr_mean = blocks[len - 1].wsum / blocks[len - 1].weight;

            if prev_mean > curr_mean {
                // Pool the two violating blocks.
                let Some(last) = blocks.pop() else { break };
                let Some(prev) = blocks.last_mut() else { break };
                prev.wsum = prev.wsum + last.wsum;
                prev.weight = prev.weight + last.weight;
                prev.last_idx = last.last_idx;
            } else {
                break;
            }
        }
    }

    // Extract breakpoints: for each block, emit the first and (if distinct)
    // last x at the pooled weighted mean.
    let mut result_x = Vec::new();
    let mut result_y = Vec::new();

    for block in &blocks {
        let mean = block.wsum / block.weight;
        let bx0 = xs[block.first_idx];
        let bx1 = xs[block.last_idx];

        if result_x.is_empty() || result_x.last().is_none_or(|&last| last != bx0) {
            result_x.push(bx0);
            result_y.push(mean);
        }
        if bx0 != bx1 {
            result_x.push(bx1);
            result_y.push(mean);
        }
    }

    (result_x, result_y)
}

// ---------------------------------------------------------------------------
// Free functions: check_increasing / isotonic_regression
// ---------------------------------------------------------------------------

/// Average (fractional) ranks of `v`, ties resolved to the mean rank of the
/// tied group — the rank convention `scipy.stats.spearmanr` uses internally
/// (`scipy.stats.rankdata` with `method='average'`).
///
/// Returned ranks are 1-based (rank 1 = smallest), matching `rankdata`. NaN is
/// ordered as greater-than-all via `total_cmp`-style fallback so the routine
/// never panics (goal.md R-APG-1).
fn average_ranks<F: Float>(v: &[F]) -> Vec<F> {
    let n = v.len();
    let mut idx: Vec<usize> = (0..n).collect();
    idx.sort_by(|&a, &b| v[a].partial_cmp(&v[b]).unwrap_or(std::cmp::Ordering::Equal));

    let mut ranks = vec![F::zero(); n];
    let mut i = 0;
    while i < n {
        // Find the extent of the tied group [i, j).
        let mut j = i + 1;
        while j < n && v[idx[j]] == v[idx[i]] {
            j += 1;
        }
        // Average of the 1-based positions i+1 .. j is (i + j + 1) / 2.
        let count = j - i;
        let sum_pos = {
            // Σ_{k=i}^{j-1} (k + 1) = count*(i+1) + (0+1+...+(count-1)).
            let mut s = F::zero();
            for k in 0..count {
                s = s + F::from(i + 1 + k).unwrap_or_else(F::zero);
            }
            s
        };
        let avg = sum_pos / F::from(count).unwrap_or_else(F::one);
        for &orig in &idx[i..j] {
            ranks[orig] = avg;
        }
        i = j;
    }
    ranks
}

/// Determine whether `y` is monotonically increasing or decreasing with respect
/// to `x`, via the sign of the Spearman rank correlation.
///
/// This is the free function `sklearn.isotonic.check_increasing`
/// (`sklearn/isotonic.py:32-98`): it computes the Spearman rho between `x` and
/// `y` and returns `rho >= 0` (`isotonic.py:76-77`: `rho, _ = spearmanr(x, y);
/// increasing_bool = rho >= 0`). The Spearman rho is the Pearson correlation of
/// the average ranks of `x` and `y`.
///
/// scikit-learn additionally emits a `UserWarning` when the 95% Fisher-transform
/// confidence interval of rho spans zero (`isotonic.py:79-96`). That branch is
/// purely advisory (it does not change the returned bool), so it is intentionally
/// omitted here — the contract is the returned direction.
///
/// Degenerate inputs return `true` (sklearn's `rho` is `NaN` for a constant
/// input, and `np.nan >= 0` is `False` in numpy — but for empty/constant data
/// the direction is conventionally treated as increasing; this only affects
/// inputs that PAVA handles identically in either direction).
#[must_use]
pub fn check_increasing<F: Float + Send + Sync + 'static>(x: &[F], y: &[F]) -> bool {
    let n = x.len();
    if n == 0 || n != y.len() {
        return true;
    }

    let rx = average_ranks(x);
    let ry = average_ranks(y);

    // Pearson correlation of the ranks.
    let nf = F::from(n).unwrap_or_else(F::one);
    let mean_x = rx.iter().fold(F::zero(), |a, &v| a + v) / nf;
    let mean_y = ry.iter().fold(F::zero(), |a, &v| a + v) / nf;

    let mut cov = F::zero();
    let mut var_x = F::zero();
    let mut var_y = F::zero();
    for i in 0..n {
        let dx = rx[i] - mean_x;
        let dy = ry[i] - mean_y;
        cov = cov + dx * dy;
        var_x = var_x + dx * dx;
        var_y = var_y + dy * dy;
    }

    // Constant ranks (no variance): rho is undefined; treat as increasing.
    if var_x <= F::zero() || var_y <= F::zero() {
        return true;
    }

    let rho = cov / (var_x.sqrt() * var_y.sqrt());
    rho >= F::zero()
}

/// Solve the isotonic regression model on the sequence `y` (the free function
/// `sklearn.isotonic.isotonic_regression`, `sklearn/isotonic.py:111-171`).
///
/// Unlike the [`IsotonicRegression`] estimator, this operates purely on the
/// **order of `y`** (there is no `X` and no `_make_unique` collapse): index `i`
/// precedes index `i+1`. For `increasing = false` the sequence is reversed,
/// pooled increasing, then reversed back (`isotonic.py:156,158,170`: `order =
/// np.s_[::-1]`). Optional per-element `sample_weight` weights the pool
/// (defaults to unit weight, `isotonic.py:159`); `y_min`/`y_max` clip the pooled
/// result to `[y_min, y_max]` afterward (`isotonic.py:163-170`, unset bounds
/// default to `∓inf`).
///
/// Returns the isotonic fit `y_` in the original index order.
#[must_use]
pub fn isotonic_regression<F: Float + Send + Sync + 'static>(
    y: &[F],
    sample_weight: Option<&[F]>,
    y_min: Option<F>,
    y_max: Option<F>,
    increasing: bool,
) -> Vec<F> {
    let n = y.len();
    if n == 0 {
        return Vec::new();
    }

    // Build the working sequence in pool order (`np.s_[:]` vs `np.s_[::-1]`).
    let mut vals: Vec<F> = Vec::with_capacity(n);
    let mut wts: Vec<F> = Vec::with_capacity(n);
    for i in 0..n {
        let src = if increasing { i } else { n - 1 - i };
        vals.push(y[src]);
        wts.push(match sample_weight {
            Some(sw) if sw.len() == n => sw[src],
            _ => F::one(),
        });
    }

    // Weighted PAV on the contiguous sequence (no X de-duplication): mirrors
    // `_inplace_contiguous_isotonic_regression`. Each block carries its weighted
    // sum, total weight, and the count of original elements it spans.
    struct Block<F> {
        wsum: F,
        weight: F,
        len: usize,
    }
    let mut blocks: Vec<Block<F>> = Vec::with_capacity(n);
    for i in 0..n {
        blocks.push(Block {
            wsum: vals[i] * wts[i],
            weight: wts[i],
            len: 1,
        });
        while blocks.len() > 1 {
            let k = blocks.len();
            let prev_mean = blocks[k - 2].wsum / blocks[k - 2].weight;
            let curr_mean = blocks[k - 1].wsum / blocks[k - 1].weight;
            if prev_mean > curr_mean {
                let Some(last) = blocks.pop() else { break };
                let Some(prev) = blocks.last_mut() else { break };
                prev.wsum = prev.wsum + last.wsum;
                prev.weight = prev.weight + last.weight;
                prev.len += last.len;
            } else {
                break;
            }
        }
    }

    // Expand the pooled block means back to per-element values (in pool order).
    let lo = y_min.unwrap_or_else(F::neg_infinity);
    let hi = y_max.unwrap_or_else(F::infinity);
    let clip = y_min.is_some() || y_max.is_some();

    let mut pooled: Vec<F> = Vec::with_capacity(n);
    for block in &blocks {
        let mut mean = block.wsum / block.weight;
        if clip {
            mean = mean.max(lo).min(hi);
        }
        for _ in 0..block.len {
            pooled.push(mean);
        }
    }

    // Undo the reversal so the result is in the original index order
    // (`isotonic.py:170`: `return y[order]`).
    if increasing {
        pooled
    } else {
        pooled.into_iter().rev().collect()
    }
}

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

impl<F: Float + Send + Sync + 'static> Fit<Array2<F>, Array1<F>> for IsotonicRegression<F> {
    type Fitted = FittedIsotonicRegression<F>;
    type Error = FerroError;

    /// Fit the isotonic regression model using PAV (equal sample weights).
    ///
    /// The input `x` must have exactly one column (univariate regression).
    ///
    /// This delegates to [`IsotonicRegression::fit_with_sample_weight`] with an
    /// all-ones weight vector. With unit weights no row is dropped (none has
    /// zero weight) and the weighted `make_unique`/PAV reduce to the plain-mean
    /// special case, so this path is byte-identical to the prior unweighted
    /// implementation.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` and `y` have different
    /// sample counts or if `x` does not have exactly one column.
    /// Returns [`FerroError::InsufficientSamples`] if there are fewer than
    /// 2 samples.
    fn fit(&self, x: &Array2<F>, y: &Array1<F>) -> Result<FittedIsotonicRegression<F>, FerroError> {
        let sample_weight = Array1::<F>::from_elem(y.len(), F::one());
        self.fit_with_sample_weight(x, y, &sample_weight)
    }
}

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

    /// Predict target values for the given feature matrix.
    ///
    /// Uses linear interpolation between breakpoints.
    ///
    /// # Errors
    ///
    /// Returns [`FerroError::ShapeMismatch`] if `x` does not have exactly
    /// one column.
    /// Returns [`FerroError::InvalidParameter`] if `out_of_bounds = Raise`
    /// and a value is outside the training range.
    fn predict(&self, x: &Array2<F>) -> Result<Array1<F>, FerroError> {
        let (n_samples, n_features) = x.dim();

        if n_features != 1 {
            return Err(FerroError::ShapeMismatch {
                expected: vec![n_samples, 1],
                actual: vec![n_samples, n_features],
                context: "IsotonicRegression requires exactly 1 feature".into(),
            });
        }

        let mut result = Array1::<F>::zeros(n_samples);
        for i in 0..n_samples {
            result[i] = self.predict_single(x[[i, 0]])?;
        }

        Ok(result)
    }
}

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

    #[test]
    fn test_monotonicity_increasing() {
        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![1.0, 4.0, 2.0, 5.0, 3.0, 7.0];

        let model = IsotonicRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        // Check monotonicity: each prediction should be >= the previous.
        for i in 1..preds.len() {
            assert!(
                preds[i] >= preds[i - 1] - 1e-10,
                "Monotonicity violated at index {i}: {} < {}",
                preds[i],
                preds[i - 1]
            );
        }
    }

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

        let model = IsotonicRegression::<f64>::new().with_increasing(false);
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        // Check monotonicity: each prediction should be <= the previous.
        for i in 1..preds.len() {
            assert!(
                preds[i] <= preds[i - 1] + 1e-10,
                "Decreasing monotonicity violated at index {i}: {} > {}",
                preds[i],
                preds[i - 1]
            );
        }
    }

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

        let model = IsotonicRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        for i in 0..4 {
            assert_relative_eq!(preds[i], y[i], epsilon = 1e-10);
        }
    }

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

        let model = IsotonicRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        // Predict at intermediate points.
        let x_new = Array2::from_shape_vec((3, 1), vec![2.0, 3.0, 4.0]).unwrap();
        let preds = fitted.predict(&x_new).unwrap();

        // Linear interpolation: at x=2, y should be 2.0; at x=4, y should be 4.0.
        assert_relative_eq!(preds[0], 2.0, epsilon = 1e-10);
        assert_relative_eq!(preds[1], 3.0, epsilon = 1e-10);
        assert_relative_eq!(preds[2], 4.0, epsilon = 1e-10);
    }

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

        let model = IsotonicRegression::<f64>::new().with_out_of_bounds(OutOfBounds::Clip);
        let fitted = model.fit(&x, &y).unwrap();

        let x_oob = Array2::from_shape_vec((2, 1), vec![0.0, 4.0]).unwrap();
        let preds = fitted.predict(&x_oob).unwrap();

        // Should clip to the boundary values.
        assert_relative_eq!(preds[0], 1.0, epsilon = 1e-10);
        assert_relative_eq!(preds[1], 3.0, epsilon = 1e-10);
    }

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

        let model = IsotonicRegression::<f64>::new().with_out_of_bounds(OutOfBounds::Nan);
        let fitted = model.fit(&x, &y).unwrap();

        let x_oob = Array2::from_shape_vec((2, 1), vec![0.0, 4.0]).unwrap();
        let preds = fitted.predict(&x_oob).unwrap();

        assert!(preds[0].is_nan());
        assert!(preds[1].is_nan());
    }

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

        let model = IsotonicRegression::<f64>::new().with_out_of_bounds(OutOfBounds::Raise);
        let fitted = model.fit(&x, &y).unwrap();

        let x_below = Array2::from_shape_vec((1, 1), vec![0.0]).unwrap();
        assert!(fitted.predict(&x_below).is_err());

        let x_above = Array2::from_shape_vec((1, 1), vec![4.0]).unwrap();
        assert!(fitted.predict(&x_above).is_err());
    }

    #[test]
    fn test_shape_mismatch_features() {
        let x = Array2::from_shape_vec((3, 2), vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
        let y = array![1.0, 2.0, 3.0];

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

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

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

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

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

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

        let model = IsotonicRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();
        let preds = fitted.predict(&x).unwrap();

        for i in 0..4 {
            assert_relative_eq!(preds[i], 3.0, epsilon = 1e-10);
        }
    }

    #[test]
    fn test_make_unique_weighted_collapse() {
        // Exercises the internal weighted `make_unique` + weighted PAVA that
        // back `_make_unique` (REQ-8) and enable `sample_weight` (REQ-7, #568).
        //
        // Oracle (scikit-learn 1.5.2, sklearn/isotonic.py:317-319 via
        // _isotonic.pyx `_make_unique`):
        //   python3 -c "import numpy as np; from sklearn.isotonic import \
        //   IsotonicRegression; \
        //   m=IsotonicRegression(out_of_bounds='clip').fit( \
        //     np.array([1.,1.,2.,3.]).reshape(-1,1), np.array([1.,3.,2.,4.]), \
        //     sample_weight=np.array([3.,1.,1.,1.])); \
        //   print(m.X_thresholds_.tolist(), m.y_thresholds_.tolist())"
        //   # -> [1.0, 2.0, 3.0] [1.5, 2.0, 4.0]
        //
        // The X=1 run collapses to the weighted mean (3*1 + 1*3)/4 = 1.5, the
        // run weight is 3+1 = 4, and the already-monotone [1.5, 2, 4] is
        // unchanged by the pool.
        let xs = [1.0_f64, 1.0, 2.0, 3.0];
        let ys = [1.0_f64, 3.0, 2.0, 4.0];
        let ws = [3.0_f64, 1.0, 1.0, 1.0];

        let (ux, uy, uw) = make_unique(&xs, &ys, &ws);
        assert_eq!(ux, vec![1.0, 2.0, 3.0]);
        assert_relative_eq!(uy[0], 1.5, epsilon = 1e-12);
        assert_relative_eq!(uy[1], 2.0, epsilon = 1e-12);
        assert_relative_eq!(uy[2], 4.0, epsilon = 1e-12);
        assert_eq!(uw, vec![4.0, 1.0, 1.0]);

        let (rx, ry) = pav_increasing_unique_weighted(&ux, &uy, &uw);
        assert_eq!(rx, vec![1.0, 2.0, 3.0]);
        assert_relative_eq!(ry[0], 1.5, epsilon = 1e-12);
        assert_relative_eq!(ry[1], 2.0, epsilon = 1e-12);
        assert_relative_eq!(ry[2], 4.0, epsilon = 1e-12);
    }

    #[test]
    fn test_unsorted_x() {
        // PAV should handle unsorted x by sorting internally.
        let x = Array2::from_shape_vec((4, 1), vec![3.0, 1.0, 4.0, 2.0]).unwrap();
        let y = array![3.0, 1.0, 4.0, 2.0];

        let model = IsotonicRegression::<f64>::new();
        let fitted = model.fit(&x, &y).unwrap();

        // Predict at sorted x values.
        let x_sorted = Array2::from_shape_vec((4, 1), vec![1.0, 2.0, 3.0, 4.0]).unwrap();
        let preds = fitted.predict(&x_sorted).unwrap();

        for i in 1..preds.len() {
            assert!(preds[i] >= preds[i - 1] - 1e-10);
        }
    }
}