limma/

lowess.rs

//! Cleveland's LOWESS scatterplot smoother: a faithful port of R's `lowess()`
//! / `clowess` (`src/library/stats/src/lowess.c`) and the `weights = NULL`
//! classic path of limma's `loessFit` (`loessFit.R`), together with the
//! prior-weight-aware [`weighted_lowess`] (`weightedLowess.R` +
//! `src/weighted_lowess.c`) used by the weighted path of `loessFit`.
//!
//! Used by the robust *trended* variance moderation in `eBayes`
//! (`fitFDistRobustly` with a covariate), which detrends the log-variances
//! against the covariate with `loessFit(z, covariate, span = 0.4)`. Unlike the
//! spline trend, a smoother has no basis-independence escape hatch: its fitted
//! values *are* the trend, so we reproduce R's `clowess` arithmetic exactly.

/// One local weighted regression of `clowess`: fit at the abscissa `xs` using
/// the points in the 1-based inclusive window `[nleft, nright]` of the sorted
/// data `(x, y)`. Tricube distance weights (multiplied by the robustness
/// weights `rw` when `userw`), then a linear fit — or a locally constant one
/// when the points are too close together to define a slope. Returns
/// `(ys, ok)`, with `ok = false` when every weight is zero. `w` is scratch of
/// length `n`. Mirrors the C `lowest`.
#[allow(clippy::too_many_arguments)]
fn lowest(
    x: &[f64],
    y: &[f64],
    n: usize,
    xs: f64,
    nleft: usize,
    nright: usize,
    w: &mut [f64],
    userw: bool,
    rw: &[f64],
) -> (f64, bool) {
    let range = x[n - 1] - x[0];
    let h = (xs - x[nleft - 1]).max(x[nright - 1] - xs);
    let h9 = 0.999 * h;
    let h1 = 0.001 * h;

    // Sum of weights, picking up all ties on the right.
    let mut a = 0.0;
    let mut j = nleft; // 1-based
    while j <= n {
        w[j - 1] = 0.0;
        let r = (x[j - 1] - xs).abs();
        if r <= h9 {
            if r <= h1 {
                w[j - 1] = 1.0;
            } else {
                let t = 1.0 - (r / h).powi(3);
                w[j - 1] = t * t * t; // tricube
            }
            if userw {
                w[j - 1] *= rw[j - 1];
            }
            a += w[j - 1];
        } else if x[j - 1] > xs {
            break;
        }
        j += 1;
    }
    let nrt = j - 1; // 1-based rightmost point used (>= nright on ties)
    if a <= 0.0 {
        return (0.0, false);
    }

    // Normalise weights to sum to one.
    for ww in w.iter_mut().take(nrt).skip(nleft - 1) {
        *ww /= a;
    }
    if h > 0.0 {
        // Weighted centre of the x values.
        let mut center = 0.0;
        for jj in nleft..=nrt {
            center += w[jj - 1] * x[jj - 1];
        }
        let mut b = xs - center;
        let mut c = 0.0;
        for jj in nleft..=nrt {
            let d = x[jj - 1] - center;
            c += w[jj - 1] * d * d;
        }
        if c.sqrt() > 0.001 * range {
            // Points spread out enough to use a slope.
            b /= c;
            for jj in nleft..=nrt {
                w[jj - 1] *= b * (x[jj - 1] - center) + 1.0;
            }
        }
    }
    let mut ys = 0.0;
    for jj in nleft..=nrt {
        ys += w[jj - 1] * y[jj - 1];
    }
    (ys, true)
}

/// Cleveland's LOWESS: smooth `y` against ascending-sorted `x`. `f` is the
/// span (fraction of points in each local fit), `nsteps` the number of
/// robustness iterations, `delta` the interpolation skip distance. Returns the
/// fitted values in the (sorted) input order. Faithful port of `clowess`.
fn clowess(x: &[f64], y: &[f64], f: f64, nsteps: usize, delta: f64) -> Vec<f64> {
    let n = x.len();
    let mut ys = vec![0.0_f64; n];
    if n < 2 {
        if n == 1 {
            ys[0] = y[0];
        }
        return ys;
    }

    // Window size: at least two, at most n points.
    let ns = (((f * n as f64) + 1e-7) as usize).clamp(2, n);

    let mut rw = vec![0.0_f64; n]; // robustness weights
    let mut res = vec![0.0_f64; n]; // residuals
    let mut w = vec![0.0_f64; n]; // scratch for `lowest`

    let mut iter = 1usize;
    while iter <= nsteps + 1 {
        let mut nleft = 1usize; // 1-based window endpoints
        let mut nright = ns;
        let mut last = 0usize; // 1-based index of previous estimated point (0 = none)
        let mut i = 1usize; // 1-based current point

        loop {
            if nright < n {
                // Move the window right while that brings the far edge closer.
                let d1 = x[i - 1] - x[nleft - 1];
                let d2 = x[nright] - x[i - 1]; // x[nright+1] (1-based)
                if d1 > d2 {
                    nleft += 1;
                    nright += 1;
                    continue;
                }
            }

            let (yi, ok) = lowest(x, y, n, x[i - 1], nleft, nright, &mut w, iter > 1, &rw);
            ys[i - 1] = if ok { yi } else { y[i - 1] };

            // Interpolate any points skipped by the delta jump.
            if last < i - 1 {
                let denom = x[i - 1] - x[last - 1];
                for j in (last + 1)..i {
                    let alpha = (x[j - 1] - x[last - 1]) / denom;
                    ys[j - 1] = alpha * ys[i - 1] + (1.0 - alpha) * ys[last - 1];
                }
            }
            last = i;

            // Skip points within delta of the last estimated abscissa; copy the
            // fit across exact ties.
            let cut = x[last - 1] + delta;
            i = last + 1;
            while i <= n {
                if x[i - 1] > cut {
                    break;
                }
                if x[i - 1] == x[last - 1] {
                    ys[i - 1] = ys[last - 1];
                    last = i;
                }
                i += 1;
            }
            i = (last + 1).max(i - 1);
            if last >= n {
                break;
            }
        }

        for k in 0..n {
            res[k] = y[k] - ys[k];
        }
        let sc: f64 = res.iter().map(|v| v.abs()).sum::<f64>() / n as f64;

        // No reweighting after the final pass.
        if iter > nsteps {
            break;
        }
        for k in 0..n {
            rw[k] = res[k].abs();
        }
        // cmad = 6 * median(|res|). `rw` currently holds |res| and is fully
        // overwritten with the bisquare weights below, so we can quickselect in
        // place — no clone, no full O(n log n) sort. `select_nth_unstable_by`
        // puts the m1-th order statistic at index m1 with every smaller value in
        // `rw[..m1]`; for even n the (m1-1)-th statistic is the max of that lower
        // partition, giving the same two central values, summed in the same
        // order, that the full sort produced.
        let m1 = n / 2;
        let cmad = {
            rw.select_nth_unstable_by(m1, |a, b| a.partial_cmp(b).unwrap());
            if n % 2 == 0 {
                let lo = rw[..m1].iter().copied().fold(f64::NEG_INFINITY, f64::max);
                3.0 * (rw[m1] + lo)
            } else {
                6.0 * rw[m1]
            }
        };
        if cmad < 1e-7 * sc {
            break; // residuals effectively zero
        }
        let c9 = 0.999 * cmad;
        let c1 = 0.001 * cmad;
        for k in 0..n {
            let r = res[k].abs();
            rw[k] = if r <= c1 {
                1.0
            } else if r <= c9 {
                let t = 1.0 - (r / cmad).powi(2);
                t * t // bisquare
            } else {
                0.0
            };
        }
        iter += 1;
    }
    ys
}

/// The `weights = NULL` classic-lowess path of `loessFit`, used directly by
/// `fitFDistRobustly` (covariate detrend) and `normalizeCyclicLoess`, and by
/// the public [`loess_fit`] wrapper when weights are absent or all-equal.
/// Returns `(fitted, residuals)` aligned to the original order of `(y, x)`;
/// non-finite observations are left as `NaN`, exactly as `loessFit` does.
pub(crate) fn loess_fit_unweighted(
    y: &[f64],
    x: &[f64],
    span: f64,
    iterations: usize,
) -> (Vec<f64>, Vec<f64>) {
    let n = y.len();
    let mut fitted = vec![f64::NAN; n];
    let mut residuals = vec![f64::NAN; n];

    let obs: Vec<usize> = (0..n)
        .filter(|&k| y[k].is_finite() && x[k].is_finite())
        .collect();
    let nobs = obs.len();
    if nobs == 0 {
        return (fitted, residuals);
    }
    if span < 1.0 / nobs as f64 {
        for &k in &obs {
            fitted[k] = y[k];
            residuals[k] = 0.0;
        }
        return (fitted, residuals);
    }

    // Sort the observations by x (stable, matching R's order()), smooth, then
    // scatter the fit back to the original positions.
    let mut order = obs.clone();
    order.sort_by(|&a, &b| x[a].partial_cmp(&x[b]).unwrap());
    let xs: Vec<f64> = order.iter().map(|&k| x[k]).collect();
    let ys: Vec<f64> = order.iter().map(|&k| y[k]).collect();
    let delta = 0.01 * (xs[nobs - 1] - xs[0]);

    let fit = clowess(&xs, &ys, span, iterations.saturating_sub(1), delta);
    for (j, &k) in order.iter().enumerate() {
        fitted[k] = fit[j];
        residuals[k] = y[k] - fit[j];
    }
    (fitted, residuals)
}

/// Seed points for `weightedLowess`: indices spaced more than `delta` apart in
/// the ascending covariate `x`, always including the first and last point.
/// Mirrors the C `find_seeds`.
fn find_seeds(x: &[f64], delta: f64) -> Vec<usize> {
    let npts = x.len();
    let mut indices = vec![0usize];
    let mut last_pt = 0usize;
    for pt in 1..npts.saturating_sub(1) {
        if x[pt] - x[last_pt] > delta {
            indices.push(pt);
            last_pt = pt;
        }
    }
    indices.push(npts - 1);
    indices
}

/// For each seed, expand a window left/right accumulating prior weights until
/// the span weight is reached, always extending by the side that adds the
/// smaller distance, then out to any ties. Returns the per-seed inclusive
/// window `(start, end)` (0-based) and the maximum in-window distance. Mirrors
/// the C `find_limits`.
fn find_limits(
    indices: &[usize],
    x: &[f64],
    w: &[f64],
    spanweight: f64,
) -> (Vec<usize>, Vec<usize>, Vec<f64>) {
    let npts = x.len();
    let num = indices.len();
    let mut start = vec![0usize; num];
    let mut end = vec![0usize; num];
    let mut dist = vec![0.0f64; num];

    for (curx, &curpt) in indices.iter().enumerate() {
        let mut left = curpt;
        let mut right = curpt;
        let mut curw = w[curpt];
        let mut at_end = curpt == npts - 1;
        let mut at_start = curpt == 0;
        let mut mdist = 0.0f64;

        while curw < spanweight && (!at_end || !at_start) {
            if at_end {
                left -= 1;
                curw += w[left];
                if left == 0 {
                    at_start = true;
                }
                let ldist = x[curpt] - x[left];
                if mdist < ldist {
                    mdist = ldist;
                }
            } else if at_start {
                right += 1;
                curw += w[right];
                if right == npts - 1 {
                    at_end = true;
                }
                let rdist = x[right] - x[curpt];
                if mdist < rdist {
                    mdist = rdist;
                }
            } else {
                let ldist = x[curpt] - x[left - 1];
                let rdist = x[right + 1] - x[curpt];
                if ldist < rdist {
                    left -= 1;
                    curw += w[left];
                    if left == 0 {
                        at_start = true;
                    }
                    if mdist < ldist {
                        mdist = ldist;
                    }
                } else {
                    right += 1;
                    curw += w[right];
                    if right == npts - 1 {
                        at_end = true;
                    }
                    if mdist < rdist {
                        mdist = rdist;
                    }
                }
            }
        }

        // Extend symmetrically across exact ties.
        while left > 0 && x[left] == x[left - 1] {
            left -= 1;
        }
        while right < npts - 1 && x[right] == x[right + 1] {
            right += 1;
        }

        start[curx] = left;
        end[curx] = right;
        dist[curx] = mdist;
    }
    (start, end, dist)
}

/// Local weighted linear regression at `x[curpt]` over the inclusive window
/// `[left, right]`, combining tricube distance weights with the prior weights
/// `w` and robustness weights `rw`. Falls back to a weighted mean when the
/// window distance or covariate variance is numerically negligible. `work` is
/// scratch of length `npts`. Mirrors the C `lowess_fit`.
#[allow(clippy::too_many_arguments)]
fn weighted_local_fit(
    x: &[f64],
    y: &[f64],
    w: &[f64],
    rw: &[f64],
    curpt: usize,
    left: usize,
    right: usize,
    dist: f64,
    work: &mut [f64],
) -> f64 {
    if dist < 1e-7 {
        let mut ymean = 0.0;
        let mut allweight = 0.0;
        for pt in left..=right {
            work[pt] = w[pt] * rw[pt];
            ymean += y[pt] * work[pt];
            allweight += work[pt];
        }
        return ymean / allweight;
    }
    let mut xmean = 0.0;
    let mut ymean = 0.0;
    let mut allweight = 0.0;
    for pt in left..=right {
        let tri = (1.0 - ((x[curpt] - x[pt]).abs() / dist).powi(3)).powi(3);
        work[pt] = tri * w[pt] * rw[pt];
        xmean += work[pt] * x[pt];
        ymean += work[pt] * y[pt];
        allweight += work[pt];
    }
    xmean /= allweight;
    ymean /= allweight;

    let mut var = 0.0;
    let mut covar = 0.0;
    for pt in left..=right {
        let temp = x[pt] - xmean;
        var += temp * temp * work[pt];
        covar += temp * (y[pt] - ymean) * work[pt];
    }
    if var < 1e-7 {
        return ymean;
    }
    let slope = covar / var;
    let intercept = ymean - slope * xmean;
    slope * x[curpt] + intercept
}

/// The robustness-iteration loop of the C `weighted_lowess`, operating on data
/// already sorted by ascending `x`. Returns the fitted values and the final
/// Tukey-biweight robustness weights, both in sorted order. The weights are
/// recomputed from the final fit's residuals but (as in the C) not themselves
/// applied to any fit.
fn weighted_lowess_core(
    x: &[f64],
    y: &[f64],
    w: &[f64],
    span: f64,
    iterations: usize,
    delta: f64,
) -> (Vec<f64>, Vec<f64>) {
    let n = x.len();
    let totalweight: f64 = w.iter().sum();
    let spanweight = totalweight * span;
    let subrange = (x[n - 1] - x[0]) / n as f64;

    let seeds = find_seeds(x, delta);
    let (frame_start, frame_end, max_dist) = find_limits(&seeds, x, w, spanweight);

    let mut fit = vec![0.0f64; n];
    let mut rob = vec![1.0f64; n];
    let mut work = vec![0.0f64; n];
    // MAD scratch reused across robustness iterations (refilled each pass).
    let mut absres = vec![0.0f64; n];
    let mut ror = vec![0usize; n];
    let mut sorted = vec![0.0f64; n];

    for _ in 0..iterations.max(1) {
        // Fit at the seeds; linearly interpolate the curve between them.
        fit[0] = weighted_local_fit(
            x,
            y,
            w,
            &rob,
            0,
            frame_start[0],
            frame_end[0],
            max_dist[0],
            &mut work,
        );
        let mut last_pt = 0usize;
        for cur_seed in 1..seeds.len() {
            let pt = seeds[cur_seed];
            fit[pt] = weighted_local_fit(
                x,
                y,
                w,
                &rob,
                pt,
                frame_start[cur_seed],
                frame_end[cur_seed],
                max_dist[cur_seed],
                &mut work,
            );
            if pt - last_pt > 1 {
                let current = x[pt] - x[last_pt];
                if current > 1e-7 * subrange {
                    let slope = (fit[pt] - fit[last_pt]) / current;
                    let intercept = fit[pt] - slope * x[pt];
                    for subpt in (last_pt + 1)..pt {
                        fit[subpt] = slope * x[subpt] + intercept;
                    }
                } else {
                    let endave = 0.5 * (fit[pt] + fit[last_pt]);
                    for f in fit.iter_mut().take(pt).skip(last_pt + 1) {
                        *f = endave;
                    }
                }
            }
            last_pt = pt;
        }

        // Weighted MAD of the absolute residuals (reusing hoisted scratch).
        let mut resid_scale = 0.0;
        for pt in 0..n {
            absres[pt] = (y[pt] - fit[pt]).abs();
            resid_scale += absres[pt];
        }
        resid_scale /= n as f64;

        for (i, r) in ror.iter_mut().enumerate() {
            *r = i;
        }
        ror.sort_by(|&a, &b| absres[a].partial_cmp(&absres[b]).unwrap());
        for pt in 0..n {
            sorted[pt] = absres[ror[pt]];
        }

        let halfweight = totalweight / 2.0;
        let mut current = 0.0;
        let mut cmad = 0.0;
        for pt in 0..n {
            current += w[ror[pt]];
            if current == halfweight {
                let next = if pt + 1 < n {
                    sorted[pt + 1]
                } else {
                    sorted[pt]
                };
                cmad = 3.0 * (sorted[pt] + next);
                break;
            } else if current > halfweight {
                cmad = 6.0 * sorted[pt];
                break;
            }
        }

        if cmad <= 1e-7 * resid_scale {
            break;
        }

        for pt in 0..n {
            rob[ror[pt]] = if sorted[pt] < cmad {
                let t = 1.0 - (sorted[pt] / cmad).powi(2);
                t * t
            } else {
                0.0
            };
        }
    }
    (fit, rob)
}

/// Fitted curve, residuals and robustness weights returned by
/// [`weighted_lowess`], all aligned to the original input order.
pub struct WeightedLowess {
    pub fitted: Vec<f64>,
    pub residuals: Vec<f64>,
    /// Final Tukey-biweight robustness weights (the `weights` component of
    /// limma's `weightedLowess`, *not* the input prior weights).
    pub weights: Vec<f64>,
    pub delta: f64,
}

/// Port of limma's `weightedLowess` (`weightedLowess.R` +
/// `src/weighted_lowess.c`): a prior-weight-aware lowess smoother of degree 1
/// with Tukey-biweight robustness iterations. `weights` defaults to all-ones
/// when `None`; `delta` (the interpolation skip distance) is chosen from `npts`
/// evenly spaced clusters when `None`, exactly as the R wrapper does. The
/// result matches the default `output.style = "loess"`: `fitted`, `residuals`
/// and robustness `weights` in the original (unsorted) order.
pub fn weighted_lowess(
    x: &[f64],
    y: &[f64],
    weights: Option<&[f64]>,
    span: f64,
    iterations: usize,
    npts: usize,
    delta: Option<f64>,
) -> WeightedLowess {
    let n = x.len();
    assert_eq!(y.len(), n, "x and y should have same length");
    let w_in: Vec<f64> = match weights {
        Some(wv) => {
            assert_eq!(wv.len(), n, "weights should have same length as x and y");
            wv.to_vec()
        }
        None => vec![1.0; n],
    };

    // Order by ascending x (stable, matching R's order()).
    let mut o: Vec<usize> = (0..n).collect();
    o.sort_by(|&a, &b| x[a].partial_cmp(&x[b]).unwrap());
    let xs: Vec<f64> = o.iter().map(|&k| x[k]).collect();
    let ys: Vec<f64> = o.iter().map(|&k| y[k]).collect();
    let ws: Vec<f64> = o.iter().map(|&k| w_in[k]).collect();

    // Choose delta from the cluster structure when not supplied.
    let delta = match delta {
        Some(d) => d,
        None if npts >= n => 0.0,
        None => {
            let mut dx: Vec<f64> = (0..n - 1).map(|i| xs[i + 1] - xs[i]).collect();
            dx.sort_by(|a, b| a.partial_cmp(b).unwrap());
            let mut cum = 0.0;
            let cumrange: Vec<f64> = dx
                .iter()
                .map(|&d| {
                    cum += d;
                    cum
                })
                .collect();
            let mut best = f64::INFINITY;
            for k in 0..npts {
                let cand = cumrange[n - 2 - k] / (npts - k) as f64;
                if cand < best {
                    best = cand;
                }
            }
            best
        }
    };

    let (fit_sorted, rob_sorted) = weighted_lowess_core(&xs, &ys, &ws, span, iterations, delta);

    // Scatter the sorted fit back to the original order.
    let mut fitted = vec![0.0f64; n];
    let mut robweights = vec![0.0f64; n];
    for (j, &k) in o.iter().enumerate() {
        fitted[k] = fit_sorted[j];
        robweights[k] = rob_sorted[j];
    }
    let residuals: Vec<f64> = (0..n).map(|i| y[i] - fitted[i]).collect();

    WeightedLowess {
        fitted,
        residuals,
        weights: robweights,
        delta,
    }
}

/// Fitted values and residuals returned by [`loess_fit`], aligned to the input
/// order with `NaN` at non-finite observations.
pub struct LoessFit {
    pub fitted: Vec<f64>,
    pub residuals: Vec<f64>,
}

/// Port of limma's `loessFit` for the default `method = "weightedLowess"`: a
/// fast univariate lowess fit allowing prior weights. Non-finite `(y, x)` pairs
/// are dropped and returned as `NaN`. With no weights — or weights that are all
/// equal when `equal_weights_as_null` — it uses the classic
/// [`loess_fit_unweighted`] path; otherwise weights are clamped to
/// `[min_weight, max_weight]` (NaN -> 0 first) and the fit uses
/// [`weighted_lowess`], dropping to a weighted straight-line fit when there are
/// too few observations to support a smoother. The `method = "locfit"` and
/// `method = "loess"` variants (which call the R packages of those names) are
/// out of scope for this pure-Rust port.
#[allow(clippy::too_many_arguments)]
pub fn loess_fit(
    y: &[f64],
    x: &[f64],
    weights: Option<&[f64]>,
    span: f64,
    iterations: usize,
    min_weight: f64,
    max_weight: f64,
    equal_weights_as_null: bool,
) -> LoessFit {
    let n = y.len();
    assert_eq!(x.len(), n, "y and x have different lengths");
    let mut fitted = vec![f64::NAN; n];
    let mut residuals = vec![f64::NAN; n];

    let obs: Vec<usize> = (0..n)
        .filter(|&k| y[k].is_finite() && x[k].is_finite())
        .collect();
    let nobs = obs.len();
    if nobs == 0 {
        return LoessFit { fitted, residuals };
    }

    // Spans below 1/nobs short-circuit to interpolation (fitted = y).
    if span < 1.0 / nobs as f64 {
        for &k in &obs {
            fitted[k] = y[k];
            residuals[k] = 0.0;
        }
        return LoessFit { fitted, residuals };
    }

    let min_weight = min_weight.max(0.0);
    let xobs: Vec<f64> = obs.iter().map(|&k| x[k]).collect();
    let yobs: Vec<f64> = obs.iter().map(|&k| y[k]).collect();

    // Resolve the prior weights over the observed points; `None` means the
    // classic unweighted path (weights absent, or all equal and treated as
    // null).
    let wobs: Option<Vec<f64>> = match weights {
        None => None,
        Some(wv) => {
            assert_eq!(wv.len(), n, "y and weights have different lengths");
            let wo: Vec<f64> = obs
                .iter()
                .map(|&k| {
                    let w = if wv[k].is_nan() { 0.0 } else { wv[k] };
                    w.max(min_weight).min(max_weight)
                })
                .collect();
            if equal_weights_as_null {
                let lo = wo.iter().copied().fold(f64::INFINITY, f64::min);
                let hi = wo.iter().copied().fold(f64::NEG_INFINITY, f64::max);
                if hi - lo < 1e-15 {
                    None
                } else {
                    Some(wo)
                }
            } else {
                Some(wo)
            }
        }
    };

    let Some(wobs) = wobs else {
        let (f, _) = loess_fit_unweighted(&yobs, &xobs, span, iterations);
        for (j, &k) in obs.iter().enumerate() {
            fitted[k] = f[j];
            residuals[k] = y[k] - f[j];
        }
        return LoessFit { fitted, residuals };
    };

    // Number of observations with positive weight (always positive).
    let nwobs = if min_weight > 0.0 {
        nobs
    } else {
        wobs.iter().filter(|&&w| w > 0.0).count()
    };

    // Too few observations to estimate a lowess curve: a single point, or a
    // weighted straight line.
    if (nwobs as f64) < 4.0 + 1.0 / span {
        if nwobs == 1 {
            let kpos = wobs.iter().position(|&w| w > 0.0).unwrap();
            let val = yobs[kpos];
            for (j, &k) in obs.iter().enumerate() {
                fitted[k] = val;
                residuals[k] = yobs[j] - val;
            }
        } else {
            let (slope, intercept) = weighted_line_fit(&xobs, &yobs, &wobs);
            for (j, &k) in obs.iter().enumerate() {
                let f = intercept + slope * xobs[j];
                fitted[k] = f;
                residuals[k] = yobs[j] - f;
            }
        }
        return LoessFit { fitted, residuals };
    }

    let wl = weighted_lowess(&xobs, &yobs, Some(&wobs), span, iterations, 200, None);
    for (j, &k) in obs.iter().enumerate() {
        fitted[k] = wl.fitted[j];
        residuals[k] = y[k] - wl.fitted[j];
    }
    LoessFit { fitted, residuals }
}

/// Weighted least squares of `y` on `[1, x]` (the `lm.wfit(cbind(1, x), y, w)`
/// fallback inside `loessFit`), returning `(slope, intercept)`. Closed form for
/// simple linear regression; zero-weight rows drop out of the normal equations.
fn weighted_line_fit(x: &[f64], y: &[f64], w: &[f64]) -> (f64, f64) {
    let mut sw = 0.0;
    let mut swx = 0.0;
    let mut swy = 0.0;
    let mut swxx = 0.0;
    let mut swxy = 0.0;
    for ((&xi, &yi), &wi) in x.iter().zip(y.iter()).zip(w.iter()) {
        sw += wi;
        swx += wi * xi;
        swy += wi * yi;
        swxx += wi * xi * xi;
        swxy += wi * xi * yi;
    }
    let denom = sw * swxx - swx * swx;
    let slope = (sw * swxy - swx * swy) / denom;
    let intercept = (swy - slope * swx) / sw;
    (slope, intercept)
}

/// R's `lowess(x, y, f, iter)`: order the points by `x`, smooth with `clowess`
/// (`delta = 0.01 * range(x)`, as in R), and return the sorted abscissae
/// together with the fitted ordinates — i.e. the `$x` and `$y` of `lowess`.
/// Assumes finite inputs (the variance-trend callers guarantee this).
fn lowess_xy(x: &[f64], y: &[f64], f: f64, iter: usize) -> (Vec<f64>, Vec<f64>) {
    let n = x.len();
    if n == 0 {
        return (Vec::new(), Vec::new());
    }
    let mut order: Vec<usize> = (0..n).collect();
    order.sort_by(|&a, &b| x[a].partial_cmp(&x[b]).unwrap());
    let xs: Vec<f64> = order.iter().map(|&k| x[k]).collect();
    let ys: Vec<f64> = order.iter().map(|&k| y[k]).collect();
    let delta = 0.01 * (xs[n - 1] - xs[0]);
    let fit = clowess(&xs, &ys, f, iter, delta);
    (xs, fit)
}

/// Collapse exact-duplicate abscissae by averaging their ordinates, matching
/// `approxfun(..., ties = list("ordered", mean))`. Input must be ascending.
fn collapse_ties(xs: &[f64], ys: &[f64]) -> (Vec<f64>, Vec<f64>) {
    let mut gx = Vec::new();
    let mut gy = Vec::new();
    let mut i = 0usize;
    while i < xs.len() {
        let mut j = i;
        let mut sum = 0.0;
        while j < xs.len() && xs[j] == xs[i] {
            sum += ys[j];
            j += 1;
        }
        gx.push(xs[i]);
        gy.push(sum / (j - i) as f64);
        i = j;
    }
    (gx, gy)
}

/// Linear interpolation over an ascending grid with `rule = 2` (values outside
/// the grid clamp to the nearest endpoint), matching R's `approx`/`approxfun`.
pub(crate) fn approx_rule2(gx: &[f64], gy: &[f64], xout: f64) -> f64 {
    let n = gx.len();
    if n == 1 {
        return gy[0];
    }
    if xout <= gx[0] {
        return gy[0];
    }
    if xout >= gx[n - 1] {
        return gy[n - 1];
    }
    let i = gx.partition_point(|&v| v <= xout); // first index with gx > xout
    let (x0, x1) = (gx[i - 1], gx[i]);
    let (y0, y1) = (gy[i - 1], gy[i]);
    y0 + (y1 - y0) * (xout - x0) / (x1 - x0)
}

/// A lowess variance trend together with its `approxfun` interpolation rule,
/// as used by `voom`/`vooma`: fit `lowess(x, y, f = span)` then build a
/// `rule = 2`, tie-averaged linear interpolator over the fitted curve.
pub(crate) struct LowessInterpolator {
    gx: Vec<f64>,
    gy: Vec<f64>,
}

impl LowessInterpolator {
    pub(crate) fn fit(x: &[f64], y: &[f64], span: f64, iter: usize) -> Self {
        let (xs, ys) = lowess_xy(x, y, span, iter);
        let (gx, gy) = collapse_ties(&xs, &ys);
        Self { gx, gy }
    }

    pub(crate) fn eval(&self, xout: f64) -> f64 {
        approx_rule2(&self.gx, &self.gy, xout)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn rel(a: f64, b: f64) -> f64 {
        ((a - b) / b).abs()
    }

    /// Reference values from R 4.6.0:
    /// `x <- 1:20; y <- sin(x/3) + (x %% 4)/10`
    /// `lowess(x, y, f = 0.4, iter = 3L)$y`.
    /// The classic-lowess path of `loessFit(y, x, span = 0.4)` reduces to this.
    #[test]
    fn lowess_matches_r_reference() {
        let x: Vec<f64> = (1..=20).map(|i| i as f64).collect();
        let y: Vec<f64> = (1..=20)
            .map(|i| {
                let xi = i as f64;
                (xi / 3.0).sin() + ((i % 4) as f64) / 10.0
            })
            .collect();
        // loessFit(y, x, span=0.4) classic path: lowess(x, y, f=0.4, iter=3).
        let (fitted, residuals) = loess_fit_unweighted(&y, &x, 0.4, 4);

        // R lowess(x, y, f=0.4, iter=3L)$y (x already ascending, so $y aligns
        // with the input order).
        let r_fitted = [
            0.596_248_207_878_296,
            0.733_833_572_576_796,
            0.865_358_555_762_325,
            0.994_091_762_946_332,
            1.012_360_782_819_75,
            0.898_791_455_642_421,
            0.731_283_317_640_518,
            0.529_277_971_327_881,
            0.269_444_204_960_129,
            -0.012_657_116_933_792_8,
            -0.266_008_842_117_912,
            -0.457_558_931_089_973,
            -0.613_194_376_883_428,
            -0.715_648_998_499_32,
            -0.676_645_471_029_55,
            -0.521_458_054_754_246,
            -0.319_116_552_594_766,
            -0.077_757_590_217_849,
            0.171_815_169_826_499,
            0.431_111_024_786_567,
        ];
        for (a, b) in fitted.iter().zip(r_fitted.iter()) {
            assert!(rel(*a, *b) < 1e-9, "fitted {a} vs {b}");
        }
        // Residuals are y - fitted by construction.
        for k in 0..x.len() {
            assert!((residuals[k] - (y[k] - fitted[k])).abs() < 1e-15);
        }
    }

    /// Spans below 1/n short-circuit to interpolation (fitted = y, resid = 0).
    #[test]
    fn lowess_tiny_span_is_identity() {
        let x: Vec<f64> = (1..=10).map(|i| i as f64).collect();
        let y: Vec<f64> = (1..=10).map(|i| (i * i) as f64).collect();
        let (fitted, residuals) = loess_fit_unweighted(&y, &x, 0.05, 4);
        for k in 0..x.len() {
            assert_eq!(fitted[k], y[k]);
            assert_eq!(residuals[k], 0.0);
        }
    }

    /// Relative match, tolerating exact-zero references (zeroed robustness
    /// weights at the planted outliers).
    fn close(a: f64, b: f64) {
        if b == 0.0 {
            assert!(a.abs() < 1e-12, "{a} vs 0");
        } else {
            assert!(((a - b) / b).abs() < 1e-9, "{a} vs {b}");
        }
    }

    /// Shared n=24 weighted scenario (two planted outliers) for cases A and B.
    fn case_ab_data() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
        let x = vec![
            0.45749338542921825,
            0.36490716953341845,
            2.1216858507761938,
            2.231_517_906_338_681,
            2.4560694953036157,
            3.003_471_605_493_78,
            3.2937879157555763,
            3.619468372607308,
            4.1084805104211775,
            4.5866320499433755,
            5.444_886_515_827_218,
            5.775_470_171_559_732,
            5.940_419_082_916_579,
            6.149_600_128_386_905,
            6.600_257_345_624_926,
            6.604648996606679,
            7.030_477_401_104_289,
            7.093_654_185_978_542,
            7.141587042811965,
            7.333_683_544_946_952,
            8.254_754_236_391_259,
            8.694_301_976_611_142,
            9.375_836_948_985_027,
            9.532_155_904_709_773,
        ];
        let y = vec![
            0.612_545_426_481_486_9,
            0.650_401_760_370_632,
            1.1542440348430094,
            1.5486987991164907,
            4.466_390_679_222_614,
            1.1903665393098684,
            0.37266251854917354,
            0.534_045_647_135_347_6,
            0.18834499131155025,
            0.464_472_256_763_577_8,
            1.0037463888693132,
            1.1052425234922978,
            1.3880139159873457,
            1.4149161956436145,
            2.061_812_038_829_897,
            2.2424560756224436,
            2.9043783394714633,
            0.073_326_112_417_279_27,
            3.0490864839234257,
            2.8575387507063152,
            3.430_262_298_625_616,
            3.501_476_667_740_035,
            3.0964170542301663,
            2.666902407201408,
        ];
        let w = vec![
            0.915_514_786_448_329_7,
            0.45188319422304635,
            0.34237423320300875,
            1.1903987498488278,
            0.669_161_357_125_267_5,
            1.6571537321433425,
            1.1805854137521237,
            1.054087870148942,
            1.3946744401939213,
            0.36591726122424006,
            1.3710096625145525,
            0.863_262_355_374_172_3,
            0.642_421_815_311_536_2,
            0.737_486_680_131_405_7,
            1.2065324763767422,
            1.064640044188127,
            1.0586248501669615,
            1.8839623192790895,
            1.046_319_722_617_045,
            1.4207539540715515,
            1.8820373674388975,
            0.693_199_012_149_125_4,
            1.9052787743508817,
            0.763_943_711_994_215_8,
        ];
        (x, y, w)
    }

    /// Reference: `weightedLowess(x, y, weights=w, span=0.3, iterations=4L,
    /// npts=200)` from limma 3.68.3 (npts>=n so delta=0; full robustness).
    #[test]
    fn weighted_lowess_case_a() {
        let (x, y, w) = case_ab_data();
        let out = weighted_lowess(&x, &y, Some(&w), 0.3, 4, 200, None);
        assert_eq!(out.delta, 0.0);

        let r_fitted = [
            0.704_865_579_518_969_9,
            0.675_018_400_307_334_6,
            1.3600628188780086,
            1.4805324495189862,
            1.3743558087666423,
            1.0285760664659196,
            0.852_100_138_874_763_8,
            0.617_496_271_481_045_2,
            0.513_999_292_241_682_6,
            0.610_689_342_452_828_6,
            0.988_412_214_777_719_5,
            1.2375884034126585,
            1.3951903052897316,
            1.5905658893023542,
            2.152258891028989,
            2.1589538558884804,
            2.6882075718940692,
            2.754_372_887_365_016,
            2.804_782_373_798_898,
            2.972630547835498,
            3.396_904_942_104_693,
            3.2853367132950937,
            3.0552758361559125,
            2.996_975_009_544_398,
        ];
        let r_rweights = [
            0.964_573_093_916_010_6,
            0.9974600786593063,
            0.830_223_645_617_545_5,
            0.980_606_432_031_510_2,
            0.0,
            0.8932277961273466,
            0.26829417855553755,
            0.971_005_534_450_722_9,
            0.604_671_627_514_471_7,
            0.912_342_577_459_221_6,
            0.9990140555563487,
            0.927_887_678_710_836_4,
            0.999_784_013_165_360_9,
            0.874_786_293_731_308_1,
            0.965_983_909_872_752_1,
            0.970_969_934_423_463_4,
            0.813_613_521_418_203_6,
            0.0,
            0.765_342_631_054_660_8,
            0.945_216_049_382_716,
            0.995_338_614_565_773_8,
            0.813_663_916_297_942_5,
            0.992_913_666_264_513_6,
            0.595_259_735_008_007_2,
        ];
        for (a, b) in out.fitted.iter().zip(r_fitted.iter()) {
            close(*a, *b);
        }
        for (a, b) in out.weights.iter().zip(r_rweights.iter()) {
            close(*a, *b);
        }
        for ((r, yy), f) in out.residuals.iter().zip(y.iter()).zip(out.fitted.iter()) {
            assert!((r - (yy - f)).abs() < 1e-15);
        }
    }

    /// Reference: same data, `span=0.5, iterations=1L` — a non-robust base fit
    /// whose returned robustness weights are computed but never applied.
    #[test]
    fn weighted_lowess_case_b() {
        let (x, y, w) = case_ab_data();
        let out = weighted_lowess(&x, &y, Some(&w), 0.5, 1, 200, None);
        assert_eq!(out.delta, 0.0);

        let r_fitted = [
            1.3164091907640776,
            1.3175754346917248,
            1.2097472947654593,
            1.2023153256037815,
            1.1902171010985594,
            1.1769344746702706,
            1.1383757235477594,
            1.004_305_077_879_562,
            0.796_818_060_092_517,
            0.774_791_789_454_579_5,
            1.031_357_331_214_907,
            1.2516144937670477,
            1.358_007_678_833_491,
            1.4738918429020829,
            1.7425837851859414,
            1.7453266670230196,
            2.0349113652577646,
            2.099_264_238_862_616,
            2.148618243528349,
            2.340_243_909_841_483,
            3.010_161_361_558_707,
            3.056683040538843,
            3.158_235_211_076_587,
            3.1633262263048074,
        ];
        let r_rweights = [
            0.879_412_603_471_390_4,
            0.8913033071647305,
            0.999_226_241_529_259_7,
            0.9700854954413779,
            0.0,
            0.999_954_675_492_639_4,
            0.858_129_979_744_656_7,
            0.945_216_049_382_716,
            0.909_151_582_734_925_9,
            0.975_954_372_722_726_4,
            0.999_808_488_820_543_1,
            0.994_624_937_295_328_6,
            0.999_773_821_675_273_3,
            0.999_126_419_320_665,
            0.974_562_958_604_831_9,
            0.9388779867911321,
            0.819_102_172_057_449_2,
            0.2346871608581525,
            0.806_674_164_798_635_2,
            0.933_905_106_727_050_2,
            0.956_155_094_856_453_8,
            0.950_916_050_426_533_5,
            0.999_040_200_929_231_7,
            0.939_048_642_187_021_6,
        ];
        for (a, b) in out.fitted.iter().zip(r_fitted.iter()) {
            close(*a, *b);
        }
        for (a, b) in out.weights.iter().zip(r_rweights.iter()) {
            close(*a, *b);
        }
    }

    /// Reference: n=80, `npts=30` so the wrapper picks delta>0 and the C path
    /// samples seeds and interpolates between them. `span=0.3, iterations=4L`.
    #[test]
    fn weighted_lowess_case_c() {
        let x = [
            0.531_976_991_333_067_4,
            0.683_757_108_636_200_4,
            0.789_968_837_052_583_7,
            0.801_603_901_199_996_5,
            1.5255942661315203,
            1.581659372895956,
            2.577_251_400_798_559,
            2.7313429210335016,
            3.1198115088045597,
            3.2957231905311346,
            3.5041399160400033,
            3.535_446_380_265_057,
            3.562_620_780_430_734,
            3.573_531_210_422_516,
            3.904986004345119,
            3.969_677_845_016_122,
            4.026_793_888_770_044,
            4.0656226221472025,
            4.074_543_830_938_637,
            4.320_158_297_196_031,
            4.373_944_881_372_154,
            4.3781809182837605,
            4.381_559_635_512_531,
            4.541_048_249_229_789,
            4.595673899166286,
            5.085_642_351_768_911,
            5.320_935_468_189_418,
            5.745_922_313_071_787,
            5.820684814825654,
            6.121_302_433_311_939,
            6.595_797_222_107_649,
            7.196_558_876_894_414,
            7.365_691_293_962_3,
            7.962_349_629_960_954,
            8.258_966_878_056_526,
            8.290_081_201_121_211,
            8.838_997_278_362_513,
            8.898_926_260_881_126,
            9.147_843_788_377_94,
            9.250_376_787_967_98,
            9.582185884937644,
            10.625_021_490_268_41,
            10.671597444452345,
            11.63403379265219,
            11.695629926398396,
            11.993930744938552,
            12.457542419433594,
            12.604068955406547,
            12.657269365154207,
            12.804370601661503,
            12.825214397162199,
            12.832566713914275,
            12.845537355169654,
            13.152610226534307,
            13.243088563904166,
            13.403024673461914,
            13.606816763058305,
            13.800439438782632,
            13.916632984764874,
            14.362269355915487,
            14.631125540472567,
            15.467019253410399,
            15.646071173250675,
            16.216_562_888_585_03,
            16.22252113185823,
            16.337_592_676_281_93,
            16.472833515144885,
            16.546852341853082,
            16.616179938428104,
            16.9334597280249,
            17.528_903_256_170_45,
            17.619_530_349_038_54,
            17.659138469025493,
            17.849509133957326,
            17.9549042833969,
            18.515_662_457_793_95,
            19.050154094584286,
            19.228_869_844_228_03,
            19.516776781529188,
            19.939716095104814,
        ];
        let y = [
            1.293428527408967,
            1.1000429951508104,
            0.985_676_209_295_923_4,
            1.078416645203697,
            1.0384497202430822,
            0.911_168_961_849_530_4,
            0.814_110_599_798_517_6,
            0.30870061707136665,
            0.360_630_501_091_990_5,
            2.653220107028643,
            0.035616391685330595,
            0.3396670541930909,
            0.1502442458597838,
            0.47763647905409656,
            0.094_209_347_547_079_92,
            0.10090502868850283,
            -0.012099541381990271,
            -0.018_456_884_031_812_94,
            0.027497693822005634,
            -0.139_844_621_530_572_6,
            0.012543821944120875,
            -0.003_995_997_315_695_388,
            -0.130_293_614_459_443_5,
            -0.46252856466785347,
            -0.19381646487771756,
            -0.35165581858660344,
            -0.320_150_797_023_739_7,
            -0.5361856364472154,
            -0.45587942674406273,
            -0.24861363699527775,
            -0.292_622_685_302_052_5,
            -0.28060451696229294,
            -0.11620781461766616,
            0.11487078623846304,
            0.26052335888772676,
            0.20269409789109122,
            0.431_968_733_064_165_2,
            0.276_224_130_937_391_6,
            0.47591651163197063,
            0.823_716_212_214_928_2,
            1.0238334201209505,
            1.7712082466756467,
            1.6076809524574995,
            2.032_190_109_082_845,
            1.7614737823535958,
            2.3284019188313416,
            2.4427408893291362,
            2.2803873792713096,
            2.3519177007491914,
            2.193_670_558_739_165,
            2.1263994596978226,
            2.3645007168762517,
            2.5249861145768056,
            2.122_272_855_834_001,
            0.299_203_443_460_228,
            2.172805960509645,
            2.3805956106200714,
            2.2635884407253357,
            2.1413513061524423,
            2.2424209481860458,
            1.8287138511852497,
            1.661_710_250_360_172,
            1.456_810_520_104_257,
            1.3638716209786697,
            1.3174847840211241,
            1.1294361324150146,
            1.6406879454416847,
            1.227_972_499_725_56,
            1.333044215401372,
            1.0429699942281332,
            0.830_352_777_457_079_7,
            0.594_504_552_950_756,
            0.864_889_768_902_684_3,
            1.0988525623839671,
            1.122_004_167_081_219,
            0.837_027_039_099_292_9,
            0.840_324_931_775_807_9,
            1.0662167330324928,
            0.7984426688798818,
            0.886_983_471_386_884_5,
        ];
        let w = [
            1.0889637747313827,
            1.4202703633345664,
            1.275_927_386_712_283,
            0.988_120_193_593_204,
            0.700_306_445_360_183_7,
            0.700_581_608_805_805_4,
            0.591_909_030_452_370_6,
            0.729152005398646,
            0.791_078_662_732_616_1,
            1.0620378700550646,
            1.2971524118911475,
            1.3690105613786727,
            0.929_895_969_340_577_7,
            0.623_054_652_707_651_3,
            1.2314715515822172,
            1.3028102195821702,
            1.3870370583608747,
            1.3259274356532842,
            1.3139970258343965,
            0.940_836_715_046_316_4,
            0.608_983_992_366_120_2,
            1.3987310056108981,
            0.797_507_378_039_881_6,
            0.788_213_320_076_465_6,
            0.765_858_371_742_069_7,
            0.582_622_988_615_185,
            0.907_211_071_578_785_8,
            1.3830148903653026,
            1.1364444366190583,
            1.4955437039025128,
            1.3725970827508718,
            1.314765295246616,
            1.106_141_550_699_249,
            0.869_393_355_911_597_6,
            0.928_705_062_018_707_4,
            0.695_141_761_098_057,
            0.962_294_715_689_495_2,
            0.610_878_398_176_282_6,
            1.4070586825255305,
            0.691_825_973_335_653_5,
            1.4032431468367577,
            1.4545864365063608,
            0.595_691_040_623_933_1,
            0.695_055_700_605_735_2,
            1.334_562_980_569_899,
            1.337_477_844_208_479,
            0.969_861_233_141_273_3,
            0.6534878071397543,
            1.439_056_930_365_041,
            1.483_104_110_462_591,
            1.316_310_929_134_488,
            0.9697276595979929,
            0.913_404_341_088_607_9,
            1.1443348934408277,
            0.999_848_530_394_956_5,
            1.0565558585803956,
            0.955_150_888_534_262_8,
            0.8411982893012464,
            1.2106719859875739,
            1.0034676706418395,
            0.781_219_038_413_837_6,
            0.896_646_452_136_337_8,
            1.325254688039422,
            0.828_120_636_753_737_9,
            1.2752281511202455,
            1.378995991544798,
            1.4765544817782938,
            0.808_399_976_463_988_4,
            1.0114127546548843,
            0.702_006_856_212_392_4,
            1.4021367505192757,
            1.4634487400762737,
            0.856_994_667_090_475_6,
            0.596_906_685_270_369,
            1.4358097377698869,
            1.249910962767899,
            1.3429545727558434,
            1.313_349_277_479_574,
            0.515_121_111_180_633_3,
            0.502_998_038_660_734_9,
        ];
        let out = weighted_lowess(&x, &y, Some(&w), 0.3, 4, 30, None);
        assert!(
            rel(out.delta, 0.593_369_432_979_009_3) < 1e-12,
            "delta {}",
            out.delta
        );

        let r_fitted = [
            1.1884821422349683,
            1.1401746805937334,
            1.1063703912870086,
            1.1026672684736833,
            0.872_240_929_345_854_8,
            0.855_951_874_762_393_3,
            0.5666943700277064,
            0.510_552_240_462_768_4,
            0.369_016_523_550_517_2,
            0.304_924_376_153_817_2,
            0.227_306_278_971_063_6,
            0.215_647_194_621_015_9,
            0.20552696445334262,
            0.20146372618232533,
            0.078_024_062_961_633_33,
            0.057_933_325_377_075_74,
            0.040_195_331_914_130_16,
            0.028_136_656_740_725_74,
            0.025_366_080_684_971_15,
            -0.050_912_098_098_715,
            -0.067_616_092_559_898_98,
            -0.068_931_638_772_633_97,
            -0.069_980_935_260_592_24,
            -0.11951181663114796,
            -0.1295277526928138,
            -0.21936636695564127,
            -0.262_508_751_761_269_3,
            -0.261_180_486_021_105_6,
            -0.260_946_821_212_595_9,
            -0.26000726252193695,
            -0.167_592_176_158_451_8,
            -0.050_584_694_681_363_06,
            0.00809634544820792,
            0.21510888843295684,
            0.33873746515295533,
            0.351_705_758_469_647,
            0.580_491_217_858_042_1,
            0.6082580657896921,
            0.723_588_826_350_821_9,
            0.771_095_359_716_515_1,
            0.924_832_205_407_716,
            1.4089433518016659,
            1.4327444779344392,
            1.9245662627937583,
            1.940_370_302_694_94,
            2.0169068927703515,
            2.135_858_147_609_854,
            2.1450068314468513,
            2.1483285077940617,
            2.157_513_074_233_857,
            2.1588144992196634,
            2.159_273_556_147_39,
            2.160083404651763,
            2.1792561271649924,
            2.1677064401616297,
            2.1472903708612345,
            2.1212760240350104,
            2.0965598171906232,
            2.065_629_550_176_702,
            1.9470028929533285,
            1.875434429262878,
            1.6244677094067717,
            1.5664639693873488,
            1.3816535905374199,
            1.3801654339673561,
            1.3514246681057358,
            1.3176463321656406,
            1.299_159_070_363_2,
            1.2818435104574863,
            1.2025983400628886,
            1.1018567306289988,
            1.0873066282823167,
            1.0809475784515046,
            1.0503837308105375,
            1.0334626295961167,
            0.943_433_383_693_995,
            0.861_571_735_904_992_1,
            0.834_199_996_583_610_7,
            0.792_602_471_709_358_3,
            0.731_495_126_462_035_3,
        ];
        for (a, b) in out.fitted.iter().zip(r_fitted.iter()) {
            close(*a, *b);
        }
        let sum: f64 = out.fitted.iter().sum();
        assert!(rel(sum, 73.941_169_255_579_42) < 1e-9, "sum {sum}");
    }

    /// Reference: n=15 with **unsorted** x — exercises the wrapper's internal
    /// ordering and the unsort of fitted/residuals back to input order.
    /// `span=0.4, iterations=3L`.
    #[test]
    fn weighted_lowess_case_d_unsorted() {
        let x = [
            2.4769279826432467,
            6.676_589_386_537_671,
            4.878_524_420_782_924,
            5.919_815_640_896_559,
            5.197_473_073_378_205,
            7.844762085005641,
            4.218_711_268_156_767,
            3.7431930266320705,
            2.252_489_222_213_626,
            2.2300906758755445,
            7.2444714065641165,
            6.980_926_280_841_231,
            0.829_540_582_373_738_3,
            0.291_089_225_560_426_7,
            1.8692466896027327,
        ];
        let y = [
            3.1588662152772407,
            5.265_610_750_585_987,
            4.479_077_038_757_614,
            4.999_896_114_465_269,
            4.528_259_545_777_066,
            5.9505529060869105,
            3.9898370326344743,
            3.837_755_420_730_662,
            2.957781766496248,
            3.0561989716819555,
            5.562_200_231_289_714,
            5.482_773_635_299_048,
            2.5100552062553017,
            2.106_936_485_716_76,
            2.938_948_707_868_278,
        ];
        let w = [
            1.3902532287407665,
            0.648_387_671_448_290_3,
            0.560_372_849_553_823_4,
            1.3949049064423888,
            0.839_275_473_868_474_2,
            1.0821476810146122,
            1.1763596059288828,
            0.8876171682029963,
            1.2749323339201508,
            0.48667634087614714,
            0.887_075_331_481_173_7,
            1.6853833251167087,
            1.1941875480115414,
            1.34794178516604,
            1.1733048296067863,
        ];
        let out = weighted_lowess(&x, &y, Some(&w), 0.4, 3, 200, None);
        assert_eq!(out.delta, 0.0);

        let r_fitted = [
            3.1610682329090687,
            5.328_131_047_276_528,
            4.401525118617398,
            4.953068229283593,
            4.574_573_132_447_472,
            5.921_881_766_698_609,
            4.045_695_517_076_606,
            3.8024012703068983,
            3.060_630_080_369_762,
            3.0500104995362616,
            5.6113672538311326,
            5.477_972_596_741_132,
            2.4294025367732424,
            2.142_734_808_562_499,
            2.907_051_913_114_72,
        ];
        let r_residuals = [
            -0.002_202_017_631_828_035,
            -0.062_520_296_690_541_6,
            0.077_551_920_140_216_17,
            0.04682788518167591,
            -0.046_313_586_670_406_08,
            0.028671139388301015,
            -0.055858484442131484,
            0.035_354_150_423_763_55,
            -0.10284831387351412,
            0.006_188_472_145_693_957,
            -0.049_167_022_541_418_25,
            0.004_801_038_557_915_582,
            0.080_652_669_482_059_29,
            -0.03579832284573925,
            0.031_896_794_753_558_22,
        ];
        for (a, b) in out.fitted.iter().zip(r_fitted.iter()) {
            close(*a, *b);
        }
        for (a, b) in out.residuals.iter().zip(r_residuals.iter()) {
            assert!((a - b).abs() < 1e-9, "{a} vs {b}");
        }
    }

    // ---- loessFit (the weightedLowess consumer) ----------------------------
    // Reference cases L1-L5 from `scratch/loessfit_ref.R` (limma 3.68.3, R
    // 4.6.0). The wrapper is called with R's defaults: min.weight=1e-5,
    // max.weight=1e5, equal.weights.as.null=TRUE.

    /// Compare against an R reference vector, treating `NaN` as the dropped
    /// non-finite observations limma returns as `NA`.
    fn close_vec(got: &[f64], want: &[f64]) {
        assert_eq!(got.len(), want.len(), "length mismatch");
        for (a, b) in got.iter().zip(want.iter()) {
            if b.is_nan() {
                assert!(a.is_nan(), "expected NA, got {a}");
            } else {
                close(*a, *b);
            }
        }
    }

    /// Residuals equal `y - fitted` in every loessFit path; verify that
    /// structurally (NaN-aware) rather than via a relative tolerance, which
    /// would amplify rounding on the near-zero residuals.
    fn check_residuals_are_y_minus_fitted(out: &LoessFit, y: &[f64]) {
        for ((r, yy), f) in out.residuals.iter().zip(y.iter()).zip(out.fitted.iter()) {
            let expected = yy - f;
            if expected.is_nan() {
                assert!(r.is_nan(), "expected NA residual, got {r}");
            } else {
                assert!((r - expected).abs() < 1e-12, "{r} vs {expected}");
            }
        }
    }

    /// L1: weighted (n=40, span=0.3, iter=4) with one non-finite covariate
    /// (x[14]) and one non-finite response (y[27]); a planted outlier at
    /// index 6 keeps its full residual. Exercises the weightedLowess path.
    #[test]
    fn loess_fit_l1_weighted_with_na() {
        let nan = f64::NAN;
        let x = [
            0.005_183_129_105_716_944,
            0.14047908363863826,
            0.16214577946811914,
            0.578_589_488_286_525,
            0.619_884_438_347_071_4,
            0.646_897_766_273_468_7,
            0.864_958_912_134_170_5,
            1.0355600714683533,
            1.2321620131842792,
            1.250_832_553_487_271,
            1.576_602_489_221_841,
            1.7511292267590761,
            2.024_130_343_925_208,
            2.0371259469538927,
            nan,
            2.7687330706976354,
            2.7724979422055185,
            2.899_750_091_601_163,
            3.276793584227562,
            3.3061099448241293,
            3.5034838342107832,
            3.637_880_280_148_238,
            3.981_494_423_933_327,
            3.991_432_855_837_047,
            4.117_796_053_178_608,
            4.407_502_717_804_164,
            4.801_340_533_886_105,
            4.817_655_358_929_187,
            4.830_177_505_500_615,
            5.0491762650199234,
            5.1060837297700346,
            5.7368572149425745,
            6.420_319_871_976_972,
            6.8036738759838045,
            7.339_874_927_420_169,
            8.126_682_038_418_949,
            8.510_418_669_320_643,
            8.806_991_728_488_356,
            9.071_829_735_767_096,
            9.548_492_254_689_336,
        ];
        let y = [
            -0.13028379272342372,
            0.16494167307661872,
            0.10534440906646064,
            0.7330729492059922,
            0.719_552_107_645_758_6,
            0.7335251729085972,
            3.896_540_372_057_267,
            0.9141200377945875,
            1.1454302634920055,
            1.0026956063690007,
            1.0944468337004631,
            1.1463798418508317,
            1.4395428355871485,
            0.985_150_160_028_773_8,
            0.887_350_237_670_555_2,
            1.0149620403673598,
            0.878_060_307_279_870_1,
            1.1285529306977484,
            0.398_293_292_094_334_6,
            0.42789453241135494,
            0.019349953543560616,
            0.255_488_015_807_366_5,
            0.230_069_869_085_910_2,
            -0.12743829720812616,
            0.17319277852114137,
            -0.14115542663394637,
            -0.473_134_585_337_793_2,
            nan,
            0.11773586387606913,
            0.10998459160176483,
            0.255_692_811_653_590_8,
            0.581_819_342_086_163_4,
            1.2570686237489412,
            1.9579856745571356,
            2.370_545_808_258_101,
            2.6969087938082077,
            2.462_926_449_216_594,
            2.4283890102220855,
            2.457_605_861_523_145,
            1.7983325143797855,
        ];
        let w = [
            0.636_493_015_708_401_8,
            0.791_828_367_765_992_9,
            1.9835953457048163,
            0.589_909_437_927_417_5,
            1.0680007479852065,
            0.42237004640046505,
            0.34347869313787666,
            1.7559089590795338,
            1.4984269458800554,
            1.372_503_425_157_629,
            1.9231961047044024,
            1.4591313774231822,
            0.664_618_161_669_932_3,
            1.1917643264634534,
            0.5399193662917241,
            1.1196884553181008,
            0.761_779_755_516_909,
            1.5054568398045376,
            1.6169831238687038,
            1.988671691226773,
            1.4473249244736508,
            1.2495072918944061,
            1.0664991330821068,
            1.4697452861350029,
            0.678_594_810_375_943_8,
            1.1991494379937648,
            1.3034791655139997,
            0.862_886_963_831_260_8,
            1.7764983414905147,
            0.817_129_175_551_235_5,
            1.2557559457607568,
            1.599_651_281_326_078,
            0.717_321_668_099_612,
            1.6017108282074333,
            1.4470027274917812,
            1.3925923728384078,
            0.542_866_303_306_073,
            0.613_774_243_392_981_6,
            1.4151965341996402,
            1.2633167420513927,
        ];
        let out = loess_fit(&y, &x, Some(&w), 0.3, 4, 1e-5, 1e5, true);
        let r_fitted = [
            0.0020154160943804374,
            0.131_892_620_268_470_8,
            0.15244602671480564,
            0.537_033_400_554_882_2,
            0.5749773054397006,
            0.599_992_286_147_306_4,
            0.820_962_205_024_607_4,
            0.922_188_649_092_287,
            0.999_131_539_751_303_6,
            1.0042360193467927,
            1.084424219541861,
            1.105226618879442,
            1.0789358441377175,
            1.0790114350851858,
            nan,
            0.858_709_557_143_218_8,
            0.856_107_782_666_883_5,
            0.8347718825091448,
            0.505_231_802_313_234_7,
            0.47837184475249783,
            0.295_909_335_228_325_1,
            0.20354399729824024,
            0.054_388_087_732_097_47,
            0.050_612_262_098_252_09,
            0.016304620106176948,
            0.041_451_857_739_025_66,
            0.10396584158117128,
            nan,
            0.1157231185509553,
            0.20740861707140956,
            0.23703826228746294,
            0.760_303_516_697_105_9,
            1.419_606_559_087_084,
            1.7266927821338722,
            2.0077981459480405,
            2.3279293123460407,
            2.475_834_436_511_091,
            2.5825278346690603,
            2.671_725_980_932_141,
            2.819_077_988_722_33,
        ];
        close_vec(&out.fitted, &r_fitted);
        check_residuals_are_y_minus_fitted(&out, &y);
        // The outlier at index 6 keeps its full residual.
        close(out.residuals[6], 3.075_578_167_032_66);
    }

    /// L2: too few obs for a smoother (n=6 < 4 + 1/0.3), so loessFit drops to a
    /// weighted straight-line fit. Exercises `weighted_line_fit`.
    #[test]
    fn loess_fit_l2_wls_branch() {
        let x = [
            1.5213841770309955,
            2.3736945423297584,
            2.6032693102024496,
            3.6165727430488914,
            4.216_155_161_848_292,
            4.967628859449178,
        ];
        let y = [
            2.1528132956168466,
            2.598_898_371_290_544,
            3.3797162340845226,
            3.511_792_998_604_96,
            3.9024791277875033,
            4.417381997642738,
        ];
        let w = [
            1.265_375_897_428_021,
            0.920_177_808_171_138_1,
            0.711_445_754_114_538_5,
            0.904_493_984_859_436_8,
            1.1457260956522077,
            1.515057343104854,
        ];
        let out = loess_fit(&y, &x, Some(&w), 0.3, 4, 1e-5, 1e5, true);
        let r_fitted = [
            2.221282436297829,
            2.765_360_343_792_229,
            2.911_910_901_433_168,
            3.558_759_788_700_577,
            3.941_507_155_238_185,
            4.421215315052069,
        ];
        close_vec(&out.fitted, &r_fitted);
        check_residuals_are_y_minus_fitted(&out, &y);
    }

    /// L3: all-equal weights collapse to NULL, so loessFit takes the classic
    /// `lowess` path. n=20, span=0.4, iter=4.
    #[test]
    fn loess_fit_l3_equal_weights_classic() {
        let x = [
            0.12413568049669266,
            0.943_929_295_986_890_8,
            1.0861205589026213,
            1.8004096411168575,
            2.0838854778558016,
            2.739089785143733,
            2.7455857265740633,
            3.097_543_876_618_147,
            3.157_202_487_811_446,
            3.4970000348985195,
            3.5675238333642483,
            3.869_830_997_660_756,
            4.138_796_932_995_319,
            5.527_887_191_623_449,
            6.197_990_983_724_594,
            6.255_103_558_301_926,
            6.7459734827280045,
            6.751_051_519_066_095,
            7.202_860_610_559_583,
            7.351_007_658_988_237,
        ];
        let y = [
            1.0298211176564178,
            0.676_610_752_805_095_6,
            0.855_910_305_007_615_9,
            -0.010864869491520768,
            -0.17715935228426227,
            -0.6484373184623734,
            -0.6695555621419208,
            -0.6411756350457184,
            -0.665_807_808_281_218_5,
            -0.677_004_584_491_505_2,
            -0.6202612231020217,
            -0.315_789_402_713_021_6,
            -0.020225008839295677,
            1.349_792_027_254_391,
            1.6439380305192044,
            1.6601868401446462,
            1.6584092984580048,
            1.8678128119816537,
            1.0507848512764517,
            1.087775447878236,
        ];
        let w = [2.0; 20];
        let out = loess_fit(&y, &x, Some(&w), 0.4, 4, 1e-5, 1e5, true);
        let r_fitted = [
            1.131_683_235_813_225,
            0.622_701_981_921_744,
            0.5340726402959507,
            0.0534601137616516,
            -0.19037724883930937,
            -0.581_662_324_318_019_5,
            -0.583_783_621_386_261_8,
            -0.646_813_245_039_454_5,
            -0.640_248_250_804_955_1,
            -0.569_616_038_400_150_6,
            -0.517_887_741_556_381_1,
            -0.30168620597235857,
            -0.094_559_894_415_800_07,
            1.262_038_052_285_814,
            1.4901527265313739,
            1.4785335054830242,
            1.3616342992329014,
            1.360_321_175_900_416,
            1.2346611633429285,
            1.1890588683115497,
        ];
        close_vec(&out.fitted, &r_fitted);
        check_residuals_are_y_minus_fitted(&out, &y);
    }

    /// L4: weights with one NA and two out-of-range values that get clamped to
    /// [1e-5, 1e5]; the 1e5-weighted point pins the curve to its y (zero
    /// residual). n=12, span=0.5, iter=3. Weighted path.
    #[test]
    fn loess_fit_l4_clamped_weights() {
        let nan = f64::NAN;
        let x = [
            0.000_812_362_879_514_694_2,
            0.19315525190904737,
            0.341_314_423_363_655_8,
            0.550_810_409_244_149_9,
            0.691_795_941_442_251_2,
            0.860_569_709_446_281_2,
            1.6972281779162586,
            2.683_141_722_343_862,
            2.6952867256477475,
            3.0456017875112593,
            3.4599765236489475,
            4.460_506_583_098_322,
        ];
        let y = [
            0.18370482724257098,
            0.063_986_626_702_510_85,
            -0.14976651732047797,
            0.1083635469696469,
            0.418_286_167_062_277_7,
            0.365_527_766_084_957_4,
            0.7026405120499789,
            1.2006499392960333,
            1.2752066988779682,
            1.5496202606111031,
            2.0578605625298567,
            1.9503358137564917,
        ];
        let w = [
            1.4987861497793347,
            0.835_285_151_377_320_3,
            nan,
            0.979_024_216_765_537_9,
            1.0222810602281243,
            0.813_572_955_550_625_9,
            1.1607185169123113,
            1e-08,
            0.676_338_533_638_045_2,
            10000000.0,
            0.755_478_185_601_532_5,
            0.796_876_986_976_712_9,
        ];
        let out = loess_fit(&y, &x, Some(&w), 0.5, 3, 1e-5, 1e5, true);
        let r_fitted = [
            0.092_715_644_709_798_5,
            0.1571141562155525,
            0.20738759663288137,
            0.27930367717161875,
            0.328_256_797_435_878_9,
            0.38783005712437973,
            0.730_831_473_478_615_7,
            1.2752055964112556,
            1.2752055966734654,
            1.5496202606111031,
            2.0578605625298567,
            1.9503358137564917,
        ];
        close_vec(&out.fitted, &r_fitted);
        check_residuals_are_y_minus_fitted(&out, &y);
        // The 1e5-weighted point (index 9) is fitted exactly: zero residual.
        close(out.residuals[9], 0.0);
    }

    /// L5: span (0.05) below 1/nobs (1/10) short-circuits to interpolation, so
    /// fitted == y and residuals are exactly zero regardless of weights.
    #[test]
    fn loess_fit_l5_tiny_span_identity() {
        let x: Vec<f64> = (1..=10).map(|i| i as f64).collect();
        let y: Vec<f64> = (1..=10).map(|i| (i as f64).powf(1.5)).collect();
        let w = vec![1.0; 10];
        let out = loess_fit(&y, &x, Some(&w), 0.05, 4, 1e-5, 1e5, true);
        for (f, yy) in out.fitted.iter().zip(y.iter()) {
            assert_eq!(f, yy);
        }
        for r in &out.residuals {
            assert_eq!(*r, 0.0);
        }
    }
}