siplot 0.3.0

//! Basic curve fitting utilities.
//!
//! Provides traits and simple implementations for curve fitting (Linear, Gaussian estimation).
//!
//! Additionally provides an iterative Levenberg-Marquardt least-squares solver
//! ([`leastsq`]) ported from silx `silx/math/fit/leastsq.py`, together with the
//! peak models (Gaussian/Lorentzian/PseudoVoigt) from
//! `silx/math/fit/functions/src/funs.c` and their initial-parameter estimators
//! mirroring `silx/math/fit/fittheories.py`.

/// Result of a curve fit.
#[derive(Debug, Clone)]
pub struct FitResult {
    /// The fitted y values for the input x values.
    pub y_fit: Vec<f64>,
    /// The parameters of the fit function.
    pub parameters: Vec<f64>,
    /// Names of the parameters.
    pub param_names: Vec<String>,
}

/// A function that can be fitted to data.
pub trait FitFunction {
    /// Name of the function.
    fn name(&self) -> &str;

    /// Fit the function to the given data.
    fn fit(&self, x: &[f64], y: &[f64]) -> Option<FitResult>;
}

/// Simple linear fit: y = m*x + c
pub struct LinearFit;

impl FitFunction for LinearFit {
    fn name(&self) -> &str {
        "Linear"
    }

    fn fit(&self, x: &[f64], y: &[f64]) -> Option<FitResult> {
        if x.len() != y.len() || x.len() < 2 {
            return None;
        }
        let n = x.len() as f64;
        let sum_x: f64 = x.iter().sum();
        let sum_y: f64 = y.iter().sum();
        let sum_xy: f64 = x.iter().zip(y.iter()).map(|(&xi, &yi)| xi * yi).sum();
        let sum_xx: f64 = x.iter().map(|&xi| xi * xi).sum();

        let denominator = n * sum_xx - sum_x * sum_x;
        if denominator.abs() < 1e-12 {
            return None;
        }

        let m = (n * sum_xy - sum_x * sum_y) / denominator;
        let c = (sum_y - m * sum_x) / n;

        let y_fit = x.iter().map(|&xi| m * xi + c).collect();

        Some(FitResult {
            y_fit,
            parameters: vec![m, c],
            param_names: vec!["Slope (m)".to_string(), "Intercept (c)".to_string()],
        })
    }
}

/// Gaussian estimation: y = A * exp(-(x - mu)^2 / (2 * sigma^2)) + bg
/// Note: This is a direct analytical estimation based on moments/peak, not an iterative L-M fit.
pub struct GaussianEstimateFit;

impl FitFunction for GaussianEstimateFit {
    fn name(&self) -> &str {
        "Gaussian (Estimate)"
    }

    fn fit(&self, x: &[f64], y: &[f64]) -> Option<FitResult> {
        if x.len() != y.len() || x.len() < 3 {
            return None;
        }

        let bg = y.iter().copied().fold(f64::INFINITY, f64::min);
        let mut max_y = f64::NEG_INFINITY;
        let mut max_idx = 0;
        for (i, &yi) in y.iter().enumerate() {
            if yi > max_y {
                max_y = yi;
                max_idx = i;
            }
        }

        let a = max_y - bg;
        let mu = x[max_idx];

        // Estimate FWHM by finding first points below half max
        let half_max = bg + a / 2.0;
        let mut left_idx = max_idx;
        while left_idx > 0 && y[left_idx] > half_max {
            left_idx -= 1;
        }
        let mut right_idx = max_idx;
        while right_idx < y.len() - 1 && y[right_idx] > half_max {
            right_idx += 1;
        }

        let fwhm = x[right_idx] - x[left_idx];
        let sigma = if fwhm > 0.0 {
            fwhm / 2.355
        } else {
            (x.last().unwrap() - x.first().unwrap()) / 4.0
        };

        let y_fit = x
            .iter()
            .map(|&xi| {
                let z = (xi - mu) / sigma;
                a * (-0.5 * z * z).exp() + bg
            })
            .collect();

        Some(FitResult {
            y_fit,
            parameters: vec![a, mu, sigma, bg],
            param_names: vec![
                "Amplitude (A)".to_string(),
                "Center (mu)".to_string(),
                "Sigma".to_string(),
                "Background".to_string(),
            ],
        })
    }
}

// ---------------------------------------------------------------------------
// Iterative Levenberg-Marquardt least-squares core.
//
// Ported from silx `silx/math/fit/leastsq.py` (leastsq / chisq_alpha_beta),
// itself a refactor of PyMca Gefit. We port the *unconstrained* path: silx's
// CFREE branch where `n_free == nparameters`, `noigno == range(n)`, and
// `derivfactor == 1`. Constraints (positivity/quoted/factor/...) are DEFERRED.
// ---------------------------------------------------------------------------

/// `LOG2`, matching the C constant in `funs.c`
/// (`#define LOG2 0.69314718055994529`, i.e. `ln(2)`). Used to convert FWHM to
/// sigma: `sigma = fwhm / (2 * sqrt(2 * LOG2))`.
pub const LOG2: f64 = std::f64::consts::LN_2;

/// `2 * sqrt(2 * LOG2)`: the FWHM/sigma conversion factor for a Gaussian.
/// silx computes `inv_two_sqrt_two_log2 = 1 / (2*sqrt(2*LOG2))` and uses
/// `sigma = fwhm * inv_two_sqrt_two_log2`.
pub fn fwhm_to_sigma_factor() -> f64 {
    2.0 * (2.0 * LOG2).sqrt()
}

/// Outputs of a successful [`leastsq`] run.
///
/// Mirrors the silx `leastsq` return tuple (`fittedpar`, `cov`, `ddict`) with
/// the unconstrained-only subset of `ddict` we need: `chisq`, `reduced_chisq`,
/// `niter`, `nfev`.
#[derive(Debug, Clone)]
pub struct LeastSqResult {
    /// Optimal parameter values minimising the weighted sum of squared
    /// residuals (silx `fittedpar`).
    pub parameters: Vec<f64>,
    /// Estimated covariance matrix of the parameters, row-major
    /// `n_param x n_param` (silx `cov0 = inv(alpha0)`). Standard errors are the
    /// square roots of the diagonal: `perr[i] = sqrt(cov[i][i])`.
    pub covariance: Vec<Vec<f64>>,
    /// Per-parameter uncertainties propagated through the applied constraints
    /// (silx `ddict["uncertainties"]` via `_get_sigma_parameters`). For an
    /// unconstrained fit this equals [`LeastSqResult::std_errors`]; with
    /// constraints it additionally carries the QUOTED `B·cos` factor and ties
    /// FACTOR/DELTA/SUM uncertainties to their reference parameter.
    pub uncertainties: Vec<f64>,
    /// The chi-square `sum( weight * (model - y)^2 )` at the optimum
    /// (silx `chisq0`).
    pub chisq: f64,
    /// Reduced chi-square `chisq / (M - n_free)` where `M` is the number of
    /// data points and `n_free` the number of fitted parameters (silx
    /// `reduced_chisq`). `None` when degrees of freedom are non-positive.
    pub reduced_chisq: Option<f64>,
    /// Number of iterations performed (silx `niter`).
    pub niter: usize,
    /// Number of model function evaluations (silx `nfev`).
    pub nfev: usize,
}

impl LeastSqResult {
    /// Per-parameter standard error: `sqrt(abs(cov[i][i]))`.
    ///
    /// Mirrors the silx docstring note "To compute one standard deviation
    /// errors use `perr = np.sqrt(np.diag(pcov))`"; `abs` guards a tiny
    /// negative diagonal from round-off, matching silx `sqrt(abs(diag(cov0)))`.
    pub fn std_errors(&self) -> Vec<f64> {
        (0..self.parameters.len())
            .map(|i| self.covariance[i][i].abs().sqrt())
            .collect()
    }
}

/// Why a [`leastsq`] call could not run / converge to a covariance.
#[derive(Debug, Clone, PartialEq, Eq)]
pub enum FitError {
    /// `xdata` and `ydata` have different lengths.
    LengthMismatch,
    /// There are no free parameters (silx `raise ValueError("No free
    /// parameters to fit")`).
    NoFreeParameters,
    /// Fewer data points than free parameters: the problem is under-determined.
    NotEnoughData,
    /// A non-finite value (NaN/inf) was found in inputs while `check_finite`
    /// is on (silx `asarray_chkfinite`).
    NonFinite,
    /// The curvature matrix `alpha0` is singular and cannot be inverted, so no
    /// covariance is available (silx `LinAlgError` from `inv(alpha0)`).
    SingularMatrix,
    /// A QUOTED constraint has equal min/max (`B == 0`), so the `sin/arcsin`
    /// reparametrisation is undefined (silx `raise ValueError("Invalid
    /// parameter limits")`).
    InvalidConstraint,
    /// A FACTOR/DELTA/SUM constraint references a parameter index outside the
    /// parameter vector.
    BadConstraintReference,
}

/// A per-parameter fit constraint for [`leastsq_constrained`].
///
/// Faithful to silx `silx.math.fit.leastsq` constraint codes (`CFREE`=0 …
/// `CIGNORED`=7). The constraint set has one entry per parameter; a parameter
/// is either *fitted* (Free/Positive/Quoted), *held* (Fixed/Ignored), or
/// *derived* from another parameter (Factor/Delta/Sum).
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum Constraint {
    /// `CFREE` — varied with no restriction.
    Free,
    /// `CPOSITIVE` — varied but kept positive (silx takes `abs(value)` rather
    /// than the commented-out square reparametrisation).
    Positive,
    /// `CQUOTED` — confined to `[min, max]` via silx's `A + B·sin(arcsin(...) +
    /// Δ)` reparametrisation (`A = (max+min)/2`, `B = (max−min)/2`). A start
    /// value outside `[min, max]` is held fixed at its start (silx warning
    /// path).
    Quoted {
        /// Lower bound (silx `min(constraints[i][1], constraints[i][2])`).
        min: f64,
        /// Upper bound (silx `max(...)`).
        max: f64,
    },
    /// `CFIXED` — held at its starting value; not fitted, gets 100 %
    /// uncertainty in the covariance.
    Fixed,
    /// `CFACTOR` — `value = factor · params[reference]` (silx `CFACTOR`).
    Factor {
        /// Index of the parameter this one is tied to.
        reference: usize,
        /// Multiplicative factor.
        factor: f64,
    },
    /// `CDELTA` — `value = delta + params[reference]` (silx `CDELTA`).
    Delta {
        /// Index of the parameter this one is tied to.
        reference: usize,
        /// Additive offset.
        delta: f64,
    },
    /// `CSUM` — `value = sum − params[reference]` (silx `CSUM`).
    Sum {
        /// Index of the parameter this one is tied to.
        reference: usize,
        /// Constant the pair sums to.
        sum: f64,
    },
    /// `CIGNORED` — set to 0 and stripped from the model call (silx `CIGNORED`;
    /// the model must accept the reduced parameter count).
    Ignored,
}

/// Apply the dependent-parameter relations to a full parameter vector (silx
/// `_get_parameters`). Free parameters pass through, Positive is `abs`, Quoted
/// passes through (its bounding is done at the step site), Factor/Delta/Sum are
/// recomputed from their reference, Ignored becomes 0. Two passes so that the
/// independent values are set before the dependent ones read them.
fn get_parameters(params: &[f64], constraints: &[Constraint]) -> Vec<f64> {
    let mut out: Vec<f64> = params
        .iter()
        .zip(constraints)
        .map(|(&p, c)| match c {
            Constraint::Positive => p.abs(),
            _ => p,
        })
        .collect();
    for (i, c) in constraints.iter().enumerate() {
        match *c {
            Constraint::Factor { reference, factor } => out[i] = factor * out[reference],
            Constraint::Delta { reference, delta } => out[i] = delta + out[reference],
            Constraint::Sum { reference, sum } => out[i] = sum - out[reference],
            Constraint::Ignored => out[i] = 0.0,
            _ => {}
        }
    }
    out
}

/// Propagate free-parameter sigmas back onto the full parameter vector through
/// the constraints (silx `_get_sigma_parameters`). `sigma0` holds one entry per
/// free parameter, in free-parameter order.
fn get_sigma_parameters(
    parameters: &[f64],
    sigma0: &[f64],
    constraints: &[Constraint],
) -> Vec<f64> {
    let mut sigma_par = vec![0.0_f64; parameters.len()];
    let mut n_free = 0usize;
    for (i, c) in constraints.iter().enumerate() {
        match *c {
            Constraint::Free | Constraint::Positive => {
                sigma_par[i] = sigma0[n_free];
                n_free += 1;
            }
            Constraint::Quoted { min, max } => {
                let pmax = min.max(max);
                let pmin = min.min(max);
                let b = 0.5 * (pmax - pmin);
                if b > 0.0 && parameters[i] < pmax && parameters[i] > pmin {
                    sigma_par[i] = (b * parameters[i].cos() * sigma0[n_free]).abs();
                    n_free += 1;
                } else {
                    sigma_par[i] = parameters[i];
                }
            }
            Constraint::Fixed => sigma_par[i] = parameters[i],
            _ => {}
        }
    }
    for (i, c) in constraints.iter().enumerate() {
        match *c {
            Constraint::Factor { reference, .. }
            | Constraint::Delta { reference, .. }
            | Constraint::Sum { reference, .. } => sigma_par[i] = sigma_par[reference],
            _ => {}
        }
    }
    sigma_par
}

/// Invert a square row-major matrix via Gauss-Jordan elimination with partial
/// pivoting. Returns `None` if the matrix is singular.
///
/// This stands in for numpy's `numpy.linalg.inv` used by silx `leastsq` (for
/// `inv(alpha)` in the LM step and `inv(alpha0)` for the covariance).
pub fn invert_matrix(m: &[Vec<f64>]) -> Option<Vec<Vec<f64>>> {
    let n = m.len();
    if n == 0 {
        return Some(Vec::new());
    }
    // Augment [ m | I ].
    let mut a: Vec<Vec<f64>> = Vec::with_capacity(n);
    for (i, row) in m.iter().enumerate() {
        if row.len() != n {
            return None;
        }
        let mut aug = row.clone();
        aug.extend((0..n).map(|j| if i == j { 1.0 } else { 0.0 }));
        a.push(aug);
    }
    for col in 0..n {
        // Partial pivot: largest magnitude in this column at/below the diagonal.
        let mut pivot = col;
        let mut best = a[col][col].abs();
        for (r, row) in a.iter().enumerate().skip(col + 1) {
            let v = row[col].abs();
            if v > best {
                best = v;
                pivot = r;
            }
        }
        if best == 0.0 {
            return None; // singular
        }
        a.swap(col, pivot);
        let pivot_val = a[col][col];
        for v in a[col].iter_mut() {
            *v /= pivot_val;
        }
        let pivot_row = a[col].clone();
        for (r, row) in a.iter_mut().enumerate() {
            if r == col {
                continue;
            }
            let factor = row[col];
            if factor != 0.0 {
                for (cell, &pv) in row.iter_mut().zip(pivot_row.iter()) {
                    *cell -= factor * pv;
                }
            }
        }
    }
    // Extract the right half.
    let inv = a
        .into_iter()
        .map(|row| row[n..].to_vec())
        .collect::<Vec<_>>();
    Some(inv)
}

/// Default `deltachi` (relative chi-square decrement, in percent) that stops
/// the LM iteration when an accepted step improves chi-square by less than
/// this. silx default is `0.001` (i.e. 0.1 %).
pub const DEFAULT_DELTACHI: f64 = 0.001;

/// Default maximum number of iterations (silx `max_iter=100`).
pub const DEFAULT_MAX_ITER: usize = 100;

/// Run an iterative Levenberg-Marquardt least-squares fit.
///
/// `model(x, params) -> y_hat` is evaluated over the whole `x` array. `p0` is
/// the initial parameter guess. `sigma` is the optional per-point uncertainty
/// used as weight (`weight = 1/sigma^2`); when `None`, every weight is 1, as in
/// silx (`sigma = numpy.ones(...)`).
///
/// Faithful to silx `leastsq` (unconstrained path): forward numerical
/// derivatives with step `delta[i] = (p[i] + (p[i]==0)) * sqrt(epsfcn)`,
/// `epsfcn = f64::EPSILON`, accept-if-chi-square-decreases with `flambda`
/// damping (start 0.001, `*10` on rejection up to 1000, `/10` on acceptance),
/// and the two-stop convergence test (`lastdeltachi < deltachi` or
/// `absdeltachi < sqrt(epsfcn)`), with the silx rule that the first iteration
/// always proceeds regardless of those limits.
pub fn leastsq<F>(
    model: F,
    xdata: &[f64],
    ydata: &[f64],
    p0: &[f64],
    sigma: Option<&[f64]>,
    max_iter: usize,
    deltachi: f64,
) -> Result<LeastSqResult, FitError>
where
    F: Fn(&[f64], &[f64]) -> Vec<f64>,
{
    if xdata.len() != ydata.len() {
        return Err(FitError::LengthMismatch);
    }
    let n_param = p0.len();
    if n_param == 0 {
        return Err(FitError::NoFreeParameters);
    }
    let m = ydata.len();
    if m < n_param {
        return Err(FitError::NotEnoughData);
    }
    // check_finite: silx asarray_chkfinite on xdata/ydata/sigma.
    if xdata.iter().chain(ydata.iter()).any(|v| !v.is_finite()) {
        return Err(FitError::NonFinite);
    }
    // weight0 = (1/sigma)^2 ; sigma==0 → divisor 1 (silx `sigma + (sigma==0)`).
    let weight0: Vec<f64> = match sigma {
        Some(s) => {
            if s.len() != m {
                return Err(FitError::LengthMismatch);
            }
            if s.iter().any(|v| !v.is_finite()) {
                return Err(FitError::NonFinite);
            }
            s.iter()
                .map(|&sv| {
                    let denom = if sv == 0.0 { 1.0 } else { sv };
                    let w = 1.0 / denom;
                    w * w
                })
                .collect()
        }
        None => vec![1.0; m],
    };

    let epsfcn = f64::EPSILON;
    let sqrt_epsfcn = epsfcn.sqrt();

    let mut fittedpar = p0.to_vec();
    let mut flambda = 0.001_f64;
    let mut iiter = max_iter as i64;
    let mut last_evaluation: Option<Vec<f64>> = None;
    let mut iteration_counter: usize = 0;
    let mut nfev: usize = 0;

    // Outputs of the most recent chisq_alpha_beta, captured for covariance.
    let mut chisq0: f64;
    let mut alpha0: Vec<Vec<f64>> = vec![vec![0.0; n_param]; n_param];

    loop {
        if iiter <= 0 {
            break;
        }
        iteration_counter += 1;

        // --- chisq_alpha_beta (unconstrained) ---
        // yfit at current parameters (reuse last_evaluation if available).
        let yfit0 = match &last_evaluation {
            Some(ev) => ev.clone(),
            None => {
                let ev = model(xdata, &fittedpar);
                nfev += 1;
                ev
            }
        };
        // delta[i] = (p[i] + (p[i]==0)) * sqrt(epsfcn)
        let delta: Vec<f64> = fittedpar
            .iter()
            .map(|&p| (p + if p == 0.0 { 1.0 } else { 0.0 }) * sqrt_epsfcn)
            .collect();
        // Forward numerical derivatives deriv[i][j] = (f(p+delta_i) - f0)/delta_i.
        let mut deriv: Vec<Vec<f64>> = Vec::with_capacity(n_param);
        for i in 0..n_param {
            let mut pwork = fittedpar.clone();
            pwork[i] = fittedpar[i] + delta[i];
            let f1 = model(xdata, &pwork);
            nfev += 1;
            let di = delta[i];
            let row: Vec<f64> = f1
                .iter()
                .zip(yfit0.iter())
                .map(|(&a, &b)| (a - b) / di)
                .collect();
            deriv.push(row);
        }
        // deltay = y - yfit ; help0 = weight * deltay
        let deltay: Vec<f64> = ydata
            .iter()
            .zip(yfit0.iter())
            .map(|(&y, &f)| y - f)
            .collect();
        let help0: Vec<f64> = weight0
            .iter()
            .zip(deltay.iter())
            .map(|(&w, &d)| w * d)
            .collect();
        // beta[i] = sum_j help0[j]*deriv[i][j]
        let mut beta = vec![0.0_f64; n_param];
        for i in 0..n_param {
            let mut s = 0.0;
            for j in 0..m {
                s += help0[j] * deriv[i][j];
            }
            beta[i] = s;
        }
        // alpha[i][k] = sum_j deriv[i][j]*weight[j]*deriv[k][j]
        let mut alpha = vec![vec![0.0_f64; n_param]; n_param];
        for i in 0..n_param {
            for k in 0..n_param {
                let mut s = 0.0;
                for j in 0..m {
                    s += deriv[i][j] * weight0[j] * deriv[k][j];
                }
                alpha[i][k] = s;
            }
        }
        // chisq = sum(help0 * deltay)
        chisq0 = help0.iter().zip(deltay.iter()).map(|(&h, &d)| h * d).sum();
        alpha0 = alpha.clone();

        // --- LM inner loop: pick a step that decreases chisq ---
        loop {
            // alpha' = alpha0 * (1 + flambda*I): only the diagonal is scaled.
            let mut alpha_lm = alpha0.clone();
            for (d, row) in alpha_lm.iter_mut().enumerate() {
                row[d] *= 1.0 + flambda;
            }
            let inv_alpha = match invert_matrix(&alpha_lm) {
                Some(inv) => inv,
                None => {
                    // Treat as a rejected step: damp harder.
                    flambda *= 10.0;
                    if flambda > 1000.0 {
                        iiter = 0;
                        break;
                    }
                    continue;
                }
            };
            // deltapar = beta · inv_alpha  (row-vector times matrix)
            // numpy: numpy.dot(beta, inv(alpha)) → sum_i beta[i]*inv[i][k]
            let mut deltapar = vec![0.0_f64; n_param];
            for (k, dp) in deltapar.iter_mut().enumerate() {
                let mut s = 0.0;
                for (i, &b) in beta.iter().enumerate() {
                    s += b * inv_alpha[i][k];
                }
                *dp = s;
            }
            let newpar: Vec<f64> = fittedpar
                .iter()
                .zip(deltapar.iter())
                .map(|(&p, &d)| p + d)
                .collect();
            let yfit = model(xdata, &newpar);
            nfev += 1;
            let chisq: f64 = weight0
                .iter()
                .zip(ydata.iter().zip(yfit.iter()))
                .map(|(&w, (&y, &f))| {
                    let r = y - f;
                    w * r * r
                })
                .sum();
            let absdeltachi = chisq0 - chisq;
            if absdeltachi < 0.0 {
                // Step worsened chi-square: reject, damp harder (silx flambda *= 10).
                flambda *= 10.0;
                if flambda > 1000.0 {
                    iiter = 0;
                    break;
                }
            } else {
                // Step improved chi-square: accept it.
                fittedpar = newpar;
                let lastdeltachi =
                    100.0 * (absdeltachi / (chisq + if chisq == 0.0 { 1.0 } else { 0.0 }));
                // silx convergence test: after the first iteration (which is
                // always allowed to proceed), stop when either the relative
                // chi-square decrement falls below `deltachi` OR the absolute
                // decrement falls below `sqrt(epsfcn)`. Both branches stop the
                // loop, so they are combined here.
                if iteration_counter >= 2 && (lastdeltachi < deltachi || absdeltachi < sqrt_epsfcn)
                {
                    iiter = 0;
                }
                // silx sets `chisq0 = chisq` here, but it is recomputed from
                // scratch at the top of every outer iteration via the
                // chisq_alpha_beta block, so persisting it has no effect.
                flambda /= 10.0;
                last_evaluation = Some(yfit);
                break;
            }
        }
        iiter -= 1;
    }

    // Covariance is inv(alpha0) (silx cov0).
    let covariance = invert_matrix(&alpha0).ok_or(FitError::SingularMatrix)?;
    let chisq_final = {
        // Recompute at the final parameters for a definite chisq value.
        let yfit = model(xdata, &fittedpar);
        nfev += 1;
        weight0
            .iter()
            .zip(ydata.iter().zip(yfit.iter()))
            .map(|(&w, (&y, &f))| {
                let r = y - f;
                w * r * r
            })
            .sum::<f64>()
    };
    let dof = m as i64 - n_param as i64;
    let reduced_chisq = if dof > 0 {
        Some(chisq_final / dof as f64)
    } else {
        None
    };

    // silx: constraints is None => uncertainties = sigma0 = sqrt(abs(diag(cov0))).
    let uncertainties: Vec<f64> = (0..n_param)
        .map(|i| covariance[i][i].abs().sqrt())
        .collect();

    Ok(LeastSqResult {
        parameters: fittedpar,
        covariance,
        uncertainties,
        chisq: chisq_final,
        reduced_chisq,
        niter: iteration_counter,
        nfev,
    })
}

/// Gather elements of `v` at the given indices (numpy `take`).
fn take(v: &[f64], indices: &[usize]) -> Vec<f64> {
    indices.iter().map(|&i| v[i]).collect()
}

/// Output of one constrained `chisq_alpha_beta` evaluation.
struct CabOut {
    /// `sum(weight * (y - yfit)^2)` at the current parameters.
    chisq: f64,
    /// Curvature matrix over the free parameters (`n_free x n_free`).
    alpha: Vec<Vec<f64>>,
    /// Gradient vector over the free parameters (`n_free`).
    beta: Vec<f64>,
    /// Number of free parameters this evaluation fitted.
    n_free: usize,
    /// Full-parameter index of each free parameter (silx `free_index`).
    free_index: Vec<usize>,
    /// Full-parameter indices passed to the model (silx `noigno`).
    noigno: Vec<usize>,
    /// Current values of the free parameters (silx `fitparam`).
    fitparam: Vec<f64>,
}

/// One constrained curvature/gradient evaluation (silx `chisq_alpha_beta` with
/// a non-None `constraints`). Builds the free-parameter set, the numerical
/// forward-difference Jacobian scaled by each parameter's `derivfactor`, and the
/// normal-equation matrices restricted to the free parameters.
#[allow(clippy::too_many_arguments)]
fn chisq_alpha_beta_constrained<F>(
    model: &F,
    parameters: &[f64],
    xdata: &[f64],
    ydata: &[f64],
    weight0: &[f64],
    constraints: &[Constraint],
    sqrt_epsfcn: f64,
    last_evaluation: Option<&[f64]>,
    nfev: &mut usize,
) -> CabOut
where
    F: Fn(&[f64], &[f64]) -> Vec<f64>,
{
    let m = ydata.len();

    // Classify parameters into the free set (silx: CFREE/CPOSITIVE always free,
    // CQUOTED free only when in bounds; others held/derived/ignored).
    let mut fitparam: Vec<f64> = Vec::new();
    let mut free_index: Vec<usize> = Vec::new();
    let mut noigno: Vec<usize> = Vec::new();
    let mut derivfactor: Vec<f64> = Vec::new();
    for (i, c) in constraints.iter().enumerate() {
        if !matches!(c, Constraint::Ignored) {
            noigno.push(i);
        }
        match *c {
            Constraint::Free => {
                fitparam.push(parameters[i]);
                derivfactor.push(1.0);
                free_index.push(i);
            }
            Constraint::Positive => {
                fitparam.push(parameters[i].abs());
                derivfactor.push(1.0);
                free_index.push(i);
            }
            Constraint::Quoted { min, max } => {
                let pmax = min.max(max);
                let pmin = min.min(max);
                if (pmax - pmin) > 0.0 && parameters[i] <= pmax && parameters[i] >= pmin {
                    let a = 0.5 * (pmax + pmin);
                    let b = 0.5 * (pmax - pmin);
                    fitparam.push(parameters[i]);
                    derivfactor.push(b * ((parameters[i] - a) / b).asin().cos());
                    free_index.push(i);
                }
                // Out of bounds: kept at its start value (not fitted).
            }
            _ => {}
        }
    }
    let n_free = fitparam.len();

    // delta[i] = (fitparam[i] + (fitparam[i]==0)) * sqrt(epsfcn).
    let delta: Vec<f64> = fitparam
        .iter()
        .map(|&p| (p + if p == 0.0 { 1.0 } else { 0.0 }) * sqrt_epsfcn)
        .collect();

    // pwork is the full parameter vector with the free values substituted in.
    let mut pwork = parameters.to_vec();
    for (i, &fi) in free_index.iter().enumerate() {
        pwork[fi] = fitparam[i];
    }

    // Base evaluation (silx f2 / yfit). We use one constraint-expanded base for
    // both the derivative reference and the residual, so the forward difference
    // is self-consistent; this matches silx's `model(*parameters)` exactly when
    // nothing is ignored and the start is in-domain (the practical case).
    let yfit: Vec<f64> = match last_evaluation {
        Some(ev) => ev.to_vec(),
        None => {
            let base_in = take(&get_parameters(&pwork, constraints), &noigno);
            let ev = model(xdata, &base_in);
            *nfev += 1;
            ev
        }
    };

    // Numerical forward derivatives for each free parameter, scaled by
    // derivfactor (CQUOTED `B·cos`, else 1).
    let mut deriv: Vec<Vec<f64>> = Vec::with_capacity(n_free);
    for i in 0..n_free {
        let fi = free_index[i];
        pwork[fi] = fitparam[i] + delta[i];
        let newpar = take(&get_parameters(&pwork, constraints), &noigno);
        let f1 = model(xdata, &newpar);
        *nfev += 1;
        let di = delta[i];
        let df = derivfactor[i];
        let row: Vec<f64> = f1
            .iter()
            .zip(yfit.iter())
            .map(|(&a, &b)| (a - b) / di * df)
            .collect();
        deriv.push(row);
        pwork[fi] = fitparam[i]; // restore
    }

    let deltay: Vec<f64> = ydata
        .iter()
        .zip(yfit.iter())
        .map(|(&y, &f)| y - f)
        .collect();
    let help0: Vec<f64> = weight0
        .iter()
        .zip(deltay.iter())
        .map(|(&w, &d)| w * d)
        .collect();
    let mut beta = vec![0.0_f64; n_free];
    for (i, b) in beta.iter_mut().enumerate() {
        let mut s = 0.0;
        for j in 0..m {
            s += help0[j] * deriv[i][j];
        }
        *b = s;
    }
    let mut alpha = vec![vec![0.0_f64; n_free]; n_free];
    for i in 0..n_free {
        for k in 0..n_free {
            let mut s = 0.0;
            for j in 0..m {
                s += deriv[i][j] * weight0[j] * deriv[k][j];
            }
            alpha[i][k] = s;
        }
    }
    let chisq = help0.iter().zip(deltay.iter()).map(|(&h, &d)| h * d).sum();

    CabOut {
        chisq,
        alpha,
        beta,
        n_free,
        free_index,
        noigno,
        fitparam,
    }
}

/// Run a constrained Levenberg-Marquardt least-squares fit (silx `leastsq` with
/// a non-`None` `constraints`).
///
/// `constraints` has one entry per parameter in `p0`. Free / Positive /
/// in-bounds Quoted parameters are varied; Fixed parameters are held at their
/// start (and receive 100 % uncertainty); Factor / Delta / Sum parameters are
/// derived from another parameter each step; Ignored parameters are set to 0 and
/// dropped from the model call (the model must accept the reduced count). When
/// every entry is [`Constraint::Free`] this matches an unconstrained fit except
/// that the covariance is recomputed at the final parameters (silx's second
/// `chisq_alpha_beta` pass).
///
/// The non-finite check, weighting, LM damping and convergence test match
/// [`leastsq`]; only the per-step parameter handling differs.
#[allow(clippy::too_many_arguments)]
pub fn leastsq_constrained<F>(
    model: F,
    xdata: &[f64],
    ydata: &[f64],
    p0: &[f64],
    constraints: &[Constraint],
    sigma: Option<&[f64]>,
    max_iter: usize,
    deltachi: f64,
) -> Result<LeastSqResult, FitError>
where
    F: Fn(&[f64], &[f64]) -> Vec<f64>,
{
    if xdata.len() != ydata.len() {
        return Err(FitError::LengthMismatch);
    }
    let n_param = p0.len();
    if n_param == 0 {
        return Err(FitError::NoFreeParameters);
    }
    if constraints.len() != n_param {
        return Err(FitError::BadConstraintReference);
    }
    let m = ydata.len();
    if m < 1 {
        return Err(FitError::NotEnoughData);
    }
    if xdata.iter().chain(ydata.iter()).any(|v| !v.is_finite()) {
        return Err(FitError::NonFinite);
    }
    // Validate Factor/Delta/Sum references and reject equal-bound Quoted.
    for c in constraints {
        match *c {
            Constraint::Factor { reference, .. }
            | Constraint::Delta { reference, .. }
            | Constraint::Sum { reference, .. } => {
                if reference >= n_param {
                    return Err(FitError::BadConstraintReference);
                }
            }
            Constraint::Quoted { min, max } => {
                if (min.max(max) - min.min(max)) == 0.0 {
                    return Err(FitError::InvalidConstraint);
                }
            }
            _ => {}
        }
    }

    let weight0: Vec<f64> = match sigma {
        Some(s) => {
            if s.len() != m {
                return Err(FitError::LengthMismatch);
            }
            if s.iter().any(|v| !v.is_finite()) {
                return Err(FitError::NonFinite);
            }
            s.iter()
                .map(|&sv| {
                    let denom = if sv == 0.0 { 1.0 } else { sv };
                    let w = 1.0 / denom;
                    w * w
                })
                .collect()
        }
        None => vec![1.0; m],
    };

    let epsfcn = f64::EPSILON;
    let sqrt_epsfcn = epsfcn.sqrt();

    // Count the initial free set so `alpha0` is sized even if no iteration runs
    // (max_iter == 0); also enforces silx's "No free parameters to fit".
    let n_free_initial = constraints
        .iter()
        .enumerate()
        .filter(|(i, c)| match **c {
            Constraint::Free | Constraint::Positive => true,
            Constraint::Quoted { min, max } => {
                let (pmax, pmin) = (min.max(max), min.min(max));
                (pmax - pmin) > 0.0 && p0[*i] <= pmax && p0[*i] >= pmin
            }
            _ => false,
        })
        .count();
    if n_free_initial == 0 {
        return Err(FitError::NoFreeParameters);
    }

    let mut fittedpar = p0.to_vec();
    let mut flambda = 0.001_f64;
    let mut iiter = max_iter as i64;
    let mut last_evaluation: Option<Vec<f64>> = None;
    let mut iteration_counter: usize = 0;
    let mut nfev: usize = 0;

    let mut alpha0: Vec<Vec<f64>> = vec![vec![0.0; n_free_initial]; n_free_initial];
    let mut n_free_final = n_free_initial;

    loop {
        if iiter <= 0 {
            break;
        }
        iteration_counter += 1;

        let cab = chisq_alpha_beta_constrained(
            &model,
            &fittedpar,
            xdata,
            ydata,
            &weight0,
            constraints,
            sqrt_epsfcn,
            last_evaluation.as_deref(),
            &mut nfev,
        );
        let chisq0 = cab.chisq;
        alpha0 = cab.alpha.clone();
        n_free_final = cab.n_free;
        let beta = &cab.beta;
        let free_index = &cab.free_index;
        let noigno = &cab.noigno;
        let fitparam = &cab.fitparam;
        if cab.n_free == 0 {
            return Err(FitError::NoFreeParameters);
        }

        loop {
            let mut alpha_lm = alpha0.clone();
            for (d, row) in alpha_lm.iter_mut().enumerate() {
                row[d] *= 1.0 + flambda;
            }
            let inv_alpha = match invert_matrix(&alpha_lm) {
                Some(inv) => inv,
                None => {
                    flambda *= 10.0;
                    if flambda > 1000.0 {
                        iiter = 0;
                        break;
                    }
                    continue;
                }
            };
            // deltapar = beta · inv_alpha (free-parameter space).
            let mut deltapar = vec![0.0_f64; cab.n_free];
            for (k, dp) in deltapar.iter_mut().enumerate() {
                let mut s = 0.0;
                for (i, &b) in beta.iter().enumerate() {
                    s += b * inv_alpha[i][k];
                }
                *dp = s;
            }
            // Rebuild the full parameter vector: Fixed/derived entries come from
            // the original p0 template, free entries take the LM step (Quoted via
            // the sin reparametrisation), then get_parameters applies the
            // dependent relations.
            let mut newpar = p0.to_vec();
            for (i, &fi) in free_index.iter().enumerate() {
                let pv = match constraints[fi] {
                    Constraint::Quoted { min, max } => {
                        let pmax = min.max(max);
                        let pmin = min.min(max);
                        let a = 0.5 * (pmax + pmin);
                        let b = 0.5 * (pmax - pmin);
                        a + b * (((fitparam[i] - a) / b).asin() + deltapar[i]).sin()
                    }
                    // Free and Positive both step additively; Positive's abs is
                    // applied by get_parameters below.
                    _ => fitparam[i] + deltapar[i],
                };
                newpar[fi] = pv;
            }
            let newpar = get_parameters(&newpar, constraints);
            let workpar = take(&newpar, noigno);
            let yfit = model(xdata, &workpar);
            nfev += 1;
            let chisq: f64 = weight0
                .iter()
                .zip(ydata.iter().zip(yfit.iter()))
                .map(|(&w, (&y, &f))| {
                    let r = y - f;
                    w * r * r
                })
                .sum();
            let absdeltachi = chisq0 - chisq;
            if absdeltachi < 0.0 {
                flambda *= 10.0;
                if flambda > 1000.0 {
                    iiter = 0;
                    break;
                }
            } else {
                fittedpar = newpar;
                let lastdeltachi =
                    100.0 * (absdeltachi / (chisq + if chisq == 0.0 { 1.0 } else { 0.0 }));
                if iteration_counter >= 2 && (lastdeltachi < deltachi || absdeltachi < sqrt_epsfcn)
                {
                    iiter = 0;
                }
                flambda /= 10.0;
                last_evaluation = Some(yfit);
                break;
            }
        }
        iiter -= 1;
    }

    // cov0 = inv(alpha0): the free-space covariance (silx).
    let cov0 = invert_matrix(&alpha0).ok_or(FitError::SingularMatrix)?;

    // Second pass: every non-fixed/ignored parameter becomes Free, recompute the
    // curvature at the final parameters and invert for the reported covariance;
    // fixed/ignored parameters get a zero row/col and a 100 % diagonal.
    let new_constraints: Vec<Constraint> = constraints
        .iter()
        .map(|c| match c {
            Constraint::Fixed | Constraint::Ignored => *c,
            _ => Constraint::Free,
        })
        .collect();
    let cab2 = chisq_alpha_beta_constrained(
        &model,
        &fittedpar,
        xdata,
        ydata,
        &weight0,
        &new_constraints,
        sqrt_epsfcn,
        last_evaluation.as_deref(),
        &mut nfev,
    );
    let mut covariance = vec![vec![0.0_f64; n_param]; n_param];
    if let Some(cov_free) = invert_matrix(&cab2.alpha) {
        for (r, &pr) in cab2.free_index.iter().enumerate() {
            for (cc, &pc) in cab2.free_index.iter().enumerate() {
                covariance[pr][pc] = cov_free[r][cc];
            }
        }
    }
    for (idx, c) in constraints.iter().enumerate() {
        if matches!(c, Constraint::Fixed | Constraint::Ignored) {
            covariance[idx][idx] = fittedpar[idx] * fittedpar[idx];
        }
    }

    // Uncertainties: sigma0 = sqrt(abs(diag(cov0))) over the free parameters,
    // propagated through the constraints (silx _get_sigma_parameters).
    let sigma0: Vec<f64> = (0..n_free_final).map(|i| cov0[i][i].abs().sqrt()).collect();
    let uncertainties = get_sigma_parameters(&fittedpar, &sigma0, constraints);

    // Final chi-square at the converged parameters.
    let workpar = take(&get_parameters(&fittedpar, constraints), &cab2.noigno);
    let yfit_final = model(xdata, &workpar);
    nfev += 1;
    let chisq_final: f64 = weight0
        .iter()
        .zip(ydata.iter().zip(yfit_final.iter()))
        .map(|(&w, (&y, &f))| {
            let r = y - f;
            w * r * r
        })
        .sum();
    let dof = m as i64 - n_free_final as i64;
    let reduced_chisq = if dof > 0 {
        Some(chisq_final / dof as f64)
    } else {
        None
    };

    Ok(LeastSqResult {
        parameters: fittedpar,
        covariance,
        uncertainties,
        chisq: chisq_final,
        reduced_chisq,
        niter: iteration_counter,
        nfev,
    })
}

// ---------------------------------------------------------------------------
// Peak models (CPU). Each evaluates a single peak + flat background:
// y(x) = peak(x; params...) + background.
//
// The peak formulas are ported byte-for-byte from
// `silx/math/fit/functions/src/funs.c` (single-peak case of the sum_* loops).
// A trailing `background` parameter is appended (constant offset) so that a
// model is fully described by one parameter vector for `leastsq`.
// ---------------------------------------------------------------------------

/// Evaluate a Gaussian peak (height parameterisation) plus flat background.
///
/// `params = [height, centroid, fwhm, background]`. Mirrors C `sum_gauss`:
/// `sigma = fwhm / (2*sqrt(2*LOG2))`, `y = height*exp(-0.5*((x-c)/sigma)^2)`,
/// with the C guard `(x-c)/sigma <= 20` skipping far-tail terms.
pub fn gaussian_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, centroid, fwhm, bg) = (params[0], params[1], params[2], params[3]);
    let sigma = fwhm / fwhm_to_sigma_factor();
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if sigma != 0.0 {
                let dhelp = (xi - centroid) / sigma;
                if dhelp <= 20.0 {
                    y += height * (-0.5 * dhelp * dhelp).exp();
                }
            }
            y
        })
        .collect()
}

/// Evaluate a Gaussian peak (area parameterisation) plus flat background.
///
/// `params = [area, centroid, fwhm, background]`. Mirrors C `sum_agauss`:
/// `sigma = fwhm/(2*sqrt(2*LOG2))`, `height = area/(sigma*sqrt(2*pi))`,
/// with the C guard `(x-c)/sigma <= 35`.
pub fn gaussian_area_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (area, centroid, fwhm, bg) = (params[0], params[1], params[2], params[3]);
    let sigma = fwhm / fwhm_to_sigma_factor();
    let sqrt2pi = (2.0 * std::f64::consts::PI).sqrt();
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if sigma != 0.0 {
                let height = area / (sigma * sqrt2pi);
                let dhelp = (xi - centroid) / sigma;
                if dhelp <= 35.0 {
                    y += height * (-0.5 * dhelp * dhelp).exp();
                }
            }
            y
        })
        .collect()
}

/// Evaluate a Lorentzian peak (height parameterisation) plus flat background.
///
/// `params = [height, centroid, fwhm, background]`. Mirrors C `sum_lorentz`:
/// `dhelp = (x-c)/(0.5*fwhm)`, `y = height/(1 + dhelp^2)`.
pub fn lorentzian_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, centroid, fwhm, bg) = (params[0], params[1], params[2], params[3]);
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if fwhm != 0.0 {
                let dhelp = (xi - centroid) / (0.5 * fwhm);
                y += height / (1.0 + dhelp * dhelp);
            }
            y
        })
        .collect()
}

/// Evaluate a pseudo-Voigt peak (height parameterisation) plus flat background.
///
/// `params = [height, centroid, fwhm, eta, background]`. Mirrors C
/// `sum_pvoigt`: `PV = eta*L + (1-eta)*G` where `L = height/(1+((x-c)/(0.5*fwhm))^2)`
/// and `G = height*exp(-0.5*((x-c)/sigma)^2)` with `sigma = fwhm/(2*sqrt(2*LOG2))`,
/// C guard `(x-c)/sigma <= 35` on the Gaussian term.
pub fn pseudo_voigt_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, centroid, fwhm, eta, bg) = (params[0], params[1], params[2], params[3], params[4]);
    let sigma = fwhm / fwhm_to_sigma_factor();
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if fwhm != 0.0 {
                // Lorentzian term.
                let dl = (xi - centroid) / (0.5 * fwhm);
                y += eta * height / (1.0 + dl * dl);
            }
            if sigma != 0.0 {
                // Gaussian term.
                let dg = (xi - centroid) / sigma;
                if dg <= 35.0 {
                    y += (1.0 - eta) * height * (-0.5 * dg * dg).exp();
                }
            }
            y
        })
        .collect()
}

// ---------------------------------------------------------------------------
// Step / slit models (non-peak fit theories).
//
// Ported from `silx/math/fit/functions/src/funs.c` (`sum_stepdown`,
// `sum_stepup`, `sum_slit`) and the pure-Python `atan_stepup`
// (`silx/math/fit/functions/functions.pyx`). As with the peak models above, a
// trailing constant `background` parameter is appended so each model is a
// complete `leastsq` target; silx keeps the baseline in a separate background
// theory.
//
// `erf`/`erfc` are not in Rust's std. They are approximated with
// Abramowitz & Stegun 7.1.26 (`|error| <= 1.5e-7`). silx calls the C library
// `erf`/`erfc` at full double precision; the difference is far below fit noise
// but is documented here rather than hidden.
// ---------------------------------------------------------------------------

/// Gaussian error function via Abramowitz & Stegun 7.1.26 (`|error| <= 1.5e-7`).
/// Exact at `x == 0` (returns `0`) and odd-symmetric, so the step models hit
/// their half-height exactly at the centre.
fn erf(x: f64) -> f64 {
    if x == 0.0 {
        return 0.0;
    }
    let sign = if x < 0.0 { -1.0 } else { 1.0 };
    let x = x.abs();
    // A&S 7.1.26 coefficients.
    const A1: f64 = 0.254829592;
    const A2: f64 = -0.284496736;
    const A3: f64 = 1.421413741;
    const A4: f64 = -1.453152027;
    const A5: f64 = 1.061405429;
    const P: f64 = 0.3275911;
    let t = 1.0 / (1.0 + P * x);
    let poly = ((((A5 * t + A4) * t + A3) * t + A2) * t + A1) * t;
    sign * (1.0 - poly * (-x * x).exp())
}

/// Complementary error function, `1 - erf(x)`.
fn erfc(x: f64) -> f64 {
    1.0 - erf(x)
}

/// The C `sum_step*` edge scale: `denom = fwhm * sqrt(2) / (2*sqrt(2*LOG2))`,
/// i.e. `sigma * sqrt(2)`. The erf argument is `(x - centre) / denom`.
fn step_denom(fwhm: f64) -> f64 {
    fwhm * std::f64::consts::SQRT_2 / fwhm_to_sigma_factor()
}

/// Evaluate a step-down (descending error-function edge) plus flat background.
///
/// `params = [height, centroid, fwhm, background]`. Mirrors C `sum_stepdown`:
/// `y = background + height * 0.5 * erfc((x - centroid) / denom)` where `denom`
/// is `step_denom`. `fwhm` is the full-width at half maximum of the step's
/// derivative (its sharpness).
pub fn stepdown_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, centroid, fwhm, bg) = (params[0], params[1], params[2], params[3]);
    let denom = step_denom(fwhm);
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if denom != 0.0 {
                y += height * 0.5 * erfc((xi - centroid) / denom);
            }
            y
        })
        .collect()
}

/// Evaluate a step-up (ascending error-function edge) plus flat background.
///
/// `params = [height, centroid, fwhm, background]`. Mirrors C `sum_stepup`:
/// `y = background + height * 0.5 * (1 + erf((x - centroid) / denom))`.
pub fn stepup_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, centroid, fwhm, bg) = (params[0], params[1], params[2], params[3]);
    let denom = step_denom(fwhm);
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if denom != 0.0 {
                y += height * 0.5 * (1.0 + erf((xi - centroid) / denom));
            }
            y
        })
        .collect()
}

/// Evaluate a slit (a rising then falling pair of edges) plus flat background.
///
/// `params = [height, position, fwhm, beamfwhm, background]`. Mirrors C
/// `sum_slit`: with `c1 = position - 0.5*fwhm`, `c2 = position + 0.5*fwhm` and
/// `denom = step_denom(beamfwhm)`,
/// `y = background + height * 0.25 * (1 + erf((x-c1)/denom)) * erfc((x-c2)/denom)`.
/// `fwhm` is the slit width; `beamfwhm` is the sharpness of its two edges.
pub fn slit_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, position, fwhm, beamfwhm, bg) =
        (params[0], params[1], params[2], params[3], params[4]);
    let denom = step_denom(beamfwhm);
    let c1 = position - 0.5 * fwhm;
    let c2 = position + 0.5 * fwhm;
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if denom != 0.0 {
                y += height * 0.25 * (1.0 + erf((xi - c1) / denom)) * erfc((xi - c2) / denom);
            }
            y
        })
        .collect()
}

/// Evaluate an arctan step-up plus flat background.
///
/// `params = [height, position, width, background]`. Mirrors Python
/// `atan_stepup`: `y = background + height * (0.5 + atan((x - position)/width)/pi)`.
/// A lower `width` yields a sharper step.
pub fn atan_stepup_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let (height, position, width, bg) = (params[0], params[1], params[2], params[3]);
    x.iter()
        .map(|&xi| {
            let mut y = bg;
            if width != 0.0 {
                y += height * (0.5 + ((xi - position) / width).atan() / std::f64::consts::PI);
            }
            y
        })
        .collect()
}

// ---------------------------------------------------------------------------
// Initial-parameter estimators.
//
// silx `estimate_height_position_fwhm` runs a peak search + strip background +
// a 4-iteration constrained micro-fit. The strip/snip background estimator and
// multi-peak search are DEFERRED; we port the single-peak analytical seed:
// background = min(y), height = max(y) - background, centroid = x[argmax],
// fwhm from the half-maximum crossing (the same shape silx ends up with for a
// single dominant peak). Area/eta conversions follow `estimate_agauss` /
// `estimate_pvoigt`.
// ---------------------------------------------------------------------------

/// Analytical single-peak seed: `(height, centroid, fwhm, background)`.
///
/// `background = min(y)`; `height = max(y) - background`; `centroid` is the
/// `x` at the maximum; `fwhm` is the width between the outermost half-maximum
/// crossings around the peak. Returns `None` if there are fewer than 3 points
/// or lengths differ.
pub fn estimate_height_position_fwhm(x: &[f64], y: &[f64]) -> Option<(f64, f64, f64, f64)> {
    if x.len() != y.len() || x.len() < 3 {
        return None;
    }
    let bg = y.iter().copied().fold(f64::INFINITY, f64::min);
    let mut max_y = f64::NEG_INFINITY;
    let mut max_idx = 0;
    for (i, &yi) in y.iter().enumerate() {
        if yi > max_y {
            max_y = yi;
            max_idx = i;
        }
    }
    let height = max_y - bg;
    let centroid = x[max_idx];
    let half_max = bg + height / 2.0;
    let mut left = max_idx;
    while left > 0 && y[left] > half_max {
        left -= 1;
    }
    let mut right = max_idx;
    while right < y.len() - 1 && y[right] > half_max {
        right += 1;
    }
    let fwhm = if right > left {
        x[right] - x[left]
    } else {
        (x[x.len() - 1] - x[0]).abs() / 4.0
    };
    let fwhm = if fwhm > 0.0 {
        fwhm
    } else {
        (x[x.len() - 1] - x[0]).abs() / 4.0
    };
    Some((height, centroid, fwhm, bg))
}

/// Seed for [`gaussian_model`]: `[height, centroid, fwhm, background]`.
pub fn estimate_gaussian(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    let (h, c, f, bg) = estimate_height_position_fwhm(x, y)?;
    Some(vec![h, c, f, bg])
}

/// Seed for [`gaussian_area_model`]: `[area, centroid, fwhm, background]`.
///
/// Area conversion mirrors silx `estimate_agauss`:
/// `area = sqrt(2*pi) * height * fwhm / (2*sqrt(2*ln2))`.
pub fn estimate_gaussian_area(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    let (h, c, f, bg) = estimate_height_position_fwhm(x, y)?;
    let area = (2.0 * std::f64::consts::PI).sqrt() * h * f / fwhm_to_sigma_factor();
    Some(vec![area, c, f, bg])
}

/// Seed for [`lorentzian_model`]: `[height, centroid, fwhm, background]`.
///
/// Same height/position/fwhm seed as Gaussian (silx `estimate_lorentz` reuses
/// `estimate_height_position_fwhm` without converting height).
pub fn estimate_lorentzian(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    let (h, c, f, bg) = estimate_height_position_fwhm(x, y)?;
    Some(vec![h, c, f, bg])
}

/// Seed for [`pseudo_voigt_model`]: `[height, centroid, fwhm, eta, background]`.
///
/// Eta seeds to 0.5, mirroring silx `estimate_pvoigt` (`newpar.append(0.5)`).
pub fn estimate_pseudo_voigt(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    let (h, c, f, bg) = estimate_height_position_fwhm(x, y)?;
    Some(vec![h, c, f, 0.5, bg])
}

/// `numpy.convolve(y, kernel, mode="valid")`: the kernel is applied reversed,
/// and the output length is `y.len() - kernel.len() + 1` (empty when `y` is
/// shorter than `kernel`).
fn convolve_valid(y: &[f64], kernel: &[f64]) -> Vec<f64> {
    let (n, m) = (y.len(), kernel.len());
    if n < m || m == 0 {
        return Vec::new();
    }
    (0..=n - m)
        .map(|k| (0..m).map(|j| y[k + j] * kernel[m - 1 - j]).sum())
        .collect()
}

/// Shared step-edge seed used by [`estimate_stepup`] / [`estimate_stepdown`].
///
/// silx convolves `y` with an edge-detecting kernel, then takes the dominant
/// peak of that derivative as the step centre and width; the height is the data
/// amplitude `max(y) - min(y)`. The derivative is *not* rescaled (silx's
/// `y_deriv *= max(y)/max(y_deriv)` is a uniform positive scale that leaves the
/// argmax and half-maximum crossings — hence centre and fwhm — unchanged).
/// Multi-step search is out of scope: the single dominant edge is used, matching
/// the peak estimators. The appended `background = min(y)`.
fn estimate_step(x: &[f64], y: &[f64], kernel: &[f64]) -> Option<Vec<f64>> {
    if x.len() != y.len() || x.len() < 3 {
        return None;
    }
    let bg = y.iter().copied().fold(f64::INFINITY, f64::min);
    let max_y = y.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    let data_amplitude = max_y - bg;

    let cutoff = kernel.len() / 2;
    let y_deriv = convolve_valid(y, kernel);
    let (center, fwhm) = if y_deriv.len() >= 3 && x.len() > 2 * cutoff {
        let x_slice = &x[cutoff..x.len() - cutoff];
        match estimate_height_position_fwhm(x_slice, &y_deriv) {
            Some((_h, c, f, _b)) => (c, f),
            None => step_fallback(x),
        }
    } else {
        step_fallback(x)
    };
    Some(vec![data_amplitude, center, fwhm, bg])
}

/// silx no-peak fallback: centre at the middle of `x`, `fwhm = FwhmPoints * dx`
/// with the silx default `FwhmPoints = 8`.
fn step_fallback(x: &[f64]) -> (f64, f64) {
    let center = x[x.len() / 2];
    let dx = if x.len() > 1 { x[1] - x[0] } else { 1.0 };
    (center, 8.0 * dx)
}

/// Seed for [`stepup_model`]: `[height, centroid, fwhm, background]`.
///
/// Mirrors silx `estimate_stepup` (edge kernel `[0.25, 0.75, 0, -0.75, -0.25]`).
pub fn estimate_stepup(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    estimate_step(x, y, &[0.25, 0.75, 0.0, -0.75, -0.25])
}

/// Seed for [`stepdown_model`]: `[height, centroid, fwhm, background]`.
///
/// Mirrors silx `estimate_stepdown` (edge kernel `[-0.25, -0.75, 0, 0.75, 0.25]`).
pub fn estimate_stepdown(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    estimate_step(x, y, &[-0.25, -0.75, 0.0, 0.75, 0.25])
}

/// Seed for [`atan_stepup_model`]: `[height, position, width, background]`.
///
/// silx uses `estimate_stepup` for the arctan step (the same edge detection;
/// the step fwhm seeds the arctan `width`).
pub fn estimate_atan_stepup(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    estimate_stepup(x, y)
}

/// Seed for [`slit_model`]: `[height, position, fwhm, beamfwhm, background]`.
///
/// Mirrors silx `estimate_slit`: seed the up- and down-edges to size the beam
/// sharpness, then take the slit centre/width from the half-maximum crossings
/// of the background-subtracted data. silx subtracts a *strip* background; that
/// filter is a separate theory (not yet ported), so the baseline here is
/// `min(y)` — the same simplification the peak estimators use. `beamfwhm` is
/// seeded as the average of the up/down edge FWHMs (silx's docstring intent;
/// its code has an index typo that reads the down-step centre instead), then
/// clamped to silx's `[range*3/n, edge_distance/10]` bounds.
pub fn estimate_slit(x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
    let up = estimate_stepup(x, y)?; // [h, center_up, fwhm_up, bg]
    let down = estimate_stepdown(x, y)?; // [h, center_down, fwhm_down, bg]
    let (center_up, fwhm_up) = (up[1], up[2]);
    let (center_down, fwhm_down) = (down[1], down[2]);
    let edge_distance = (center_down - center_up).abs();

    let bg = y.iter().copied().fold(f64::INFINITY, f64::min);
    let y_minus_bg: Vec<f64> = y.iter().map(|&yi| yi - bg).collect();
    let height = y_minus_bg.iter().copied().fold(f64::NEG_INFINITY, f64::max);

    // Slit centre/width from the half-maximum crossings of (y - background).
    let threshold = 0.5 * height;
    let first = y_minus_bg.iter().position(|&v| v >= threshold)?;
    let last = y_minus_bg.iter().rposition(|&v| v >= threshold)?;
    let position = (x[first] + x[last]) / 2.0;
    let fwhm = x[last] - x[first];

    // Beam sharpness: average of the edge FWHMs, clamped to silx's bounds.
    let mut beamfwhm = 0.5 * (fwhm_up + fwhm_down);
    beamfwhm = beamfwhm.min(edge_distance / 10.0);
    let xmin = x.iter().copied().fold(f64::INFINITY, f64::min);
    let xmax = x.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    beamfwhm = beamfwhm.max((xmax - xmin) * 3.0 / x.len() as f64);

    Some(vec![height, position, fwhm, beamfwhm, bg])
}

// ---------------------------------------------------------------------------
// Iterative fit models exposed through the FitFunction trait, and fit range.
// ---------------------------------------------------------------------------

/// Which fit model (peak, step, or slit) an [`IterativeFit`] fits.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum PeakModel {
    /// Gaussian, height parameterisation: `[height, centroid, fwhm, bg]`.
    Gaussian,
    /// Gaussian, area parameterisation: `[area, centroid, fwhm, bg]`.
    GaussianArea,
    /// Lorentzian, height parameterisation: `[height, centroid, fwhm, bg]`.
    Lorentzian,
    /// Pseudo-Voigt: `[height, centroid, fwhm, eta, bg]`.
    PseudoVoigt,
    /// Step down (descending erf edge): `[height, centroid, fwhm, bg]`.
    StepDown,
    /// Step up (ascending erf edge): `[height, centroid, fwhm, bg]`.
    StepUp,
    /// Slit (rising then falling edges): `[height, position, fwhm, beamfwhm, bg]`.
    Slit,
    /// Arctan step up: `[height, position, width, bg]`.
    AtanStepUp,
}

impl PeakModel {
    /// Display name for this model.
    pub fn name(self) -> &'static str {
        match self {
            PeakModel::Gaussian => "Gaussian",
            PeakModel::GaussianArea => "Gaussian (Area)",
            PeakModel::Lorentzian => "Lorentzian",
            PeakModel::PseudoVoigt => "Pseudo-Voigt",
            PeakModel::StepDown => "Step Down",
            PeakModel::StepUp => "Step Up",
            PeakModel::Slit => "Slit",
            PeakModel::AtanStepUp => "Arctan Step Up",
        }
    }

    /// Parameter names for this model, in parameter-vector order.
    pub fn param_names(self) -> Vec<String> {
        let owned = |s: &str| s.to_string();
        match self {
            PeakModel::Gaussian => vec![
                owned("Height"),
                owned("Center"),
                owned("FWHM"),
                owned("Background"),
            ],
            PeakModel::GaussianArea => vec![
                owned("Area"),
                owned("Center"),
                owned("FWHM"),
                owned("Background"),
            ],
            PeakModel::Lorentzian => vec![
                owned("Height"),
                owned("Center"),
                owned("FWHM"),
                owned("Background"),
            ],
            PeakModel::PseudoVoigt => vec![
                owned("Height"),
                owned("Center"),
                owned("FWHM"),
                owned("Eta"),
                owned("Background"),
            ],
            PeakModel::StepDown | PeakModel::StepUp => vec![
                owned("Height"),
                owned("Center"),
                owned("FWHM"),
                owned("Background"),
            ],
            PeakModel::Slit => vec![
                owned("Height"),
                owned("Center"),
                owned("FWHM"),
                owned("BeamFWHM"),
                owned("Background"),
            ],
            PeakModel::AtanStepUp => vec![
                owned("Height"),
                owned("Center"),
                owned("Width"),
                owned("Background"),
            ],
        }
    }

    /// Evaluate this model over `x` with the given parameter vector.
    pub fn eval(self, x: &[f64], params: &[f64]) -> Vec<f64> {
        match self {
            PeakModel::Gaussian => gaussian_model(x, params),
            PeakModel::GaussianArea => gaussian_area_model(x, params),
            PeakModel::Lorentzian => lorentzian_model(x, params),
            PeakModel::PseudoVoigt => pseudo_voigt_model(x, params),
            PeakModel::StepDown => stepdown_model(x, params),
            PeakModel::StepUp => stepup_model(x, params),
            PeakModel::Slit => slit_model(x, params),
            PeakModel::AtanStepUp => atan_stepup_model(x, params),
        }
    }

    /// Estimate an initial parameter vector for this model from the data.
    pub fn estimate(self, x: &[f64], y: &[f64]) -> Option<Vec<f64>> {
        match self {
            PeakModel::Gaussian => estimate_gaussian(x, y),
            PeakModel::GaussianArea => estimate_gaussian_area(x, y),
            PeakModel::Lorentzian => estimate_lorentzian(x, y),
            PeakModel::PseudoVoigt => estimate_pseudo_voigt(x, y),
            PeakModel::StepDown => estimate_stepdown(x, y),
            PeakModel::StepUp => estimate_stepup(x, y),
            PeakModel::Slit => estimate_slit(x, y),
            PeakModel::AtanStepUp => estimate_atan_stepup(x, y),
        }
    }
}

/// Outcome of an iterative peak fit: the [`FitResult`] plus the solver
/// diagnostics needed for a results table (errors + reduced chi-square).
#[derive(Debug, Clone)]
pub struct IterativeFitResult {
    /// The fitted curve and parameters (compatible with the simple fitters).
    pub fit: FitResult,
    /// Full solver output (covariance, chi-square, iteration counts).
    pub solver: LeastSqResult,
}

impl IterativeFitResult {
    /// Per-parameter standard errors (`sqrt(diag(covariance))`).
    pub fn std_errors(&self) -> Vec<f64> {
        self.solver.std_errors()
    }

    /// Reduced chi-square, if degrees of freedom were positive.
    pub fn reduced_chisq(&self) -> Option<f64> {
        self.solver.reduced_chisq
    }
}

/// An iterative (Levenberg-Marquardt) peak fitter for one [`PeakModel`].
///
/// Estimates initial parameters with [`PeakModel::estimate`], then refines them
/// with [`leastsq`]. The [`FitFunction`] impl returns the refined [`FitResult`];
/// use [`IterativeFit::fit_full`] to also obtain the covariance / chi-square.
pub struct IterativeFit {
    /// The peak model fitted by this instance.
    pub model: PeakModel,
    /// Maximum LM iterations (defaults to [`DEFAULT_MAX_ITER`]).
    pub max_iter: usize,
    /// Relative chi-square stop threshold (defaults to [`DEFAULT_DELTACHI`]).
    pub deltachi: f64,
}

impl IterativeFit {
    /// Create an iterative fitter for `model` with silx default iteration
    /// controls.
    pub fn new(model: PeakModel) -> Self {
        Self {
            model,
            max_iter: DEFAULT_MAX_ITER,
            deltachi: DEFAULT_DELTACHI,
        }
    }

    /// Fit and return the full solver diagnostics (covariance, chi-square).
    pub fn fit_full(&self, x: &[f64], y: &[f64]) -> Option<IterativeFitResult> {
        let p0 = self.model.estimate(x, y)?;
        let model = self.model;
        let solver = leastsq(
            |xx, pp| model.eval(xx, pp),
            x,
            y,
            &p0,
            None,
            self.max_iter,
            self.deltachi,
        )
        .ok()?;
        let y_fit = self.model.eval(x, &solver.parameters);
        let fit = FitResult {
            y_fit,
            parameters: solver.parameters.clone(),
            param_names: self.model.param_names(),
        };
        Some(IterativeFitResult { fit, solver })
    }
}

impl FitFunction for IterativeFit {
    fn name(&self) -> &str {
        self.model.name()
    }

    fn fit(&self, x: &[f64], y: &[f64]) -> Option<FitResult> {
        self.fit_full(x, y).map(|r| r.fit)
    }
}

/// Fit `model` to only the data points whose `x` falls within `[xmin, xmax]`
/// (inclusive), mirroring silx `FitWidget` xmin/xmax range restriction.
///
/// Points outside the range are dropped before fitting, so they cannot
/// influence the fitted parameters. `xmin`/`xmax` may be given in any order.
/// Returns `None` if fewer than 3 points remain in range.
pub fn fit_in_range(
    xs: &[f64],
    ys: &[f64],
    xmin: f64,
    xmax: f64,
    model: &IterativeFit,
) -> Option<IterativeFitResult> {
    if xs.len() != ys.len() {
        return None;
    }
    let (lo, hi) = if xmin <= xmax {
        (xmin, xmax)
    } else {
        (xmax, xmin)
    };
    let mut xr = Vec::new();
    let mut yr = Vec::new();
    for (&xi, &yi) in xs.iter().zip(ys.iter()) {
        if xi >= lo && xi <= hi {
            xr.push(xi);
            yr.push(yi);
        }
    }
    if xr.len() < 3 {
        return None;
    }
    model.fit_full(&xr, &yr)
}

// ---------------------------------------------------------------------------
// Multi-peak (simultaneous) Gaussian fitting.
//
// Ports silx `fittheories.estimate_height_position_fwhm` + `functions.sum_gauss`
// (C `sum_gauss`): locate N peaks with `peak_search`, seed one (height, centre,
// FWHM) triple per peak, and fit them all at once with the constrained solver.
// This composes `peak_search` and `leastsq_constrained` into the multi-peak fit
// that `IterativeFit` (single peak) does not provide.
// ---------------------------------------------------------------------------

/// Sum of `N` Gaussians (silx `functions.sum_gauss` / C `sum_gauss`).
///
/// `params` is a flat vector of `(height, centroid, fwhm)` triples — one per
/// peak, with **no** background term. `y(x) = Σ_k height_k · exp(-0.5·((x −
/// centroid_k)/sigma_k)²)`, `sigma = fwhm / (2·√(2·ln2))`, with the C far-tail
/// guard skipping terms where `(x − centroid)/sigma > 20`. A trailing partial
/// triple (length not a multiple of 3) is ignored, as silx requires a multiple
/// of 3.
pub fn multi_gaussian_model(x: &[f64], params: &[f64]) -> Vec<f64> {
    let inv = 1.0 / fwhm_to_sigma_factor();
    let mut y = vec![0.0_f64; x.len()];
    for triple in params.chunks_exact(3) {
        let (height, centroid, fwhm) = (triple[0], triple[1], triple[2]);
        let sigma = fwhm * inv;
        if sigma == 0.0 {
            continue;
        }
        for (yi, &xi) in y.iter_mut().zip(x.iter()) {
            let dhelp = (xi - centroid) / sigma;
            if dhelp <= 20.0 {
                *yi += height * (-0.5 * dhelp * dhelp).exp();
            }
        }
    }
    y
}

/// Estimate initial parameters and fit constraints for a multi-peak Gaussian fit
/// (silx `estimate_height_position_fwhm`, default config).
///
/// Locates peaks with [`peak_search`](crate::core::peaks::peak_search)
/// (`search_fwhm` floored at 3, `sensitivity` floored at 1; falls back to the
/// global maximum when none are found, silx `ForcePeakPresence`). Seeds each
/// peak as `(y[peak], x[peak], 5·|xspan|/n)`, refines the seeds with a quick
/// 4-iteration constrained fit (height/FWHM positive, centre quoted to ±½ of the
/// `search_fwhm`-sample x-width), and returns the refined seeds together with the
/// final per-parameter constraints (height & FWHM `Positive`, centre `Free` —
/// silx default `PositiveHeightAreaFlag`/`PositiveFwhmFlag` on, `QuotedPositionFlag`
/// off). Background removal (silx `StripBackgroundFlag`, off by default) is left
/// to the caller via [`crate::core::background`]. Returns `None` for empty or
/// length-mismatched input, or when no peak can be seeded.
pub fn estimate_multi_gaussian(
    x: &[f64],
    y: &[f64],
    search_fwhm: f64,
    sensitivity: f64,
) -> Option<(Vec<f64>, Vec<Constraint>)> {
    let npoints = y.len();
    if npoints == 0 || x.len() != npoints {
        return None;
    }
    let search_fwhm = search_fwhm.max(3.0);
    let search_sens = sensitivity.max(1.0);

    let found = crate::core::peaks::peak_search(y, search_fwhm, search_sens);
    let peaks: Vec<usize> = if found.is_empty() {
        // silx ForcePeakPresence: use the (first) global maximum.
        let maxv = y.iter().copied().fold(f64::NEG_INFINITY, f64::max);
        match y.iter().position(|&v| v == maxv) {
            Some(p) => vec![p],
            None => return None,
        }
    } else {
        found.iter().map(|p| p.index).collect()
    };
    if peaks.is_empty() {
        return None;
    }

    // silx seeds FWHM as 5 sampling intervals (in x units).
    let sig = 5.0 * (x[npoints - 1] - x[0]).abs() / npoints as f64;
    let mut param: Vec<f64> = Vec::with_capacity(peaks.len() * 3);
    let mut index_largest = 0usize;
    let mut height_largest = f64::NEG_INFINITY;
    for (k, &pi) in peaks.iter().enumerate() {
        let height = y[pi];
        // silx zeroes a near-zero position only for the first peak.
        let pos = if k == 0 && x[pi].abs() < 1.0e-16 {
            0.0
        } else {
            x[pi]
        };
        param.push(height);
        param.push(pos);
        param.push(sig);
        if height > height_largest {
            height_largest = height;
            index_largest = k;
        }
    }
    let _ = index_largest; // used by silx SameFwhmFlag (off by default).

    // Preliminary constraints for the quick refine: height & FWHM positive,
    // centre quoted around its seed.
    let sf = search_fwhm as usize;
    let (fwhmx, use_fwhmx) = if x.len() > sf {
        ((x[sf] - x[0]).abs(), true)
    } else {
        (0.0, false)
    };
    let xmin = x.iter().copied().fold(f64::INFINITY, f64::min);
    let xmax = x.iter().copied().fold(f64::NEG_INFINITY, f64::max);
    let mut prelim: Vec<Constraint> = Vec::with_capacity(param.len());
    for k in 0..peaks.len() {
        let pos = param[3 * k + 1];
        prelim.push(Constraint::Positive);
        if use_fwhmx && fwhmx > 0.0 {
            prelim.push(Constraint::Quoted {
                min: pos - 0.5 * fwhmx,
                max: pos + 0.5 * fwhmx,
            });
        } else if xmax > xmin {
            prelim.push(Constraint::Quoted {
                min: xmin,
                max: xmax,
            });
        } else {
            prelim.push(Constraint::Free);
        }
        prelim.push(Constraint::Positive);
    }

    // Quick 4-iteration refine (silx max_iter=4); fall back to the raw seeds.
    let fittedpar = leastsq_constrained(
        multi_gaussian_model,
        x,
        y,
        &param,
        &prelim,
        None,
        4,
        DEFAULT_DELTACHI,
    )
    .map(|r| r.parameters)
    .unwrap_or(param);

    // Final constraints (default config): height & FWHM positive, centre free.
    let mut cons: Vec<Constraint> = Vec::with_capacity(fittedpar.len());
    for _ in 0..peaks.len() {
        cons.push(Constraint::Positive);
        cons.push(Constraint::Free);
        cons.push(Constraint::Positive);
    }

    Some((fittedpar, cons))
}

/// Fit a sum of Gaussians to `(x, y)`, discovering the peaks automatically
/// (silx `FitManager` multi-peak Gaussian fit).
///
/// Seeds the peaks with [`estimate_multi_gaussian`] then runs the full
/// constrained Levenberg-Marquardt fit ([`leastsq_constrained`]) over all peaks
/// simultaneously. The returned [`LeastSqResult::parameters`] is the flat
/// `(height, centroid, fwhm)` triple vector accepted by [`multi_gaussian_model`].
/// Returns `None` when no peak can be seeded or the solver fails. Background is
/// the caller's responsibility (see [`crate::core::background`]).
pub fn fit_multi_gaussian(
    x: &[f64],
    y: &[f64],
    search_fwhm: f64,
    sensitivity: f64,
    max_iter: usize,
    deltachi: f64,
) -> Option<LeastSqResult> {
    let (seeds, cons) = estimate_multi_gaussian(x, y, search_fwhm, sensitivity)?;
    leastsq_constrained(
        multi_gaussian_model,
        x,
        y,
        &seeds,
        &cons,
        None,
        max_iter,
        deltachi,
    )
    .ok()
}

/// Run the multi-peak Gaussian fit ([`fit_multi_gaussian`]) and package it as an
/// [`IterativeFitResult`] for display.
///
/// The fitted curve is [`multi_gaussian_model`] evaluated over the located
/// peaks; the parameter vector is the flat `(height, centre, fwhm)` triples;
/// the per-peak names are `Height i` / `Center i` / `FWHM i` (1-based); and the
/// solver covariance / chi-square come straight from the constrained solve.
/// Returns `None` when no peak is seeded or the solver fails.
pub fn fit_multi_gaussian_full(
    x: &[f64],
    y: &[f64],
    search_fwhm: f64,
    sensitivity: f64,
    max_iter: usize,
    deltachi: f64,
) -> Option<IterativeFitResult> {
    let solver = fit_multi_gaussian(x, y, search_fwhm, sensitivity, max_iter, deltachi)?;
    let y_fit = multi_gaussian_model(x, &solver.parameters);
    let mut param_names = Vec::with_capacity(solver.parameters.len());
    for peak in 0..solver.parameters.len() / 3 {
        let i = peak + 1;
        param_names.push(format!("Height {i}"));
        param_names.push(format!("Center {i}"));
        param_names.push(format!("FWHM {i}"));
    }
    let fit = FitResult {
        y_fit,
        parameters: solver.parameters.clone(),
        param_names,
    };
    Some(IterativeFitResult { fit, solver })
}

/// Fit a single [`PeakModel`] from explicit initial parameters `p0` under
/// per-parameter [`Constraint`]s (silx FitWidget parameter table → editable
/// values + constraint codes feeding `leastsq_constrained`).
///
/// `p0` and `constraints` must each have exactly one entry per model parameter
/// (the order of [`PeakModel::param_names`]). Returns `None` on a length
/// mismatch or solver failure. This is the single owner of the
/// constrained-fit-from-`p0` path; [`fit_peak_constrained`] estimates `p0` and
/// delegates here.
pub fn fit_peak_from(
    model: PeakModel,
    x: &[f64],
    y: &[f64],
    p0: &[f64],
    constraints: &[Constraint],
    max_iter: usize,
    deltachi: f64,
) -> Option<IterativeFitResult> {
    if constraints.len() != p0.len() {
        return None;
    }
    let solver = leastsq_constrained(
        |xx, pp| model.eval(xx, pp),
        x,
        y,
        p0,
        constraints,
        None,
        max_iter,
        deltachi,
    )
    .ok()?;
    let y_fit = model.eval(x, &solver.parameters);
    let fit = FitResult {
        y_fit,
        parameters: solver.parameters.clone(),
        param_names: model.param_names(),
    };
    Some(IterativeFitResult { fit, solver })
}

/// Fit a single [`PeakModel`] under per-parameter [`Constraint`]s, with the
/// initial parameters estimated from the data ([`PeakModel::estimate`]).
///
/// `constraints` must have exactly one entry per model parameter. Returns `None`
/// on estimate/solver failure or a constraint-count mismatch. Delegates to
/// [`fit_peak_from`] once the estimate is in hand.
pub fn fit_peak_constrained(
    model: PeakModel,
    x: &[f64],
    y: &[f64],
    constraints: &[Constraint],
    max_iter: usize,
    deltachi: f64,
) -> Option<IterativeFitResult> {
    let p0 = model.estimate(x, y)?;
    fit_peak_from(model, x, y, &p0, constraints, max_iter, deltachi)
}

/// Outcome of [`fit_peak_with_background`]: the peak fit on the
/// background-subtracted residual, the estimated background curve, and the
/// total displayed curve.
#[derive(Debug, Clone)]
pub struct BackgroundPeakFit {
    /// The peak fit on the background-subtracted residual. `peak.fit.parameters`
    /// are in the [`PeakModel`]'s own parameterisation; per-parameter errors come
    /// from [`IterativeFitResult::std_errors`].
    pub peak: IterativeFitResult,
    /// The background curve estimated from the data and sampled at the fit `x`
    /// (held fixed during the peak fit), length `x.len()`.
    pub background: Vec<f64>,
    /// The total curve to display, `background[i] + peak.fit.y_fit[i]`.
    pub total: Vec<f64>,
}

/// Fit `model` on top of an estimated `background`, mirroring silx `FitManager`'s
/// background-then-peak workflow.
///
/// The background is estimated from the data with [`Background`] (silx
/// `estimate_*`: the strip/snip filters, or a polynomial least-squares-fitted to
/// the strip background — silx `EstimatePolyOnStrip = True`), subtracted from
/// `y`, and the [`PeakModel`] is fitted to the residual. The background is then
/// added back to produce [`BackgroundPeakFit::total`].
///
/// Returns `None` on length mismatch, empty input, or when the peak fit fails.
///
/// # Deviation from silx
///
/// silx's background theories are `is_background=True` `FitTheory`s whose
/// parameters are concatenated with the peak parameters into ONE `leastsq`, so
/// the analytic-background coefficients (Constant / Linear / Polynomial) are
/// refined *simultaneously* with the peak. Here the background is estimated once
/// and held fixed while the peak is fitted on the residual. A simultaneous
/// refinement is blocked by siplot's single-peak [`PeakModel`]s baking in their
/// own trailing constant-background parameter: a free analytic-background
/// constant would be collinear with it (singular covariance). Removing that
/// baked-in constant is a [`PeakModel`] redesign tracked separately.
///
/// [`Background`]: crate::core::background::Background
pub fn fit_peak_with_background(
    model: PeakModel,
    background: crate::core::background::Background,
    x: &[f64],
    y: &[f64],
    max_iter: usize,
    deltachi: f64,
) -> Option<BackgroundPeakFit> {
    if x.is_empty() || x.len() != y.len() {
        return None;
    }
    let bg = background.compute(x, y);
    let residual: Vec<f64> = y.iter().zip(&bg).map(|(&yi, &bi)| yi - bi).collect();
    let fitter = IterativeFit {
        model,
        max_iter,
        deltachi,
    };
    let peak = fitter.fit_full(x, &residual)?;
    let total: Vec<f64> = bg
        .iter()
        .zip(&peak.fit.y_fit)
        .map(|(&bi, &fi)| bi + fi)
        .collect();
    Some(BackgroundPeakFit {
        peak,
        background: bg,
        total,
    })
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::core::peaks::DEFAULT_PEAK_SENSITIVITY;

    /// Synthetic noiseless Gaussian sampled on a grid.
    fn synth_gaussian(xs: &[f64], height: f64, center: f64, fwhm: f64, bg: f64) -> Vec<f64> {
        gaussian_model(xs, &[height, center, fwhm, bg])
    }

    fn linspace(a: f64, b: f64, n: usize) -> Vec<f64> {
        (0..n)
            .map(|i| a + (b - a) * (i as f64) / ((n - 1) as f64))
            .collect()
    }

    #[test]
    fn invert_identity() {
        let id = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
        let inv = invert_matrix(&id).unwrap();
        assert_eq!(inv, id);
    }

    #[test]
    fn invert_known_2x2() {
        // [[4,7],[2,6]] inverse = [[0.6,-0.7],[-0.2,0.4]]
        let m = vec![vec![4.0, 7.0], vec![2.0, 6.0]];
        let inv = invert_matrix(&m).unwrap();
        let expected = [[0.6, -0.7], [-0.2, 0.4]];
        for i in 0..2 {
            for j in 0..2 {
                assert!((inv[i][j] - expected[i][j]).abs() < 1e-12);
            }
        }
    }

    #[test]
    fn invert_singular_returns_none() {
        let m = vec![vec![1.0, 2.0], vec![2.0, 4.0]];
        assert!(invert_matrix(&m).is_none());
    }

    #[test]
    fn leastsq_recovers_noiseless_line_exactly() {
        // Model: y = a*x + b, params [a, b]. Noiseless data with a=2.5, b=-1.0.
        let xs = linspace(-5.0, 5.0, 21);
        let (a_true, b_true) = (2.5, -1.0);
        let ys: Vec<f64> = xs.iter().map(|&x| a_true * x + b_true).collect();
        let model = |x: &[f64], p: &[f64]| x.iter().map(|&xi| p[0] * xi + p[1]).collect::<Vec<_>>();
        let res = leastsq(
            model,
            &xs,
            &ys,
            &[0.0, 0.0],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (res.parameters[0] - a_true).abs() < 1e-6,
            "slope {} vs {}",
            res.parameters[0],
            a_true
        );
        assert!(
            (res.parameters[1] - b_true).abs() < 1e-6,
            "intercept {} vs {}",
            res.parameters[1],
            b_true
        );
        // Noiseless → chisq essentially zero.
        assert!(res.chisq < 1e-12, "chisq {}", res.chisq);
    }

    #[test]
    fn leastsq_converges_on_noisy_gaussian() {
        // Synthetic gaussian + small deterministic "noise" so the test is
        // reproducible. height=10, center=2, fwhm=1.5, bg=1.
        let xs = linspace(-3.0, 7.0, 101);
        let clean = synth_gaussian(&xs, 10.0, 2.0, 1.5, 1.0);
        // Deterministic pseudo-noise: small sinusoidal perturbation.
        let ys: Vec<f64> = clean
            .iter()
            .enumerate()
            .map(|(i, &c)| c + 0.05 * ((i as f64) * 0.7).sin())
            .collect();
        let fit = IterativeFit::new(PeakModel::Gaussian)
            .fit_full(&xs, &ys)
            .expect("fit should succeed");
        let p = &fit.fit.parameters;
        assert!((p[0] - 10.0).abs() < 0.2, "height {}", p[0]);
        assert!((p[1] - 2.0).abs() < 0.05, "center {}", p[1]);
        assert!((p[2] - 1.5).abs() < 0.1, "fwhm {}", p[2]);
        assert!((p[3] - 1.0).abs() < 0.1, "bg {}", p[3]);
        // Reduced chi-square (sigma=1) is on the order of the perturbation
        // variance, not enormous.
        let rc = fit.reduced_chisq().unwrap();
        assert!(rc < 0.01, "reduced chisq {}", rc);
    }

    #[test]
    fn gaussian_model_recovers_own_peak() {
        let xs = linspace(0.0, 20.0, 201);
        let ys = synth_gaussian(&xs, 5.0, 8.0, 2.0, 0.5);
        let fit = IterativeFit::new(PeakModel::Gaussian)
            .fit_full(&xs, &ys)
            .unwrap();
        let p = &fit.fit.parameters;
        assert!((p[0] - 5.0).abs() < 1e-3, "height {}", p[0]);
        assert!((p[1] - 8.0).abs() < 1e-3, "center {}", p[1]);
        assert!((p[2] - 2.0).abs() < 1e-3, "fwhm {}", p[2]);
        assert!((p[3] - 0.5).abs() < 1e-3, "bg {}", p[3]);
        // Noiseless fit → reduced chisq near 0.
        assert!(fit.reduced_chisq().unwrap() < 1e-6);
    }

    #[test]
    fn gaussian_area_model_recovers_own_peak() {
        // Build data from the area model with a known area.
        let xs = linspace(0.0, 20.0, 201);
        let area = 12.0;
        let ys = gaussian_area_model(&xs, &[area, 9.0, 2.5, 0.2]);
        let fit = IterativeFit::new(PeakModel::GaussianArea)
            .fit_full(&xs, &ys)
            .unwrap();
        let p = &fit.fit.parameters;
        assert!((p[0] - area).abs() < 1e-2, "area {}", p[0]);
        assert!((p[1] - 9.0).abs() < 1e-3, "center {}", p[1]);
        assert!((p[2] - 2.5).abs() < 1e-3, "fwhm {}", p[2]);
        assert!((p[3] - 0.2).abs() < 1e-3, "bg {}", p[3]);
        assert!(fit.reduced_chisq().unwrap() < 1e-6);
    }

    #[test]
    fn lorentzian_model_recovers_own_peak() {
        let xs = linspace(0.0, 20.0, 201);
        let ys = lorentzian_model(&xs, &[7.0, 11.0, 3.0, 1.0]);
        let fit = IterativeFit::new(PeakModel::Lorentzian)
            .fit_full(&xs, &ys)
            .unwrap();
        let p = &fit.fit.parameters;
        assert!((p[0] - 7.0).abs() < 1e-2, "height {}", p[0]);
        assert!((p[1] - 11.0).abs() < 1e-3, "center {}", p[1]);
        assert!((p[2] - 3.0).abs() < 1e-2, "fwhm {}", p[2]);
        assert!((p[3] - 1.0).abs() < 1e-2, "bg {}", p[3]);
        assert!(fit.reduced_chisq().unwrap() < 1e-6);
    }

    #[test]
    fn pseudo_voigt_model_recovers_own_peak() {
        let xs = linspace(0.0, 20.0, 301);
        let ys = pseudo_voigt_model(&xs, &[6.0, 10.0, 2.0, 0.4, 0.5]);
        let fit = IterativeFit::new(PeakModel::PseudoVoigt)
            .fit_full(&xs, &ys)
            .unwrap();
        let p = &fit.fit.parameters;
        assert!((p[0] - 6.0).abs() < 5e-2, "height {}", p[0]);
        assert!((p[1] - 10.0).abs() < 1e-2, "center {}", p[1]);
        assert!((p[2] - 2.0).abs() < 5e-2, "fwhm {}", p[2]);
        assert!((p[3] - 0.4).abs() < 5e-2, "eta {}", p[3]);
        assert!((p[4] - 0.5).abs() < 5e-2, "bg {}", p[4]);
        assert!(fit.reduced_chisq().unwrap() < 1e-4);
    }

    #[test]
    fn pseudo_voigt_eta_limits_match_gauss_and_lorentz() {
        // eta=0 → pure Gaussian; eta=1 → pure Lorentzian (same height/center/fwhm).
        let xs = linspace(0.0, 10.0, 51);
        let g = gaussian_model(&xs, &[3.0, 5.0, 2.0, 0.0]);
        let pv_g = pseudo_voigt_model(&xs, &[3.0, 5.0, 2.0, 0.0, 0.0]);
        for (a, b) in g.iter().zip(pv_g.iter()) {
            assert!((a - b).abs() < 1e-12);
        }
        let l = lorentzian_model(&xs, &[3.0, 5.0, 2.0, 0.0]);
        let pv_l = pseudo_voigt_model(&xs, &[3.0, 5.0, 2.0, 1.0, 0.0]);
        for (a, b) in l.iter().zip(pv_l.iter()) {
            assert!((a - b).abs() < 1e-12);
        }
    }

    #[test]
    fn fit_in_range_ignores_outside_points() {
        // A clean gaussian inside [4, 12]; outside the range we plant a wildly
        // different curve. If out-of-range points were used, the fit would be
        // pulled away from the true peak.
        let xs = linspace(0.0, 20.0, 201);
        let in_range: Vec<f64> = xs
            .iter()
            .map(|&x| {
                if (4.0..=12.0).contains(&x) {
                    // true gaussian
                    let sigma = 2.0 / fwhm_to_sigma_factor();
                    let d = (x - 8.0) / sigma;
                    5.0 * (-0.5 * d * d).exp() + 0.5
                } else {
                    // garbage outside the range
                    100.0 + 50.0 * x
                }
            })
            .collect();
        let fitter = IterativeFit::new(PeakModel::Gaussian);
        let res = fit_in_range(&xs, &in_range, 4.0, 12.0, &fitter).unwrap();
        let p = &res.fit.parameters;
        assert!((p[1] - 8.0).abs() < 0.05, "center pulled to {}", p[1]);
        assert!((p[2] - 2.0).abs() < 0.1, "fwhm {}", p[2]);
        assert!((p[0] - 5.0).abs() < 0.2, "height {}", p[0]);
    }

    #[test]
    fn fit_in_range_reversed_bounds_equivalent() {
        let xs = linspace(0.0, 20.0, 201);
        let ys = synth_gaussian(&xs, 4.0, 10.0, 2.0, 0.3);
        let fitter = IterativeFit::new(PeakModel::Gaussian);
        let a = fit_in_range(&xs, &ys, 6.0, 14.0, &fitter).unwrap();
        let b = fit_in_range(&xs, &ys, 14.0, 6.0, &fitter).unwrap();
        for (pa, pb) in a.fit.parameters.iter().zip(b.fit.parameters.iter()) {
            assert!((pa - pb).abs() < 1e-12);
        }
    }

    #[test]
    fn std_errors_from_covariance_diagonal() {
        // Construct a LeastSqResult with a known covariance and verify the
        // error extraction (sqrt of the diagonal).
        let res = LeastSqResult {
            parameters: vec![1.0, 2.0, 3.0],
            covariance: vec![
                vec![4.0, 0.1, 0.0],
                vec![0.1, 9.0, 0.2],
                vec![0.0, 0.2, 16.0],
            ],
            uncertainties: vec![2.0, 3.0, 4.0],
            chisq: 0.0,
            reduced_chisq: Some(0.0),
            niter: 1,
            nfev: 1,
        };
        let errs = res.std_errors();
        assert!((errs[0] - 2.0).abs() < 1e-12);
        assert!((errs[1] - 3.0).abs() < 1e-12);
        assert!((errs[2] - 4.0).abs() < 1e-12);
    }

    #[test]
    fn std_errors_guard_negative_diagonal() {
        // A tiny negative diagonal (round-off) must not produce NaN; abs first.
        let res = LeastSqResult {
            parameters: vec![1.0],
            covariance: vec![vec![-1e-15]],
            uncertainties: vec![0.0],
            chisq: 0.0,
            reduced_chisq: None,
            niter: 0,
            nfev: 0,
        };
        let e = res.std_errors();
        assert!(e[0].is_finite() && e[0] >= 0.0);
    }

    #[test]
    fn leastsq_length_mismatch_errors() {
        let r = leastsq(
            |x: &[f64], _p: &[f64]| x.to_vec(),
            &[1.0, 2.0, 3.0],
            &[1.0, 2.0],
            &[0.0],
            None,
            10,
            DEFAULT_DELTACHI,
        );
        assert_eq!(r.unwrap_err(), FitError::LengthMismatch);
    }

    #[test]
    fn leastsq_rejects_nonfinite() {
        let r = leastsq(
            |x: &[f64], p: &[f64]| x.iter().map(|&xi| p[0] * xi).collect::<Vec<_>>(),
            &[1.0, f64::NAN, 3.0],
            &[1.0, 2.0, 3.0],
            &[1.0],
            None,
            10,
            DEFAULT_DELTACHI,
        );
        assert_eq!(r.unwrap_err(), FitError::NonFinite);
    }

    #[test]
    fn estimate_seeds_are_close() {
        let xs = linspace(0.0, 20.0, 201);
        let ys = synth_gaussian(&xs, 5.0, 8.0, 2.0, 0.5);
        let (h, c, f, bg) = estimate_height_position_fwhm(&xs, &ys).unwrap();
        assert!((h - 5.0).abs() < 0.5, "height seed {}", h);
        assert!((c - 8.0).abs() < 0.2, "center seed {}", c);
        assert!((f - 2.0).abs() < 0.5, "fwhm seed {}", f);
        assert!((bg - 0.5).abs() < 0.1, "bg seed {}", bg);
    }

    // --- Constrained leastsq (silx constraint codes) -----------------------

    // Straight line y = p[0]*x + p[1].
    fn line(x: &[f64], p: &[f64]) -> Vec<f64> {
        x.iter().map(|&xi| p[0] * xi + p[1]).collect()
    }
    // Constant y = p[0].
    fn constant(x: &[f64], p: &[f64]) -> Vec<f64> {
        vec![p[0]; x.len()]
    }

    #[test]
    fn constrained_all_free_matches_unconstrained() {
        let xs = linspace(-5.0, 5.0, 21);
        let ys: Vec<f64> = xs.iter().map(|&x| 2.5 * x - 1.0).collect();
        let free = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[0.0, 0.0],
            &[Constraint::Free, Constraint::Free],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let plain = leastsq(
            line,
            &xs,
            &ys,
            &[0.0, 0.0],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!((free.parameters[0] - plain.parameters[0]).abs() < 1e-6);
        assert!((free.parameters[1] - plain.parameters[1]).abs() < 1e-6);
        assert!((free.parameters[0] - 2.5).abs() < 1e-6);
        assert!((free.parameters[1] + 1.0).abs() < 1e-6);
    }

    #[test]
    fn constrained_fixed_holds_parameter() {
        let xs = linspace(-5.0, 5.0, 21);
        let ys: Vec<f64> = xs.iter().map(|&x| 2.5 * x - 1.0).collect();
        // b fixed at its true value -1.0; only a is fitted and must recover 2.5.
        let res = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[0.0, -1.0],
            &[Constraint::Free, Constraint::Fixed],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert_eq!(res.parameters[1], -1.0, "fixed b must not move");
        assert!(
            (res.parameters[0] - 2.5).abs() < 1e-6,
            "a {}",
            res.parameters[0]
        );
    }

    #[test]
    fn constrained_fixed_gets_full_uncertainty() {
        let xs = linspace(-5.0, 5.0, 21);
        let ys: Vec<f64> = xs.iter().map(|&x| 2.5 * x - 1.0).collect();
        let res = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[0.0, -1.0],
            &[Constraint::Free, Constraint::Fixed],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        // silx: fixed parameter gets covariance diag = value^2 and uncertainty = value.
        assert_eq!(res.uncertainties[1], res.parameters[1]);
        assert!((res.covariance[1][1] - res.parameters[1] * res.parameters[1]).abs() < 1e-12);
    }

    #[test]
    fn constrained_positive_enforces_and_recovers() {
        let xs = linspace(0.0, 10.0, 21);
        // Negative target: the abs reparametrisation keeps the parameter
        // non-negative, so it can never reach -3 and the residual stays large.
        let neg: Vec<f64> = vec![-3.0; xs.len()];
        let r_neg = leastsq_constrained(
            constant,
            &xs,
            &neg,
            &[1.0],
            &[Constraint::Positive],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            r_neg.parameters[0] >= 0.0,
            "positive violated: {}",
            r_neg.parameters[0]
        );
        assert!(
            r_neg.chisq > 1.0,
            "constraint should prevent fitting the negative target, chisq {}",
            r_neg.chisq
        );
        // Positive target: behaves like a normal fit, recovering 4.
        let pos: Vec<f64> = vec![4.0; xs.len()];
        let r_pos = leastsq_constrained(
            constant,
            &xs,
            &pos,
            &[1.0],
            &[Constraint::Positive],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (r_pos.parameters[0] - 4.0).abs() < 0.1,
            "recover {}",
            r_pos.parameters[0]
        );
    }

    #[test]
    fn constrained_quoted_clamps_to_bounds() {
        let xs = linspace(0.0, 10.0, 21);
        // Target 10 is above the [0,5] bound: the fit saturates near 5.
        let high: Vec<f64> = vec![10.0; xs.len()];
        let r_hi = leastsq_constrained(
            constant,
            &xs,
            &high,
            &[2.5],
            &[Constraint::Quoted { min: 0.0, max: 5.0 }],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (0.0..=5.0).contains(&r_hi.parameters[0]),
            "out of bounds: {}",
            r_hi.parameters[0]
        );
        assert!(
            r_hi.parameters[0] > 4.5,
            "did not saturate near 5: {}",
            r_hi.parameters[0]
        );
        // Target 3 is inside the bound and is recovered.
        let mid: Vec<f64> = vec![3.0; xs.len()];
        let r_mid = leastsq_constrained(
            constant,
            &xs,
            &mid,
            &[2.5],
            &[Constraint::Quoted { min: 0.0, max: 5.0 }],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (r_mid.parameters[0] - 3.0).abs() < 0.05,
            "recover {}",
            r_mid.parameters[0]
        );
    }

    #[test]
    fn constrained_factor_ties_parameters() {
        let xs = linspace(-5.0, 5.0, 21);
        // y = a*x + b with b = 2*a; true a = 3 => b = 6.
        let ys: Vec<f64> = xs.iter().map(|&x| 3.0 * x + 6.0).collect();
        let res = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[1.0, 0.0],
            &[
                Constraint::Free,
                Constraint::Factor {
                    reference: 0,
                    factor: 2.0,
                },
            ],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (res.parameters[0] - 3.0).abs() < 1e-4,
            "a {}",
            res.parameters[0]
        );
        assert!(
            (res.parameters[1] - 6.0).abs() < 1e-4,
            "b {}",
            res.parameters[1]
        );
        assert!(
            (res.parameters[1] - 2.0 * res.parameters[0]).abs() < 1e-9,
            "tie broken"
        );
    }

    #[test]
    fn constrained_delta_ties_parameters() {
        let xs = linspace(-5.0, 5.0, 21);
        // b = a + 5; true a = 2 => b = 7.
        let ys: Vec<f64> = xs.iter().map(|&x| 2.0 * x + 7.0).collect();
        let res = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[0.0, 0.0],
            &[
                Constraint::Free,
                Constraint::Delta {
                    reference: 0,
                    delta: 5.0,
                },
            ],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (res.parameters[0] - 2.0).abs() < 1e-4,
            "a {}",
            res.parameters[0]
        );
        assert!(
            (res.parameters[1] - res.parameters[0] - 5.0).abs() < 1e-9,
            "tie broken"
        );
    }

    #[test]
    fn constrained_sum_ties_parameters() {
        let xs = linspace(-5.0, 5.0, 21);
        // b = 10 - a; true a = 4 => b = 6.
        let ys: Vec<f64> = xs.iter().map(|&x| 4.0 * x + 6.0).collect();
        let res = leastsq_constrained(
            line,
            &xs,
            &ys,
            &[0.0, 0.0],
            &[
                Constraint::Free,
                Constraint::Sum {
                    reference: 0,
                    sum: 10.0,
                },
            ],
            None,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert!(
            (res.parameters[0] - 4.0).abs() < 1e-4,
            "a {}",
            res.parameters[0]
        );
        assert!(
            (res.parameters[0] + res.parameters[1] - 10.0).abs() < 1e-9,
            "tie broken"
        );
    }

    #[test]
    fn constrained_rejects_bad_spec() {
        let xs = linspace(0.0, 4.0, 5);
        let ys = vec![1.0; 5];
        // constraints length mismatch.
        assert_eq!(
            leastsq_constrained(constant, &xs, &ys, &[1.0], &[], None, 10, DEFAULT_DELTACHI)
                .unwrap_err(),
            FitError::BadConstraintReference
        );
        // equal Quoted bounds (B == 0).
        assert_eq!(
            leastsq_constrained(
                constant,
                &xs,
                &ys,
                &[1.0],
                &[Constraint::Quoted { min: 5.0, max: 5.0 }],
                None,
                10,
                DEFAULT_DELTACHI,
            )
            .unwrap_err(),
            FitError::InvalidConstraint
        );
        // Factor referencing a non-existent parameter.
        assert_eq!(
            leastsq_constrained(
                line,
                &xs,
                &ys,
                &[1.0, 0.0],
                &[
                    Constraint::Free,
                    Constraint::Factor {
                        reference: 9,
                        factor: 2.0
                    }
                ],
                None,
                10,
                DEFAULT_DELTACHI,
            )
            .unwrap_err(),
            FitError::BadConstraintReference
        );
        // No free parameters at all.
        assert_eq!(
            leastsq_constrained(
                line,
                &xs,
                &ys,
                &[1.0, 2.0],
                &[Constraint::Fixed, Constraint::Fixed],
                None,
                10,
                DEFAULT_DELTACHI,
            )
            .unwrap_err(),
            FitError::NoFreeParameters
        );
    }

    // --- Multi-peak (simultaneous) Gaussian fit ----------------------------

    fn grid(n: usize) -> Vec<f64> {
        (0..n).map(|i| i as f64).collect()
    }

    /// Find the fitted triple whose centre is nearest `target`.
    fn nearest_peak(params: &[f64], target: f64) -> [f64; 3] {
        params
            .chunks_exact(3)
            .min_by(|a, b| {
                (a[1] - target)
                    .abs()
                    .partial_cmp(&(b[1] - target).abs())
                    .unwrap()
            })
            .map(|t| [t[0], t[1], t[2]])
            .unwrap()
    }

    #[test]
    fn multi_gaussian_model_is_sum_of_single_gaussians() {
        let xs = grid(100);
        // Two peaks; compare against two single gaussian_model calls (bg = 0).
        let a = gaussian_model(&xs, &[100.0, 30.0, 8.0, 0.0]);
        let b = gaussian_model(&xs, &[60.0, 70.0, 5.0, 0.0]);
        let sum = multi_gaussian_model(&xs, &[100.0, 30.0, 8.0, 60.0, 70.0, 5.0]);
        for i in 0..xs.len() {
            assert!((sum[i] - (a[i] + b[i])).abs() < 1e-9, "mismatch at {i}");
        }
    }

    #[test]
    fn fit_multi_gaussian_recovers_two_peaks() {
        let xs = grid(100);
        let mut ys = gaussian_model(&xs, &[100.0, 30.0, 8.0, 0.0]);
        for (yi, g) in ys
            .iter_mut()
            .zip(gaussian_model(&xs, &[80.0, 70.0, 6.0, 0.0]))
        {
            *yi += g;
        }
        let res = fit_multi_gaussian(
            &xs,
            &ys,
            8.0,
            DEFAULT_PEAK_SENSITIVITY,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .expect("multi-peak fit should succeed");
        assert!(res.parameters.len() >= 6, "expected >=2 peaks");
        let p1 = nearest_peak(&res.parameters, 30.0);
        let p2 = nearest_peak(&res.parameters, 70.0);
        assert!((p1[1] - 30.0).abs() < 1.0, "centre1 {}", p1[1]);
        assert!((p1[0] - 100.0).abs() < 5.0, "height1 {}", p1[0]);
        assert!((p1[2] - 8.0).abs() < 1.0, "fwhm1 {}", p1[2]);
        assert!((p2[1] - 70.0).abs() < 1.0, "centre2 {}", p2[1]);
        assert!((p2[0] - 80.0).abs() < 5.0, "height2 {}", p2[0]);
        assert!((p2[2] - 6.0).abs() < 1.0, "fwhm2 {}", p2[2]);
    }

    #[test]
    fn fit_multi_gaussian_full_packages_names_errors_and_curve() {
        let xs = grid(100);
        let mut ys = gaussian_model(&xs, &[100.0, 30.0, 8.0, 0.0]);
        for (yi, g) in ys
            .iter_mut()
            .zip(gaussian_model(&xs, &[80.0, 70.0, 6.0, 0.0]))
        {
            *yi += g;
        }
        let ir = fit_multi_gaussian_full(
            &xs,
            &ys,
            8.0,
            DEFAULT_PEAK_SENSITIVITY,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .expect("multi-peak fit should succeed");
        let n = ir.fit.parameters.len();
        assert!(n >= 6 && n.is_multiple_of(3), "param count {n}");
        // Names, values, and errors all line up for the results table.
        assert_eq!(ir.fit.param_names.len(), n);
        assert_eq!(ir.std_errors().len(), n);
        assert_eq!(ir.fit.y_fit.len(), xs.len());
        // Per-peak naming is 1-based Height/Center/FWHM triples.
        assert_eq!(ir.fit.param_names[0], "Height 1");
        assert_eq!(ir.fit.param_names[1], "Center 1");
        assert_eq!(ir.fit.param_names[2], "FWHM 1");
        // The packaged curve equals the model evaluated at the fitted params.
        assert_eq!(ir.fit.y_fit, multi_gaussian_model(&xs, &ir.fit.parameters));
        // Both peaks recovered.
        assert!((nearest_peak(&ir.fit.parameters, 30.0)[1] - 30.0).abs() < 1.0);
        assert!((nearest_peak(&ir.fit.parameters, 70.0)[1] - 70.0).abs() < 1.0);
    }

    #[test]
    fn fit_peak_constrained_all_free_recovers_peak() {
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[60.0, 45.0, 7.0, 5.0]);
        let cons = vec![Constraint::Free; 4];
        let ir = fit_peak_constrained(
            PeakModel::Gaussian,
            &xs,
            &ys,
            &cons,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let p = &ir.fit.parameters;
        assert!((p[1] - 45.0).abs() < 1.0, "centre {}", p[1]);
        assert!((p[2] - 7.0).abs() < 1.0, "fwhm {}", p[2]);
        assert_eq!(ir.fit.param_names, PeakModel::Gaussian.param_names());
    }

    #[test]
    fn fit_peak_constrained_fixed_holds_param_at_estimate() {
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[60.0, 45.0, 7.0, 5.0]);
        // The fixed parameter must stay bit-identical to its estimate.
        let p0 = PeakModel::Gaussian.estimate(&xs, &ys).unwrap();
        let mut cons = vec![Constraint::Free; 4];
        cons[1] = Constraint::Fixed; // hold the centre
        let ir = fit_peak_constrained(
            PeakModel::Gaussian,
            &xs,
            &ys,
            &cons,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert_eq!(ir.fit.parameters[1], p0[1]);
    }

    #[test]
    fn fit_peak_from_uses_explicit_p0_and_recovers() {
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[60.0, 45.0, 7.0, 5.0]);
        // A deliberately-off initial guess; LM should still converge.
        let p0 = [20.0, 40.0, 12.0, 0.0];
        let cons = vec![Constraint::Free; 4];
        let ir = fit_peak_from(
            PeakModel::Gaussian,
            &xs,
            &ys,
            &p0,
            &cons,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let p = &ir.fit.parameters;
        assert!((p[1] - 45.0).abs() < 1.0, "centre {}", p[1]);
        assert!((p[2] - 7.0).abs() < 1.0, "fwhm {}", p[2]);
        // fit_peak_constrained estimates p0 then delegates here, so the two
        // agree when started from the estimate.
        let est = PeakModel::Gaussian.estimate(&xs, &ys).unwrap();
        let via_estimate = fit_peak_from(
            PeakModel::Gaussian,
            &xs,
            &ys,
            &est,
            &cons,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let direct = fit_peak_constrained(
            PeakModel::Gaussian,
            &xs,
            &ys,
            &cons,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        assert_eq!(via_estimate.fit.parameters, direct.fit.parameters);
    }

    #[test]
    fn fit_peak_constrained_rejects_count_mismatch() {
        let xs = grid(50);
        let ys = gaussian_model(&xs, &[60.0, 25.0, 7.0, 0.0]);
        // Gaussian has 4 parameters; a 3-entry constraint vector is rejected.
        assert!(
            fit_peak_constrained(
                PeakModel::Gaussian,
                &xs,
                &ys,
                &[Constraint::Free; 3],
                DEFAULT_MAX_ITER,
                DEFAULT_DELTACHI,
            )
            .is_none()
        );
    }

    #[test]
    fn fit_multi_gaussian_single_peak() {
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[50.0, 45.0, 7.0, 0.0]);
        let res = fit_multi_gaussian(
            &xs,
            &ys,
            7.0,
            DEFAULT_PEAK_SENSITIVITY,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let p = nearest_peak(&res.parameters, 45.0);
        assert!((p[0] - 50.0).abs() < 2.0, "height {}", p[0]);
        assert!((p[1] - 45.0).abs() < 0.5, "centre {}", p[1]);
        assert!((p[2] - 7.0).abs() < 0.5, "fwhm {}", p[2]);
    }

    #[test]
    fn estimate_multi_gaussian_seeds_height_and_position() {
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[50.0, 45.0, 7.0, 0.0]);
        let (seeds, cons) =
            estimate_multi_gaussian(&xs, &ys, 7.0, DEFAULT_PEAK_SENSITIVITY).unwrap();
        assert_eq!(seeds.len() % 3, 0);
        assert_eq!(seeds.len(), cons.len());
        // Final constraints: height Positive, centre Free, FWHM Positive.
        assert_eq!(cons[0], Constraint::Positive);
        assert_eq!(cons[1], Constraint::Free);
        assert_eq!(cons[2], Constraint::Positive);
    }

    #[test]
    fn estimate_multi_gaussian_rejects_empty_and_mismatch() {
        assert!(estimate_multi_gaussian(&[], &[], 5.0, 2.5).is_none());
        assert!(estimate_multi_gaussian(&[0.0, 1.0], &[1.0], 5.0, 2.5).is_none());
    }

    // --- Non-peak theories: erf, step/slit/atan models, estimators ---------

    #[test]
    fn erf_matches_known_values() {
        // Reference values to full precision; A&S 7.1.26 is good to ~1.5e-7.
        assert_eq!(erf(0.0), 0.0);
        assert!((erf(0.5) - 0.520_499_877_813_046_5).abs() < 1e-6);
        assert!((erf(1.0) - 0.842_700_792_949_714_9).abs() < 1e-6);
        assert!((erf(2.0) - 0.995_322_265_018_952_7).abs() < 1e-6);
        // Odd symmetry.
        assert!((erf(-1.0) + erf(1.0)).abs() < 1e-12);
    }

    #[test]
    fn erfc_is_one_minus_erf() {
        for &x in &[-2.0, -0.3, 0.0, 0.7, 1.5] {
            assert!((erfc(x) - (1.0 - erf(x))).abs() < 1e-15);
        }
    }

    #[test]
    fn convolve_valid_matches_numpy_reversed_kernel() {
        // numpy reverses the kernel: out[k] = sum_j y[k+j] * kernel[m-1-j].
        // [1,0,0] picks y[k+2], so valid output is [y[2], y[3]].
        assert_eq!(
            convolve_valid(&[1.0, 2.0, 3.0, 4.0], &[1.0, 0.0, 0.0]),
            vec![3.0, 4.0]
        );
        // Step-up edge kernel on a ramp: single 'valid' output = 2.5.
        let out = convolve_valid(&[1.0, 2.0, 3.0, 4.0, 5.0], &[0.25, 0.75, 0.0, -0.75, -0.25]);
        assert_eq!(out.len(), 1);
        assert!((out[0] - 2.5).abs() < 1e-12);
        // Kernel longer than data → empty.
        assert!(convolve_valid(&[1.0], &[1.0, 2.0]).is_empty());
    }

    #[test]
    fn stepup_model_half_height_at_center_and_asymptotes() {
        let p = [10.0, 0.0, 4.0, 2.0]; // height, center, fwhm, bg
        // erf(0)=0 exactly → exactly bg + height/2 at the centre.
        assert_eq!(stepup_model(&[0.0], &p)[0], 7.0);
        // Asymptotes: bg on the far left, bg+height on the far right.
        assert!((stepup_model(&[-1000.0], &p)[0] - 2.0).abs() < 1e-9);
        assert!((stepup_model(&[1000.0], &p)[0] - 12.0).abs() < 1e-9);
        // Monotone increasing.
        let xs = linspace(-20.0, 20.0, 81);
        let ys = stepup_model(&xs, &p);
        assert!(ys.windows(2).all(|w| w[1] >= w[0]));
    }

    #[test]
    fn stepdown_model_half_height_at_center_and_asymptotes() {
        let p = [10.0, 0.0, 4.0, 2.0];
        assert_eq!(stepdown_model(&[0.0], &p)[0], 7.0); // erfc(0)=1
        assert!((stepdown_model(&[-1000.0], &p)[0] - 12.0).abs() < 1e-9);
        assert!((stepdown_model(&[1000.0], &p)[0] - 2.0).abs() < 1e-9);
        let xs = linspace(-20.0, 20.0, 81);
        let ys = stepdown_model(&xs, &p);
        assert!(ys.windows(2).all(|w| w[1] <= w[0]));
    }

    #[test]
    fn atan_stepup_model_half_height_at_center_and_monotonic() {
        let p = [10.0, 0.0, 4.0, 2.0]; // height, position, width, bg
        assert_eq!(atan_stepup_model(&[0.0], &p)[0], 7.0); // atan(0)=0
        assert!((atan_stepup_model(&[-1.0e8], &p)[0] - 2.0).abs() < 1e-6);
        assert!((atan_stepup_model(&[1.0e8], &p)[0] - 12.0).abs() < 1e-6);
        let xs = linspace(-50.0, 50.0, 101);
        let ys = atan_stepup_model(&xs, &p);
        assert!(ys.windows(2).all(|w| w[1] >= w[0]));
    }

    #[test]
    fn slit_model_is_symmetric_and_localised() {
        let p = [10.0, 0.0, 10.0, 2.0, 1.0]; // height, position, fwhm, beamfwhm, bg
        // Symmetric about the position.
        for &d in &[1.0, 3.0, 7.0, 12.0] {
            let a = slit_model(&[d], &p)[0];
            let b = slit_model(&[-d], &p)[0];
            assert!((a - b).abs() < 1e-12, "slit asymmetric at {d}: {a} vs {b}");
        }
        // Inside the slit it sits above the background; far outside it is bg.
        assert!(slit_model(&[0.0], &p)[0] > 5.0);
        assert!((slit_model(&[1000.0], &p)[0] - 1.0).abs() < 1e-9);
        assert!((slit_model(&[-1000.0], &p)[0] - 1.0).abs() < 1e-9);
    }

    #[test]
    fn estimate_stepup_recovers_center_height_bg() {
        let xs = linspace(0.0, 100.0, 101);
        let ys = stepup_model(&xs, &[8.0, 40.0, 6.0, 3.0]);
        let s = estimate_stepup(&xs, &ys).unwrap();
        assert!((s[0] - 8.0).abs() < 0.5, "height seed {}", s[0]);
        assert!((s[1] - 40.0).abs() < 3.0, "center seed {}", s[1]);
        assert!((s[3] - 3.0).abs() < 0.5, "bg seed {}", s[3]);
    }

    #[test]
    fn estimate_stepdown_recovers_center_height_bg() {
        let xs = linspace(0.0, 100.0, 101);
        let ys = stepdown_model(&xs, &[8.0, 40.0, 6.0, 3.0]);
        let s = estimate_stepdown(&xs, &ys).unwrap();
        assert!((s[0] - 8.0).abs() < 0.5, "height seed {}", s[0]);
        assert!((s[1] - 40.0).abs() < 3.0, "center seed {}", s[1]);
        assert!((s[3] - 3.0).abs() < 0.5, "bg seed {}", s[3]);
    }

    #[test]
    fn estimate_slit_centers_on_the_slit() {
        let xs = linspace(0.0, 100.0, 201);
        let ys = slit_model(&xs, &[10.0, 50.0, 20.0, 3.0, 2.0]);
        let s = estimate_slit(&xs, &ys).unwrap();
        assert_eq!(s.len(), 5);
        assert!((s[1] - 50.0).abs() < 2.0, "position seed {}", s[1]);
        assert!(s[3] > 0.0, "beamfwhm seed must be positive: {}", s[3]);
    }

    #[test]
    fn peak_model_step_variants_delegate() {
        assert_eq!(PeakModel::StepUp.name(), "Step Up");
        assert_eq!(PeakModel::StepDown.name(), "Step Down");
        assert_eq!(PeakModel::Slit.name(), "Slit");
        assert_eq!(PeakModel::AtanStepUp.name(), "Arctan Step Up");
        assert_eq!(
            PeakModel::StepUp.param_names(),
            vec!["Height", "Center", "FWHM", "Background"]
        );
        assert_eq!(PeakModel::Slit.param_names().len(), 5);
        assert_eq!(
            PeakModel::AtanStepUp.param_names(),
            vec!["Height", "Center", "Width", "Background"]
        );
        // eval delegates to the matching model function.
        let xs = linspace(-5.0, 5.0, 11);
        let p = [3.0, 0.0, 2.0, 1.0];
        assert_eq!(PeakModel::StepUp.eval(&xs, &p), stepup_model(&xs, &p));
        assert_eq!(PeakModel::StepDown.eval(&xs, &p), stepdown_model(&xs, &p));
        assert_eq!(
            PeakModel::AtanStepUp.eval(&xs, &p),
            atan_stepup_model(&xs, &p)
        );
    }

    #[test]
    fn iterative_fit_recovers_stepup() {
        let xs = linspace(0.0, 100.0, 101);
        let truth = [8.0, 40.0, 6.0, 3.0];
        let ys = stepup_model(&xs, &truth);
        let fit = IterativeFit::new(PeakModel::StepUp)
            .fit_full(&xs, &ys)
            .unwrap();
        let p = &fit.fit.parameters;
        assert!((p[0] - truth[0]).abs() < 1.0, "height {}", p[0]);
        assert!((p[1] - truth[1]).abs() < 1.0, "center {}", p[1]);
        assert!((p[3] - truth[3]).abs() < 1.0, "bg {}", p[3]);
    }

    // --- Peak-on-background fit (silx background-then-peak workflow) --------

    #[test]
    fn fit_peak_with_background_none_is_byte_identical_to_plain_fit() {
        use crate::core::background::Background;
        let xs = grid(100);
        let ys = gaussian_model(&xs, &[100.0, 50.0, 8.0, 0.0]);
        let bgfit = fit_peak_with_background(
            PeakModel::Gaussian,
            Background::None,
            &xs,
            &ys,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        let plain = IterativeFit::new(PeakModel::Gaussian)
            .fit_full(&xs, &ys)
            .unwrap();
        // No background: residual == y exactly, so the peak fit and the total
        // curve match the plain iterative fit bit-for-bit.
        assert!(bgfit.background.iter().all(|&b| b == 0.0));
        assert_eq!(bgfit.peak.fit.parameters, plain.fit.parameters);
        assert_eq!(bgfit.total, bgfit.peak.fit.y_fit);
        assert_eq!(bgfit.total, plain.fit.y_fit);
    }

    #[test]
    fn fit_peak_with_background_constant_offset_recovers_peak() {
        use crate::core::background::Background;
        let xs = grid(100);
        // Gaussian sitting on a +25 constant pedestal (the gaussian model's
        // trailing parameter is its own constant background).
        let ys = gaussian_model(&xs, &[100.0, 50.0, 8.0, 25.0]);
        let bgfit = fit_peak_with_background(
            PeakModel::Gaussian,
            Background::Constant,
            &xs,
            &ys,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        // Constant background = min(y) ≈ the 25 pedestal (tails approach it).
        assert!(
            (bgfit.background[0] - 25.0).abs() < 1.0,
            "background {}",
            bgfit.background[0]
        );
        assert!(bgfit.background.iter().all(|&b| b == bgfit.background[0]));
        // Centre recovered, and the reconstructed total tracks the data.
        let p = &bgfit.peak.fit.parameters;
        assert!((p[1] - 50.0).abs() < 2.0, "centre {}", p[1]);
        let max_err = ys
            .iter()
            .zip(&bgfit.total)
            .map(|(&d, &t)| (d - t).abs())
            .fold(0.0_f64, f64::max);
        assert!(max_err < 5.0, "max |data - total| = {max_err}");
    }

    #[test]
    fn fit_peak_with_background_linear_recovers_peak_on_slope() {
        use crate::core::background::Background;
        let xs = grid(120);
        // Gaussian on a rising line y = 0.5*x + 10.
        let base = line(&xs, &[0.5, 10.0]);
        let peak = gaussian_model(&xs, &[80.0, 60.0, 8.0, 0.0]);
        let ys: Vec<f64> = base.iter().zip(&peak).map(|(&b, &p)| b + p).collect();
        let bgfit = fit_peak_with_background(
            PeakModel::Gaussian,
            Background::Linear,
            &xs,
            &ys,
            DEFAULT_MAX_ITER,
            DEFAULT_DELTACHI,
        )
        .unwrap();
        // The fitted background follows the positive slope.
        assert!(
            *bgfit.background.last().unwrap() > bgfit.background[0],
            "background not rising: {} -> {}",
            bgfit.background[0],
            bgfit.background.last().unwrap()
        );
        // Peak centre and width recovered on the residual.
        let p = &bgfit.peak.fit.parameters;
        assert!((p[1] - 60.0).abs() < 3.0, "centre {}", p[1]);
        assert!((p[2] - 8.0).abs() < 4.0, "fwhm {}", p[2]);
        // Total tracks the data (peak height 80 → allow a modest residual).
        let max_err = ys
            .iter()
            .zip(&bgfit.total)
            .map(|(&d, &t)| (d - t).abs())
            .fold(0.0_f64, f64::max);
        assert!(max_err < 12.0, "max |data - total| = {max_err}");
    }

    #[test]
    fn fit_peak_with_background_rejects_bad_input() {
        use crate::core::background::Background;
        let xs = grid(10);
        let ys = constant(&xs, &[1.0]);
        // Length mismatch.
        assert!(
            fit_peak_with_background(
                PeakModel::Gaussian,
                Background::None,
                &xs,
                &ys[..5],
                DEFAULT_MAX_ITER,
                DEFAULT_DELTACHI,
            )
            .is_none()
        );
        // Empty input.
        assert!(
            fit_peak_with_background(
                PeakModel::Gaussian,
                Background::None,
                &[],
                &[],
                DEFAULT_MAX_ITER,
                DEFAULT_DELTACHI,
            )
            .is_none()
        );
    }
}