limma-rust 0.1.0

//! normexp convolution model and background correction.
//!
//! Pure-Rust port of limma's `background-normexp.R`, `background.R` and the
//! saddlepoint optimiser in `src/normexp.c`. Implements:
//!
//! - [`normexp_signal`] — exact expected signal under the normal+exponential
//!   convolution model (`normexp.signal`).
//! - [`normexp_fit_saddle`] — `normexp.fit(method = "saddle")`, the limma
//!   default: minimise the saddlepoint approximation to `-2 log L` by
//!   Nelder-Mead ([`nmmin`], a faithful port of R's `nmmin`).
//! - [`background_correct_matrix`] — `backgroundCorrect.matrix` for the
//!   `none`/`subtract`/`half`/`minimum`/`movingmin`/`edwards`/`normexp` methods.

use crate::special::ln_norm_cdf;
use ndarray::Array2;
use std::f64::consts::PI;

/// `ln(2*pi)`.
const LN_2PI: f64 = 1.837_877_066_409_345_3;

/// Fitted parameters from [`normexp_fit_saddle`].
#[derive(Clone, Debug)]
pub struct NormexpFit {
    /// `(mu, log(sigma), log(alpha))`.
    pub par: [f64; 3],
    /// Minimised value of `-2 log L` (saddlepoint approximation).
    pub m2loglik: f64,
}

/// Background-correction method for [`background_correct_matrix`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum BackgroundMethod {
    None,
    Subtract,
    Half,
    Minimum,
    MovingMin,
    Edwards,
    Normexp,
}

/// `normexp.signal(par, x)`: expected signal given foreground `x` under the
/// normal(`mu`,`sigma`) + exponential(`alpha`) convolution model, where
/// `par = (mu, log(sigma), log(alpha))`.
pub fn normexp_signal(par: &[f64; 3], x: &[f64]) -> Vec<f64> {
    let mu = par[0];
    let sigma = par[1].exp();
    let sigma2 = sigma * sigma;
    let alpha = par[2].exp();
    assert!(alpha > 0.0, "alpha must be positive");
    assert!(sigma > 0.0, "sigma must be positive");

    let mut signal: Vec<f64> = x
        .iter()
        .map(|&xi| {
            let mu_sf = xi - mu - sigma2 / alpha;
            let z = mu_sf / sigma;
            // dnorm(0, mean=mu.sf, sd=sigma, log=TRUE)
            let log_dnorm0 = -0.5 * LN_2PI - sigma.ln() - 0.5 * z * z;
            // pnorm(0, mean=mu.sf, sd=sigma, lower.tail=FALSE, log.p=TRUE) = ln Phi(z)
            let log_pupper = ln_norm_cdf(z);
            mu_sf + sigma2 * (log_dnorm0 - log_pupper).exp()
        })
        .collect();

    // If any (non-NaN) signal is negative, floor every non-NaN signal at 1e-6.
    let any_neg = signal.iter().any(|&s| !s.is_nan() && s < 0.0);
    if any_neg {
        for s in signal.iter_mut() {
            if !s.is_nan() {
                *s = s.max(1e-6);
            }
        }
    }
    signal
}

/// `normexp.fit(x, method = "saddle")`. Returns `(mu, log(sigma), log(alpha))`
/// and the minimised saddlepoint `-2 log L`.
pub fn normexp_fit_saddle(x: &[f64]) -> NormexpFit {
    let xv: Vec<f64> = x.iter().copied().filter(|v| !v.is_nan()).collect();
    assert!(
        xv.len() >= 4,
        "Not enough data: need at least 4 non-missing corrected intensities"
    );

    // Starting values for mu, sigma^2, alpha.
    let q = quantile_type7(&xv, &[0.0, 0.05, 0.10, 1.0]);
    if q[0] == q[3] {
        return NormexpFit {
            par: [q[0], f64::NEG_INFINITY, f64::NEG_INFINITY],
            m2loglik: f64::NAN,
        };
    }
    let mu = if q[1] > q[0] {
        q[1]
    } else if q[2] > q[0] {
        q[2]
    } else {
        q[0] + 0.05 * (q[3] - q[0])
    };
    let below: Vec<f64> = xv.iter().copied().filter(|&v| v < mu).collect();
    let sigma2 = below.iter().map(|&v| (v - mu) * (v - mu)).sum::<f64>() / below.len() as f64;
    let mut alpha = xv.iter().sum::<f64>() / xv.len() as f64 - mu;
    if alpha <= 0.0 {
        alpha = 1e-6;
    }
    let par0 = [mu, 0.5 * sigma2.ln(), alpha.ln()];

    // Minimise the saddlepoint -2logL by Nelder-Mead, matching limma's
    // fit_saddle_nelder_mead: alpha=1, beta=0.5, gamma=2, intol=sqrt(eps),
    // abstol=-Inf, maxit=500.
    let (par, m2loglik) = nmmin(
        &par0,
        |p| normexp_m2loglik_saddle(p, &xv),
        -1e308,
        1.490_116e-08,
        1.0,
        0.5,
        2.0,
        500,
    );
    NormexpFit {
        par: [par[0], par[1], par[2]],
        m2loglik,
    }
}

/// Saddlepoint approximation to `-2 log L` as a function of
/// `par = (mu, log(sigma), log(alpha))`. Direct port of
/// `normexp_m2loglik_saddle` in `src/normexp.c`.
fn normexp_m2loglik_saddle(par: &[f64], x: &[f64]) -> f64 {
    let mu = par[0];
    let sigma = par[1].exp();
    let sigma2 = sigma * sigma;
    let alpha = par[2].exp();
    let alpha2 = alpha * alpha;
    let alpha3 = alpha * alpha2;
    let alpha4 = alpha2 * alpha2;
    let n = x.len();
    let c2 = sigma2 * alpha;

    let mut upperbound = vec![0.0; n];
    let mut theta = vec![0.0; n];
    let mut has_converged = vec![false; n];

    for i in 0..n {
        let err = x[i] - mu;
        let upperbound1 = ((err - alpha) / (alpha * err.abs())).max(0.0);
        let upperbound2 = err / sigma2;
        upperbound[i] = upperbound1.min(upperbound2);
        let c1 = -sigma2 - err * alpha;
        let c0 = -alpha + err;
        let theta_quadratic = (-c1 - (c1 * c1 - 4.0 * c0 * c2).sqrt()) / (2.0 * c2);
        theta[i] = theta_quadratic.min(upperbound[i]);
    }

    // Globally convergent Newton iteration per point.
    let mut j = 0;
    let mut n_converged = 0usize;
    loop {
        j += 1;
        for i in 0..n {
            if has_converged[i] {
                continue;
            }
            let omat = 1.0 - alpha * theta[i];
            let dk = mu + sigma2 * theta[i] + alpha / omat;
            let ddk = sigma2 + alpha2 / (omat * omat);
            let delta = (x[i] - dk) / ddk;
            theta[i] += delta;
            if j == 1 {
                theta[i] = theta[i].min(upperbound[i]);
            }
            if delta.abs() < 1e-10 {
                has_converged[i] = true;
                n_converged += 1;
            }
        }
        if n_converged == n || j > 50 {
            break;
        }
    }

    let mut loglik = 0.0;
    for i in 0..n {
        let omat = 1.0 - alpha * theta[i];
        let omat2 = omat * omat;
        let k1 = mu * theta[i] + 0.5 * sigma2 * theta[i] * theta[i] - omat.ln();
        let k2 = sigma2 + alpha2 / omat2;
        let mut logf = -0.5 * (2.0 * PI * k2).ln() - x[i] * theta[i] + k1;
        let k3 = 2.0 * alpha3 / (omat * omat2);
        let k4 = 6.0 * alpha4 / (omat2 * omat2);
        logf += k4 / (8.0 * k2 * k2) - (5.0 * k3 * k3) / (24.0 * k2 * k2 * k2);
        loglik += logf;
    }
    -2.0 * loglik
}

/// Faithful port of R's `nmmin` (`src/appl/optim.c`) Nelder-Mead minimiser.
/// Returns the best vertex and its objective value. Index-based loops and the
/// argument list mirror the C source so the port stays auditable against it.
#[allow(clippy::too_many_arguments, clippy::needless_range_loop)]
fn nmmin(
    start: &[f64],
    objective: impl Fn(&[f64]) -> f64,
    abstol: f64,
    intol: f64,
    refl: f64,
    contract: f64,
    extend: f64,
    maxit: i32,
) -> (Vec<f64>, f64) {
    const BIG: f64 = 1.0e35;
    let n = start.len();
    let n1 = n + 1; // number of simplex vertices; row index of f-values
    let c = n + 2; // scratch (centroid) column index is c-1

    // P[row][col]: rows 0..n-1 hold coords, row n1-1 holds f-values;
    // columns 0..n1-1 are vertices, column c-1 is the centroid scratch.
    let mut p = vec![vec![0.0_f64; c]; n1];
    let mut bvec = start.to_vec();

    let f0 = objective(&bvec);
    let mut funcount = 1i32;
    let convtol = intol * (f0.abs() + intol);
    p[n1 - 1][0] = f0;
    for i in 0..n {
        p[i][0] = bvec[i];
    }

    let mut l = 1usize;
    let mut size = 0.0;
    let mut step = 0.0;
    for i in 0..n {
        if 0.1 * bvec[i].abs() > step {
            step = 0.1 * bvec[i].abs();
        }
    }
    if step == 0.0 {
        step = 0.1;
    }

    let mut vl;
    let mut vh;
    let mut h;
    let mut shrinkfail = false;

    'outer: loop {
        // BUILD the simplex around the current point bvec.
        for j in 2..=n1 {
            for i in 0..n {
                p[i][j - 1] = bvec[i];
            }
            let mut trystep = step;
            while p[j - 2][j - 1] == bvec[j - 2] {
                p[j - 2][j - 1] = bvec[j - 2] + trystep;
                trystep *= 10.0;
            }
            size += trystep;
        }
        let mut oldsize = size;
        let mut calcvert = true;

        loop {
            if calcvert {
                for j in 0..n1 {
                    if j + 1 != l {
                        for i in 0..n {
                            bvec[i] = p[i][j];
                        }
                        let mut fv = objective(&bvec);
                        if !fv.is_finite() {
                            fv = BIG;
                        }
                        funcount += 1;
                        p[n1 - 1][j] = fv;
                    }
                }
                calcvert = false;
            }

            // Find lowest (L) and highest (H) vertices.
            vl = p[n1 - 1][l - 1];
            vh = vl;
            h = l;
            for j in 1..=n1 {
                if j != l {
                    let fj = p[n1 - 1][j - 1];
                    if fj < vl {
                        l = j;
                        vl = fj;
                    }
                    if fj > vh {
                        h = j;
                        vh = fj;
                    }
                }
            }

            if vh > vl + convtol && vl > abstol {
                // Centroid of all vertices except H, into column c-1.
                for i in 0..n {
                    let mut temp = -p[i][h - 1];
                    for j in 0..n1 {
                        temp += p[i][j];
                    }
                    p[i][c - 1] = temp / n as f64;
                }
                // Reflection.
                for i in 0..n {
                    bvec[i] = (1.0 + refl) * p[i][c - 1] - refl * p[i][h - 1];
                }
                let mut vr = objective(&bvec);
                if !vr.is_finite() {
                    vr = BIG;
                }
                funcount += 1;

                if vr < vl {
                    // Try extension.
                    p[n1 - 1][c - 1] = vr;
                    for i in 0..n {
                        let fe = extend * bvec[i] + (1.0 - extend) * p[i][c - 1];
                        p[i][c - 1] = bvec[i];
                        bvec[i] = fe;
                    }
                    let mut fe = objective(&bvec);
                    if !fe.is_finite() {
                        fe = BIG;
                    }
                    funcount += 1;
                    if fe < vr {
                        for i in 0..n {
                            p[i][h - 1] = bvec[i];
                        }
                        p[n1 - 1][h - 1] = fe;
                    } else {
                        for i in 0..n {
                            p[i][h - 1] = p[i][c - 1];
                        }
                        p[n1 - 1][h - 1] = vr;
                    }
                } else {
                    if vr < vh {
                        // Replace worst with the reflection (lo-reduction).
                        for i in 0..n {
                            p[i][h - 1] = bvec[i];
                        }
                        p[n1 - 1][h - 1] = vr;
                    }
                    // Contraction toward the (possibly updated) worst vertex.
                    for i in 0..n {
                        bvec[i] = (1.0 - contract) * p[i][h - 1] + contract * p[i][c - 1];
                    }
                    let mut fc = objective(&bvec);
                    if !fc.is_finite() {
                        fc = BIG;
                    }
                    funcount += 1;
                    if fc < p[n1 - 1][h - 1] {
                        for i in 0..n {
                            p[i][h - 1] = bvec[i];
                        }
                        p[n1 - 1][h - 1] = fc;
                    } else if vr >= vh {
                        // Shrink toward the best vertex L.
                        calcvert = true;
                        size = 0.0;
                        for j in 0..n1 {
                            if j + 1 != l {
                                for i in 0..n {
                                    p[i][j] = contract * (p[i][j] - p[i][l - 1]) + p[i][l - 1];
                                    size += (p[i][j] - p[i][l - 1]).abs();
                                }
                            }
                        }
                        if size < oldsize {
                            shrinkfail = false;
                            oldsize = size;
                        } else {
                            shrinkfail = true;
                        }
                    }
                }
            }

            if !(vh > vl + convtol && vl > abstol && !shrinkfail && funcount <= maxit) {
                break;
            }
        }

        // Rebuild around the best vertex if a shrink stalled (degenerate
        // simplex); otherwise we have converged or hit maxit.
        if shrinkfail && funcount <= maxit && vh > vl + convtol && vl > abstol {
            for i in 0..n {
                bvec[i] = p[i][l - 1];
            }
            shrinkfail = false;
            continue 'outer;
        }
        break 'outer;
    }

    let fmin = p[n1 - 1][l - 1];
    let best: Vec<f64> = (0..n).map(|i| p[i][l - 1]).collect();
    (best, fmin)
}

/// `backgroundCorrect.matrix(E, Eb, method, offset)`. `eb` is the optional
/// background matrix; when absent, the methods that require it fall back to
/// `none` (as limma does). `offset` is added after correction iff nonzero.
pub fn background_correct_matrix(
    e: &Array2<f64>,
    eb: Option<&Array2<f64>>,
    method: BackgroundMethod,
    offset: f64,
) -> Array2<f64> {
    use BackgroundMethod::*;

    let method = if eb.is_none() {
        match method {
            Subtract | Half | Minimum | MovingMin | Edwards => None,
            other => other,
        }
    } else {
        method
    };

    let mut out = match method {
        None => e.clone(),
        Subtract => e - eb.unwrap(),
        Half => (e - eb.unwrap()).mapv(|v| v.max(0.5)),
        Minimum => {
            let mut m = e - eb.unwrap();
            for mut col in m.columns_mut() {
                let lo = col
                    .iter()
                    .copied()
                    .filter(|&v| v >= 1e-18)
                    .fold(f64::INFINITY, f64::min);
                if lo.is_finite() {
                    for v in col.iter_mut() {
                        if *v < 1e-18 {
                            *v = lo / 2.0;
                        }
                    }
                }
            }
            m
        }
        MovingMin => e - &ma3x3_min(eb.unwrap()),
        Edwards => edwards(e, eb.unwrap()),
        Normexp => {
            let eb_sub = match eb {
                Some(bg) => e - bg,
                Option::None => e.clone(),
            };
            let ncol = eb_sub.ncols();
            // Each column's saddle-point MLE + signal is an independent, fairly
            // heavy optimisation. Compute them (in parallel under the `parallel`
            // feature, serially otherwise) and scatter into the output in column
            // order. The per-column arithmetic is untouched, so the result is
            // bit-identical regardless of feature or thread count.
            let solve = |j: usize| -> Vec<f64> {
                let x: Vec<f64> = eb_sub.column(j).to_vec();
                let fit = normexp_fit_saddle(&x);
                normexp_signal(&fit.par, &x)
            };
            #[cfg(feature = "parallel")]
            let cols: Vec<Vec<f64>> = {
                use rayon::prelude::*;
                (0..ncol).into_par_iter().map(solve).collect()
            };
            #[cfg(not(feature = "parallel"))]
            let cols: Vec<Vec<f64>> = (0..ncol).map(solve).collect();
            // Every cell is overwritten below, so reuse `eb_sub` rather than
            // cloning it.
            let mut m = eb_sub;
            for (j, sig) in cols.into_iter().enumerate() {
                for (i, s) in sig.into_iter().enumerate() {
                    m[[i, j]] = s;
                }
            }
            m
        }
    };

    if offset != 0.0 {
        out.mapv_inplace(|v| v + offset);
    }
    out
}

/// Edwards (2003) log-linear interpolation for dull spots.
fn edwards(e: &Array2<f64>, eb: &Array2<f64>) -> Array2<f64> {
    let sub = e - eb;
    let mut out = sub.clone();
    for (j, col) in sub.columns().into_iter().enumerate() {
        let d: Vec<f64> = col.to_vec();
        let frac = d.iter().filter(|&&v| v < 1e-16).count() as f64 / d.len() as f64;
        let prob = (frac * 1.1).min(1.0);
        let delta = quantile_type7(&d, &[prob])[0];
        for i in 0..d.len() {
            let s = sub[[i, j]];
            out[[i, j]] = if s < delta {
                delta * (1.0 - (eb[[i, j]] + delta) / e[[i, j]]).exp()
            } else {
                s
            };
        }
    }
    out
}

/// 3x3 moving minimum over a matrix, border-aware (NA padding dropped), as in
/// limma's `ma3x3.matrix(x, FUN = min)`.
fn ma3x3_min(x: &Array2<f64>) -> Array2<f64> {
    let (nr, nc) = x.dim();
    let mut out = Array2::<f64>::zeros((nr, nc));
    for r in 0..nr {
        for col in 0..nc {
            let mut m = f64::INFINITY;
            for dr in -1i64..=1 {
                for dc in -1i64..=1 {
                    let rr = r as i64 + dr;
                    let cc = col as i64 + dc;
                    if rr >= 0 && rr < nr as i64 && cc >= 0 && cc < nc as i64 {
                        let v = x[[rr as usize, cc as usize]];
                        if !v.is_nan() && v < m {
                            m = v;
                        }
                    }
                }
            }
            out[[r, col]] = m;
        }
    }
    out
}

/// Type-7 sample quantiles (R's `quantile` default).
fn quantile_type7(x: &[f64], probs: &[f64]) -> Vec<f64> {
    let mut s: Vec<f64> = x.iter().copied().filter(|v| !v.is_nan()).collect();
    s.sort_by(|a, b| a.partial_cmp(b).unwrap());
    let n = s.len();
    probs
        .iter()
        .map(|&p| {
            if n == 0 {
                return f64::NAN;
            }
            if n == 1 {
                return s[0];
            }
            let hpos = (n as f64 - 1.0) * p;
            let lo = hpos.floor() as usize;
            let hi = (lo + 1).min(n - 1);
            s[lo] + (hpos - lo as f64) * (s[hi] - s[lo])
        })
        .collect()
}

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

    // The 24-value foreground vector used by scratch/normexp_ref.R.
    fn xvec() -> Vec<f64> {
        vec![
            2.1, 3.5, 2.8, 10.2, 5.6, 4.1, 8.9, 3.3, 6.7, 12.5, 2.9, 4.8, 7.2, 3.1, 9.4, 5.0, 4.4,
            6.1, 3.8, 11.0, 2.5, 5.9, 7.7, 4.6,
        ]
    }

    fn emat() -> Array2<f64> {
        Array2::from_shape_vec((12, 2), {
            let x = xvec();
            // column-major matrix(x, 12, 2) -> row-major (12,2) needs transpose order.
            let mut v = vec![0.0; 24];
            for i in 0..12 {
                v[i * 2] = x[i];
                v[i * 2 + 1] = x[12 + i];
            }
            v
        })
        .unwrap()
    }

    fn ebmat() -> Array2<f64> {
        let col_major = [
            1.0, 4.0, 2.0, 3.0, 1.5, 4.5, 2.0, 3.5, 1.0, 2.0, 3.0, 5.0, 2.0, 3.5, 1.0, 2.0, 1.0,
            2.0, 4.0, 1.5, 3.0, 1.0, 2.0, 5.0,
        ];
        let mut v = vec![0.0; 24];
        for i in 0..12 {
            v[i * 2] = col_major[i];
            v[i * 2 + 1] = col_major[12 + i];
        }
        Array2::from_shape_vec((12, 2), v).unwrap()
    }

    fn col_major(m: &Array2<f64>) -> Vec<f64> {
        let (nr, nc) = m.dim();
        let mut v = Vec::with_capacity(nr * nc);
        for j in 0..nc {
            for i in 0..nr {
                v.push(m[[i, j]]);
            }
        }
        v
    }

    fn assert_close(a: &[f64], b: &[f64], tol: f64) {
        assert_eq!(a.len(), b.len(), "length mismatch");
        for (k, (&x, &y)) in a.iter().zip(b.iter()).enumerate() {
            assert!(
                (x - y).abs() <= tol + tol * y.abs(),
                "index {k}: got {x}, want {y} (diff {})",
                (x - y).abs()
            );
        }
    }

    #[test]
    fn normexp_fit_saddle_matches_r() {
        let fit = normexp_fit_saddle(&xvec());
        // Reference: normexp.fit(x, "saddle").
        assert!((fit.par[0] - 2.367_538_459_991_54).abs() < 1e-5);
        assert!((fit.par[1] - -1.030_352_153_905_64).abs() < 1e-5);
        assert!((fit.par[2] - 1.219_722_868_953_23).abs() < 1e-5);
        assert!((fit.m2loglik - 111.423_810_138_321).abs() < 1e-5);
    }

    #[test]
    fn normexp_signal_matches_r() {
        let par = [
            2.367_538_459_991_54,
            -1.030_352_153_905_64,
            1.219_722_868_953_23,
        ];
        let sig = normexp_signal(&par, &xvec());
        let want = [
            0.198163450107075,
            1.09613822556942,
            0.484025139654586,
            7.79484935368556,
            3.19484935368556,
            1.69485115668388,
            6.49484935368557,
            0.901027384489598,
            4.29484935368557,
            10.0948493536856,
            0.554205422821398,
            2.39484935370929,
            4.79484935368557,
            0.716807903958157,
            6.99484935368557,
            2.59484935368604,
            1.99484937706265,
            3.69484935368556,
            1.39491795493823,
            8.59484935368556,
            0.322091491194509,
            3.49484935368557,
            5.29484935368557,
            2.19484935455685,
        ];
        assert_close(&sig, &want, 1e-6);
    }

    #[test]
    fn background_correct_normexp_offsets_match_r() {
        let e = emat();
        let bc0 = background_correct_matrix(&e, None, BackgroundMethod::Normexp, 0.0);
        let want0 = [
            1.99403319184388e-07,
            1.40000019940332,
            0.700000199403319,
            8.10000019940332,
            3.50000019940332,
            2.00000019940332,
            6.80000019940332,
            1.20000019940332,
            4.60000019940332,
            10.4000001994033,
            0.800000199403319,
            2.70000019940332,
            3.38853053495406,
            0.614450546313025,
            5.58655635912811,
            1.43189569191768,
            1.07932597446681,
            2.32320650063882,
            0.821622365173466,
            7.18655615956755,
            0.492789993991968,
            2.142117281233,
            3.88691863546145,
            1.18562095914263,
        ];
        assert_close(&col_major(&bc0), &want0, 1e-6);

        let bc16 = background_correct_matrix(&e, None, BackgroundMethod::Normexp, 16.0);
        let want16: Vec<f64> = want0.iter().map(|v| v + 16.0).collect();
        assert_close(&col_major(&bc16), &want16, 1e-6);
    }

    #[test]
    fn background_correct_eb_methods_match_r() {
        let e = emat();
        let eb = ebmat();

        let sub = background_correct_matrix(&e, Some(&eb), BackgroundMethod::Subtract, 0.0);
        assert_close(
            &col_major(&sub),
            &[
                1.1, -0.5, 0.8, 7.2, 4.1, -0.4, 6.9, -0.2, 5.7, 10.5, -0.1, -0.2, 5.2, -0.4, 8.4,
                3.0, 3.4, 4.1, -0.2, 9.5, -0.5, 4.9, 5.7, -0.4,
            ],
            1e-12,
        );

        let half = background_correct_matrix(&e, Some(&eb), BackgroundMethod::Half, 0.0);
        assert_close(
            &col_major(&half),
            &[
                1.1, 0.5, 0.8, 7.2, 4.1, 0.5, 6.9, 0.5, 5.7, 10.5, 0.5, 0.5, 5.2, 0.5, 8.4, 3.0,
                3.4, 4.1, 0.5, 9.5, 0.5, 4.9, 5.7, 0.5,
            ],
            1e-12,
        );

        let mn = background_correct_matrix(&e, Some(&eb), BackgroundMethod::Minimum, 0.0);
        assert_close(
            &col_major(&mn),
            &[
                1.1, 0.4, 0.8, 7.2, 4.1, 0.4, 6.9, 0.4, 5.7, 10.5, 0.4, 0.4, 5.2, 1.5, 8.4, 3.0,
                3.4, 4.1, 1.5, 9.5, 1.5, 4.9, 5.7, 1.5,
            ],
            1e-12,
        );

        let ed = background_correct_matrix(&e, Some(&eb), BackgroundMethod::Edwards, 0.0);
        assert_close(
            &col_major(&ed),
            &[
                1.1,
                0.558422539017665,
                0.808880852322533,
                7.2,
                4.1,
                0.604489443607101,
                6.9,
                0.597824798774299,
                5.7,
                10.5,
                0.593157962311005,
                0.657982551965439,
                5.2,
                1.00197356183341,
                8.4,
                3.00530848233458,
                3.4,
                4.1,
                1.29360494333044,
                9.5,
                0.73912631360648,
                4.9,
                5.7,
                1.43478914370254,
            ],
            1e-9,
        );

        let mm = background_correct_matrix(&e, Some(&eb), BackgroundMethod::MovingMin, 0.0);
        assert_close(
            &col_major(&mm),
            &[
                1.1, 2.5, 1.8, 9.2, 4.6, 3.1, 7.4, 2.3, 5.7, 11.5, 1.9, 2.8, 6.2, 2.1, 8.4, 4.0,
                3.4, 5.1, 2.3, 10.0, 1.5, 4.9, 6.7, 2.6,
            ],
            1e-12,
        );

        let nx = background_correct_matrix(&e, Some(&eb), BackgroundMethod::Normexp, 0.0);
        let want_nx = [
            1.60000003144471,
            3.1444714387814e-08,
            1.30000003144471,
            7.70000003144471,
            4.60000003144471,
            0.100000031444714,
            7.40000003144471,
            0.300000031444714,
            6.20000003144471,
            11.0000000314447,
            0.400000031444714,
            0.300000031444714,
            5.70004586914821,
            0.100045869148213,
            8.90004586914821,
            3.50004586914821,
            3.90004586914821,
            4.60004586914821,
            0.300045869148213,
            10.0000458691482,
            4.58691482126749e-05,
            5.40004586914821,
            6.20004586914821,
            0.100045869148212,
        ];
        assert_close(&col_major(&nx), &want_nx, 1e-6);
    }
}