nereids_fitting/

poisson.rs

//! Poisson-likelihood fitting of counts.
//!
//! Minimizes half the Poisson deviance plus a Gaussian prior on each block
//! `θₚ` of measured parameters,
//!
//! ```text
//! D(θ) = Σᵢ [μᵢ(θ) − yᵢ + yᵢ ln(yᵢ / μᵢ(θ))] + ½ Σₚ (θₚ − mₚ)ᵀ Cₚ⁻¹ (θₚ − mₚ)
//! ```
//!
//! within box bounds, by projected Levenberg–Marquardt steps in the Fisher
//! metric, for models with an analytical Jacobian.  A prior adds `Cₚ⁻¹` to the
//! information, as SAMMY's Bayes update `M′⁻¹ = M⁻¹ + GᵀV⁻¹G` adds the data to
//! a prior covariance `M` (SAMMY manual, `docs/tex/fitting-procedure.tex`,
//! eq. `bayes-eq-m-prime`).

use faer::dyn_stack::{MemBuffer, MemStack};
use faer::linalg::svd::{ComputeSvdVectors, svd, svd_scratch};
use faer::{Mat, Par, Side};

use crate::error::FittingError;
use crate::lm::{FitModel, FlatMatrix};
use crate::parameters::{FitParameter, ParameterSet};

/// Configuration for the Poisson solvers.
#[derive(Debug, Clone)]
pub struct PoissonConfig {
    /// Maximum number of steps.
    pub max_iter: usize,
    /// Whether to compute the covariance and error bars after convergence.
    pub compute_covariance: bool,
}

impl Default for PoissonConfig {
    fn default() -> Self {
        Self {
            max_iter: 200,
            compute_covariance: true,
        }
    }
}

/// A Gaussian prior on free parameters: a measurement of them made
/// elsewhere, with mean `m` and covariance `C = LLᵀ`.
#[derive(Debug, Clone, PartialEq)]
pub struct Prior {
    pub(crate) parameters: Vec<usize>,
    pub(crate) mean: Vec<f64>,
    pub(crate) factor: Vec<f64>,
}

impl Prior {
    /// A measurement `mean ± sd` of the parameter at index `parameter` in the
    /// [`ParameterSet`]; [`poisson_fit`] refuses a mean that is not finite or
    /// an sd that is not finite and positive.
    pub fn measured(parameter: usize, mean: f64, sd: f64) -> Self {
        Self {
            parameters: vec![parameter],
            mean: vec![mean],
            factor: vec![sd],
        }
    }

    /// A joint measurement of the parameters at `parameters` in the
    /// [`ParameterSet`], with `mean` and `covariance` in that order.
    ///
    /// # Errors
    /// `FittingError::LengthMismatch` if `mean` or `covariance` does not
    /// match `parameters`; `FittingError::InvalidConfig` if `parameters` is
    /// empty or lists a parameter twice, `mean` is not finite, or `covariance`
    /// is not finite, not symmetric to 1e-12 of `√(Cᵢᵢ Cⱼⱼ)`, or not positive
    /// definite with the condition number of its correlation matrix at most
    /// 1e10; `FittingError::EvaluationFailed` if the eigendecomposition of
    /// that matrix fails.
    pub fn correlated(
        parameters: &[usize],
        mean: &[f64],
        covariance: &FlatMatrix,
    ) -> Result<Self, FittingError> {
        let k = parameters.len();
        if k == 0 {
            return Err(FittingError::InvalidConfig(
                "a prior needs at least one parameter".into(),
            ));
        }
        for (expected, actual, field) in [
            (k, mean.len(), "prior mean"),
            (k, covariance.nrows, "prior covariance rows"),
            (k, covariance.ncols, "prior covariance columns"),
            (k * k, covariance.data.len(), "prior covariance entries"),
        ] {
            if actual != expected {
                return Err(FittingError::LengthMismatch {
                    expected,
                    actual,
                    field,
                });
            }
        }
        if !mean.iter().all(|m| m.is_finite())
            || (0..k).any(|i| parameters[..i].contains(&parameters[i]))
        {
            return Err(FittingError::InvalidConfig(format!(
                "a prior needs distinct parameters and a finite mean; got {parameters:?} and \
                 {mean:?}"
            )));
        }
        let sd: Vec<f64> = (0..k).map(|i| covariance.get(i, i).sqrt()).collect();
        if !(covariance.data.iter().all(|c| c.is_finite())
            && sd.iter().all(|&s| s > 0.0)
            && symmetric(k, |i, j| covariance.get(i, j)))
        {
            return Err(FittingError::InvalidConfig(format!(
                "a prior covariance must be finite and symmetric with a positive diagonal; \
                 got {covariance:?}"
            )));
        }
        let correlation = Mat::from_fn(k, k, |i, j| {
            if i == j {
                1.0
            } else {
                covariance.get(i, j) / sd[i] / sd[j]
            }
        });
        let eigenvalues = correlation
            .self_adjoint_eigenvalues(Side::Lower)
            .map_err(|e| FittingError::EvaluationFailed(format!("{e:?}")))?;
        let (smallest, largest) = (eigenvalues[0], eigenvalues[k - 1]);
        if !(smallest > 0.0 && largest <= MAX_PRIOR_CONDITION * smallest) {
            return Err(FittingError::InvalidConfig(format!(
                "a prior covariance must be positive definite with a correlation condition \
                 number within {MAX_PRIOR_CONDITION:e}; its correlation eigenvalues are \
                 {eigenvalues:?}"
            )));
        }
        Ok(Self::factored(parameters, mean, covariance))
    }

    pub(crate) fn factored(parameters: &[usize], mean: &[f64], covariance: &FlatMatrix) -> Self {
        let k = parameters.len();
        let sd: Vec<f64> = (0..k).map(|i| covariance.get(i, i).sqrt()).collect();
        let correlation = |i: usize, j: usize| {
            if i == j {
                1.0
            } else {
                covariance.get(i, j) / sd[i] / sd[j]
            }
        };
        let mut factor = vec![0.0; k * k];
        for i in 0..k {
            for j in 0..=i {
                let dot: f64 = (0..j).map(|l| factor[i * k + l] * factor[j * k + l]).sum();
                factor[i * k + j] = if i == j {
                    (correlation(i, i) - dot).sqrt()
                } else {
                    (correlation(i, j) - dot) / factor[j * k + j]
                };
            }
        }
        for i in 0..k {
            for j in 0..=i {
                factor[i * k + j] *= sd[i];
            }
        }
        Self {
            parameters: parameters.to_vec(),
            mean: mean.to_vec(),
            factor,
        }
    }

    pub(crate) fn covariance(&self, i: usize, j: usize) -> f64 {
        let k = self.parameters.len();
        (0..=i.min(j))
            .map(|l| self.factor[i * k + l] * self.factor[j * k + l])
            .sum()
    }

    pub(crate) fn is_valid(&self) -> bool {
        let k = self.parameters.len();
        self.mean.iter().all(|m| m.is_finite())
            && self.factor.iter().all(|f| f.is_finite())
            && (0..k).all(|j| self.factor[j * k + j] > 0.0)
    }

    pub(crate) fn whiten(&self, v: &mut [f64]) {
        let k = self.parameters.len();
        for i in 0..k {
            let dot: f64 = (0..i).map(|j| self.factor[i * k + j] * v[j]).sum();
            v[i] = (v[i] - dot) / self.factor[i * k + i];
        }
    }

    fn residual(&self, values: &[f64]) -> Vec<f64> {
        let mut difference: Vec<f64> = self
            .parameters
            .iter()
            .zip(&self.mean)
            .map(|(&p, m)| values[p] - m)
            .collect();
        self.whiten(&mut difference);
        difference
    }

    fn inverse_factor(&self) -> Vec<f64> {
        let k = self.parameters.len();
        let mut inverse = vec![0.0; k * k];
        for j in 0..k {
            let mut column = vec![0.0; k];
            column[j] = 1.0;
            self.whiten(&mut column);
            for (i, c) in column.iter().enumerate() {
                inverse[i * k + j] = *c;
            }
        }
        inverse
    }
}

pub(crate) fn symmetric(k: usize, entry: impl Fn(usize, usize) -> f64) -> bool {
    (0..k).all(|i| {
        (0..i).all(|j| {
            let scale = (entry(i, i) * entry(j, j)).abs().sqrt();
            (entry(j, i) - entry(i, j)).abs() <= 1e-12 * scale
        })
    })
}

/// Result of [`poisson_fit`].
#[derive(Debug, Clone)]
pub struct PoissonResult {
    /// Half the Poisson deviance of the counts at `params`, without the
    /// priors.
    pub deviance: f64,
    /// Number of steps taken.
    pub iterations: usize,
    /// Whether the expected-information Newton decrement fell below 1e-6
    /// inside the bounds: within about 0.0014 standard errors of a minimum
    /// where the deviance is quadratic there; see [`poisson_fit`].
    pub converged: bool,
    /// Final parameter values (all parameters, including fixed).
    pub params: Vec<f64>,
    /// Covariance of the free parameters; the rows and columns of a
    /// parameter without an error bar are NaN.  `None` when the fit did not
    /// converge or covariance computation is disabled.
    pub covariance: Option<FlatMatrix>,
    /// Standard error of each free parameter; `None` for a parameter on a
    /// bound or with a component along a direction the data do not
    /// determine, larger than the SVD's rounding for that parameter.  `None`
    /// overall when `covariance` is.
    pub uncertainties: Option<Vec<Option<f64>>>,
    /// Whether each free parameter ended on one of its bounds.
    pub on_bound: Vec<bool>,
    /// Each count's leverage: the diagonal of the hat matrix of the free
    /// parameters off their bounds.  The leverages sum to the number of
    /// directions in those parameters the counts and priors determine, less
    /// the share the priors take.  `None` when `covariance` is `None` or its
    /// decomposition fails.
    pub leverage: Option<Vec<f64>>,
    /// The free parameters' Gaussian without their bounds; see [`Unbounded`].
    /// `None` when `covariance` is, when the decomposition over every free
    /// parameter fails, or when a bin predicted zero adds slope to the
    /// gradient of a parameter on a bound.
    pub unbounded: Option<Unbounded>,
}

/// The expected-information (Fisher scoring) expansion of the objective at
/// the fit over every free parameter, bounds left out: covariance `F⁻¹`, `F`
/// the expected information of all free parameters, and mean `θ̂ − F⁻¹g`, `g`
/// the objective's gradient at the fit, which the convergence test makes
/// negligible off the bounds.  Both are taken over the directions the counts
/// and priors determine; the mean and the rows and columns of a parameter
/// with a component along a direction they do not determine are NaN.  With
/// none on a bound the covariance is [`PoissonResult::covariance`].
#[derive(Debug, Clone)]
pub struct Unbounded {
    pub mean: Vec<f64>,
    pub covariance: FlatMatrix,
}

/// The expected-information Newton decrement below which [`poisson_fit`]
/// reports a minimum.
pub const NEWTON_DECREMENT_TOL: f64 = 1e-6;

const DEGENERATE_EIGENVALUE: f64 = 1e-12;

const MAX_PRIOR_CONDITION: f64 = 1e-2 / DEGENERATE_EIGENVALUE;

const MAX_REJECTIONS: usize = 60;

const INITIAL_DAMPING: f64 = 1e-3;

const DAMPING_FACTOR: f64 = 10.0;

/// `obs·ln(obs/mean) + mean − obs`, by C. Loader's `bd0` ("Fast and
/// accurate computation of binomial probabilities", 2000): a series in
/// `v = (obs − mean)/(obs + mean)` when `|obs − mean| < 0.1 (obs + mean)`.
/// For `obs ≥ 0` and `mean > 0`; 0 when both are 0, and +∞ when only `mean`
/// is.
pub fn half_deviance(obs: f64, mean: f64) -> f64 {
    let half_sum = obs / 2.0 + mean / 2.0;
    let difference = obs - mean;
    if difference.abs() < 0.2 * half_sum {
        let v = difference / 2.0 / half_sum;
        let mut sum = difference * v;
        let mut term = obs * (2.0 * v);
        let mut j = 1.0;
        loop {
            term *= v * v;
            let next = sum + term / (2.0 * j + 1.0);
            if next == sum {
                return sum;
            }
            sum = next;
            j += 1.0;
        }
    } else if obs == 0.0 {
        mean
    } else {
        obs * (obs.ln() - mean.ln()) + mean - obs
    }
}

fn deviance(y_obs: &[f64], y_model: &[f64]) -> f64 {
    y_obs
        .iter()
        .zip(y_model)
        .map(|(&obs, &mean)| {
            if mean > 0.0 || (mean == 0.0 && obs == 0.0) {
                half_deviance(obs, mean)
            } else {
                f64::INFINITY
            }
        })
        .sum()
}

fn objective(y_obs: &[f64], y_model: &[f64], priors: &[Prior], params: &ParameterSet) -> f64 {
    let values = params.all_values();
    deviance(y_obs, y_model)
        + priors
            .iter()
            .map(|p| 0.5 * p.residual(&values).iter().map(|w| w.powi(2)).sum::<f64>())
            .sum::<f64>()
}

struct Linearization {
    weighted: FlatMatrix,
    residual: Vec<f64>,
    zero_slope: Vec<f64>,
    gradient: Vec<f64>,
}

fn linearize(
    model: &dyn FitModel,
    params: &ParameterSet,
    free: &[usize],
    priors: &[Prior],
    y_obs: &[f64],
    y_model: &[f64],
) -> Result<Linearization, FittingError> {
    let mut weighted = model
        .analytical_jacobian(&params.all_values(), free, y_model)
        .ok_or_else(|| {
            FittingError::InvalidConfig(
                "poisson_fit needs a model with an analytical Jacobian".into(),
            )
        })?;
    for (expected, actual, field) in [
        (y_obs.len(), weighted.nrows, "analytical Jacobian rows"),
        (free.len(), weighted.ncols, "analytical Jacobian columns"),
    ] {
        if expected != actual {
            return Err(FittingError::LengthMismatch {
                expected,
                actual,
                field,
            });
        }
    }
    let root: Vec<f64> = y_model.iter().map(|mean| mean.sqrt()).collect();
    let mut zero_slope = vec![0.0; weighted.ncols];
    for (i, &r) in root.iter().enumerate() {
        for (j, z) in zero_slope.iter_mut().enumerate() {
            let slope = weighted.get(i, j);
            if r > 0.0 {
                *weighted.get_mut(i, j) = slope / r;
            } else {
                *z += slope;
                *weighted.get_mut(i, j) = 0.0;
            }
        }
    }
    let mut residual: Vec<f64> = y_obs
        .iter()
        .zip(y_model)
        .zip(&root)
        .map(|((&obs, &mean), &r)| if r > 0.0 { (mean - obs) / r } else { 0.0 })
        .collect();
    let mut row = weighted.nrows;
    weighted.nrows += priors.iter().map(|p| p.parameters.len()).sum::<usize>();
    weighted.data.resize(weighted.nrows * weighted.ncols, 0.0);
    let values = params.all_values();
    for prior in priors {
        let k = prior.parameters.len();
        let inverse = prior.inverse_factor();
        for (i, w) in prior.residual(&values).into_iter().enumerate() {
            for (j, parameter) in prior.parameters.iter().enumerate().take(i + 1) {
                let col = free
                    .iter()
                    .position(|index| index == parameter)
                    .expect("priors are on free parameters");
                *weighted.get_mut(row, col) = inverse[i * k + j];
            }
            residual.push(w);
            row += 1;
        }
    }
    let gradient = zero_slope
        .iter()
        .enumerate()
        .map(|(j, z)| {
            z + (0..weighted.nrows)
                .map(|i| weighted.get(i, j) * residual[i])
                .sum::<f64>()
        })
        .collect();
    Ok(Linearization {
        weighted,
        residual,
        zero_slope,
        gradient,
    })
}

fn column_scale(weighted: &FlatMatrix, col: usize) -> (f64, f64) {
    let peak = (0..weighted.nrows).fold(0.0_f64, |m, i| m.max(weighted.get(i, col).abs()));
    let length = (0..weighted.nrows)
        .map(|i| (weighted.get(i, col) / peak).powi(2))
        .sum::<f64>()
        .sqrt();
    (peak, length)
}

struct Decomposition {
    columns: Vec<usize>,
    largest: Vec<f64>,
    length: Vec<f64>,
    rows: usize,
    singular: Vec<f64>,
    left: Mat<f64>,
    right: Vec<Vec<f64>>,
}

impl Decomposition {
    fn new(weighted: &FlatMatrix, columns: &[usize]) -> Option<Self> {
        let (mut kept, mut largest, mut length) = (vec![], vec![], vec![]);
        for &col in columns {
            let (peak, norm) = column_scale(weighted, col);
            if peak > 0.0 {
                kept.push(col);
                largest.push(peak);
                length.push(norm);
            }
        }
        let (rows, k) = (weighted.nrows, kept.len());
        let scaled = Mat::from_fn(rows, k, |row, i| {
            weighted.get(row, kept[i]) / largest[i] / length[i]
        });
        let mut values = Mat::<f64>::zeros(rows.min(k), rows.min(k));
        let mut left = Mat::<f64>::zeros(rows, rows.min(k));
        let mut v = Mat::<f64>::zeros(k, k);
        let params = Default::default();
        svd(
            scaled.as_ref(),
            values.diagonal_mut(),
            Some(left.as_mut()),
            Some(v.as_mut()),
            Par::Seq,
            MemStack::new(&mut MemBuffer::new(svd_scratch::<f64>(
                rows,
                k,
                ComputeSvdVectors::Thin,
                ComputeSvdVectors::Full,
                Par::Seq,
                params,
            ))),
            params,
        )
        .ok()?;
        Some(Self {
            columns: kept,
            largest,
            length,
            rows,
            singular: (0..k)
                .map(|j| if j < rows { values[(j, j)] } else { 0.0 })
                .collect(),
            left,
            right: (0..k)
                .map(|j| (0..k).map(|i| v[(i, j)]).collect())
                .collect(),
        })
    }

    fn unscale(&self, i: usize, value: f64) -> f64 {
        value / self.length[i] / self.largest[i]
    }

    fn spanned(&self) -> impl Iterator<Item = usize> + '_ {
        let largest = self.singular.iter().fold(0.0_f64, |m, &s| m.max(s));
        let rank_floor = f64::EPSILON * self.rows.max(self.singular.len()) as f64 * largest;
        (0..self.singular.len()).filter(move |&k| self.singular[k] > rank_floor)
    }

    fn along(&self, k: usize, linear: &Linearization) -> f64 {
        let zero_part: f64 = self
            .columns
            .iter()
            .enumerate()
            .map(|(i, &col)| self.right[k][i] * self.unscale(i, linear.zero_slope[col]))
            .sum();
        let residual_part: f64 = (0..self.rows)
            .map(|row| self.left[(row, k)] * linear.residual[row])
            .sum();
        residual_part + zero_part / self.singular[k]
    }

    fn newton_decrement(&self, linear: &Linearization) -> f64 {
        0.5 * self
            .spanned()
            .map(|k| self.along(k, linear).powi(2))
            .sum::<f64>()
    }

    fn is_determined(&self, k: usize) -> bool {
        self.singular[k].powi(2) >= DEGENERATE_EIGENVALUE
    }

    fn determined(&self) -> impl Iterator<Item = usize> + '_ {
        (0..self.singular.len()).filter(|&k| self.is_determined(k))
    }

    fn step(
        &self,
        directions: impl Iterator<Item = usize>,
        linear: &Linearization,
        n_free: usize,
        damping: f64,
    ) -> Vec<f64> {
        let mut direction = vec![0.0; n_free];
        for k in directions {
            let s = self.singular[k];
            let coefficient = self.along(k, linear) * s / (s * s + damping);
            for (i, &col) in self.columns.iter().enumerate() {
                direction[col] += self.unscale(i, self.right[k][i] * coefficient);
            }
        }
        direction
    }

    fn leverage(&self, counts: usize) -> Vec<f64> {
        (0..counts)
            .map(|row| self.spanned().map(|k| self.left[(row, k)].powi(2)).sum())
            .collect()
    }

    fn error_bars(&self, n_free: usize) -> (FlatMatrix, Vec<Option<f64>>) {
        let n = self.columns.len();
        let determined: Vec<bool> = (0..n).map(|k| self.is_determined(k)).collect();
        let largest = self.singular.iter().fold(0.0_f64, |m, &s| m.max(s));
        let resolved: Vec<bool> = (0..n)
            .map(|i| {
                let sensitivity: f64 = (0..n)
                    .filter(|&k| determined[k])
                    .map(|k| self.right[k][i].abs() / self.singular[k])
                    .sum();
                let rounding = f64::EPSILON * self.rows.max(n) as f64 * largest * sensitivity;
                (0..n)
                    .filter(|&k| !determined[k])
                    .map(|k| self.right[k][i].powi(2))
                    .sum::<f64>()
                    <= rounding.powi(2)
            })
            .collect();
        let variance = |i: usize, j: usize| -> f64 {
            (0..n)
                .filter(|&k| determined[k])
                .map(|k| self.right[k][i] * self.right[k][j] / self.singular[k].powi(2))
                .sum()
        };
        let (mut covariance, mut errors) = withheld(n_free);
        for i in (0..n).filter(|&i| resolved[i]) {
            for j in (0..n).filter(|&j| resolved[j]) {
                *covariance.get_mut(self.columns[i], self.columns[j]) =
                    self.unscale(j, self.unscale(i, variance(i, j)));
            }
            errors[self.columns[i]] = Some(self.unscale(i, variance(i, i).sqrt()));
        }
        (covariance, errors)
    }
}

/// The inverses of the expected information `information` of `rows` counts,
/// split as [`poisson_fit`] splits its directions into those its steps take
/// and those its covariance keeps.  In the information scaled to a unit
/// `diagonal`, its own or, for an information with other parameters profiled
/// out, the one it had before, `spanned` keeps every eigen-direction whose
/// eigenvalue is above the eigendecomposition's rounding, `ε·max(rows, k)`
/// times the largest, `k` the number of parameters with information, and `determined` keeps those of them whose eigenvalue is
/// also at least 1e-12; where that rounding exceeds 1e-12, it is the cutoff
/// of both.  A parameter with no information has a zero row and
/// column in both.  `resolved[i]` is whether parameter `i` has information
/// and, beyond rounding, no component along a direction `determined` leaves
/// out; `undetermined` holds those directions, in the parameters' units.
#[derive(Debug, Clone)]
pub struct InformationInverse {
    pub spanned: FlatMatrix,
    pub determined: FlatMatrix,
    pub resolved: Vec<bool>,
    pub undetermined: Vec<Vec<f64>>,
}

/// [`InformationInverse`] of `information`, the expected information of
/// `rows` counts, scaled to `diagonal`.
///
/// # Errors
/// [`FittingError::InvalidConfig`] if `information` is not square, does not
/// hold its entries, has an entry that is not finite, or is not symmetric to
/// 1e-12 of `√(dᵢdⱼ)`, or if `diagonal` is not one finite non-negative entry
/// per row, zero only on a row of zeros; [`FittingError::EvaluationFailed`] if its eigendecomposition
/// fails.
pub fn information_inverse(
    information: &FlatMatrix,
    diagonal: &[f64],
    rows: usize,
) -> Result<InformationInverse, FittingError> {
    let n = information.nrows;
    if information.ncols != n
        || information.data.len() != n * n
        || !information.data.iter().all(|x| x.is_finite())
        || diagonal.len() != n
        || !diagonal.iter().all(|d| d.is_finite() && *d >= 0.0)
        || (0..n).any(|i| diagonal[i] == 0.0 && (0..n).any(|j| information.get(i, j) != 0.0))
        || !(0..n).all(|i| {
            (0..i).all(|j| {
                (information.get(i, j) - information.get(j, i)).abs()
                    <= 1e-12 * (diagonal[i] * diagonal[j]).sqrt()
            })
        })
    {
        return Err(FittingError::InvalidConfig(format!(
            "an expected information is a finite symmetric square matrix with a finite \
             non-negative diagonal to scale it by; got {} by {} with {} entries and {} \
             diagonal entries",
            information.nrows,
            information.ncols,
            information.data.len(),
            diagonal.len()
        )));
    }
    let kept: Vec<usize> = (0..n).filter(|&i| diagonal[i] > 0.0).collect();
    let k = kept.len();
    if k == 0 {
        return Ok(InformationInverse {
            spanned: FlatMatrix::zeros(n, n),
            determined: FlatMatrix::zeros(n, n),
            resolved: vec![false; n],
            undetermined: Vec::new(),
        });
    }
    let scale: Vec<f64> = kept.iter().map(|&i| diagonal[i].sqrt()).collect();
    let eigen = Mat::from_fn(k, k, |a, b| {
        information.get(kept[a], kept[b]) / (scale[a] * scale[b])
    })
    .self_adjoint_eigen(Side::Lower)
    .map_err(|e| FittingError::EvaluationFailed(format!("{e:?}")))?;
    let (values, vectors) = (eigen.S().column_vector(), eigen.U());
    let singular: Vec<f64> = (0..k).map(|d| values[d].max(0.0).sqrt()).collect();
    let largest = singular.iter().fold(0.0_f64, |m, &s| m.max(s));
    let rank_floor = f64::EPSILON * rows.max(k) as f64 * largest.powi(2);
    let inverse = |keep: &dyn Fn(usize) -> bool| {
        let mut inverse = FlatMatrix::zeros(n, n);
        for a in 0..k {
            for b in 0..k {
                *inverse.get_mut(kept[a], kept[b]) = (0..k)
                    .filter(|&d| keep(d))
                    .map(|d| vectors[(a, d)] * vectors[(b, d)] / values[d])
                    .sum::<f64>()
                    / (scale[a] * scale[b]);
            }
        }
        inverse
    };
    let is_spanned = |d: usize| values[d] > rank_floor;
    let is_determined = |d: usize| is_spanned(d) && values[d] >= DEGENERATE_EIGENVALUE;
    let mut resolved = vec![false; n];
    for a in 0..k {
        let sensitivity: f64 = (0..k)
            .filter(|&d| is_determined(d))
            .map(|d| vectors[(a, d)].abs() / singular[d])
            .sum();
        let rounding = f64::EPSILON * rows.max(k) as f64 * largest * sensitivity;
        resolved[kept[a]] = (0..k)
            .filter(|&d| !is_determined(d))
            .map(|d| vectors[(a, d)].powi(2))
            .sum::<f64>()
            <= rounding.powi(2);
    }
    let undetermined = (0..k)
        .filter(|&d| !is_determined(d))
        .map(|d| {
            let mut direction = vec![0.0; n];
            for a in 0..k {
                direction[kept[a]] = vectors[(a, d)] / scale[a];
            }
            direction
        })
        .collect();
    Ok(InformationInverse {
        spanned: inverse(&is_spanned),
        determined: inverse(&is_determined),
        resolved,
        undetermined,
    })
}

fn withheld(n_free: usize) -> (FlatMatrix, Vec<Option<f64>>) {
    let mut covariance = FlatMatrix::zeros(n_free, n_free);
    covariance.data.fill(f64::NAN);
    (covariance, vec![None; n_free])
}

fn on_bound(param: &FitParameter) -> bool {
    param.value == param.lower || param.value == param.upper
}

fn held_by_bound(param: &FitParameter, gradient: f64) -> bool {
    (param.value == param.lower && gradient > 0.0) || (param.value == param.upper && gradient < 0.0)
}

/// Fit `params` to the counts `y_obs` by minimizing half the Poisson
/// deviance plus `½(θₚ − mₚ)ᵀCₚ⁻¹(θₚ − mₚ)` for each of the `priors`, within
/// the parameter bounds.
///
/// The model must provide an analytical Jacobian at every point the fit
/// visits.  A trial predicting a negative count, or zero where something was
/// counted, is rejected; a bin predicted zero with no counts adds its slope
/// to the gradient and nothing to the information.  Each step is the
/// Levenberg–Marquardt step `(F + λI)⁻¹g` in coordinates scaled to unit
/// Fisher information, projected onto the box;
/// `λ` is divided by 10 after a step that lowers the objective, and
/// multiplied by 10, to at least its starting 1e-3, before retrying one that
/// does not (D. W. Marquardt, J. Soc. Indust. Appl. Math. 11, 431–441,
/// 1963).  A free parameter on its bound whose gradient points out of the
/// box is held there for the step and left out of the convergence test.
///
/// The fit has converged when the Newton decrement `½ gᵀF⁺g` over the
/// parameters not held by a bound (on it, gradient pointing out),
/// `F = Jᵀ diag(1/μ) J + Σₚ Cₚ⁻¹` the expected information, each `Cₚ⁻¹` on
/// its prior's parameters, is below 1e-6: the
/// quadratic model built from `F` predicts less than 1e-6 of further
/// decrease.  Where the deviance is quadratic near the minimum this puts the
/// fit within about 0.0014 standard errors of it, and within 0.01 where the
/// expected information overstates the curvature, as at one count per bin.
/// Where it is not quadratic — weakly determined parameters on a curved or
/// flat likelihood ridge, a minimum at infinity, or a parameter that adds
/// counts to a bin predicted near zero — the fit can stop short of the
/// minimum, or at another local minimum, and the error bars are not standard
/// errors.  A minimum where a bin's prediction reaches zero inside the box is
/// reported unconverged.
///
/// It stops unconverged when no step lowers the objective, when the Jacobian
/// is not finite, when the start predicts a negative count or zero where
/// something was counted, or after `config.max_iter` steps.  A trial whose
/// prediction has a different length from `y_obs` is rejected.
///
/// # Errors
/// `FittingError::EmptyData` if `y_obs` is empty;
/// `FittingError::InvalidConfig` if an observation is negative or not
/// finite, a parameter value is not finite, a free parameter's bounds are
/// inverted, NaN or admit no finite value, a prior is not on free
/// parameters, two priors are on one parameter, a prior's mean is not finite
/// or a measured sd not finite and positive, or the model returns no
/// analytical Jacobian;
/// `FittingError::LengthMismatch` if the prediction or the Jacobian does not
/// match `y_obs` and the free parameters; the model's error if it fails at
/// the start.
pub fn poisson_fit(
    model: &dyn FitModel,
    y_obs: &[f64],
    priors: &[Prior],
    params: &mut ParameterSet,
    config: &PoissonConfig,
) -> Result<PoissonResult, FittingError> {
    if let Some((bin, &obs)) = y_obs
        .iter()
        .enumerate()
        .find(|(_, obs)| !(obs.is_finite() && **obs >= 0.0))
    {
        return Err(FittingError::InvalidConfig(format!(
            "observed counts must be finite and non-negative, got {obs} in bin {bin}"
        )));
    }
    if y_obs.is_empty() {
        return Err(FittingError::EmptyData);
    }
    if let Some(p) = params.params.iter().find(|p| {
        !p.value.is_finite()
            || (!p.fixed
                && (p.lower.is_nan()
                    || p.upper.is_nan()
                    || p.lower > p.upper
                    || p.lower == f64::INFINITY
                    || p.upper == f64::NEG_INFINITY))
    }) {
        return Err(FittingError::InvalidConfig(format!(
            "parameter {} = {} with bounds [{}, {}]",
            p.name, p.value, p.lower, p.upper
        )));
    }
    let covered: Vec<usize> = priors.iter().flat_map(|p| p.parameters.clone()).collect();
    if let Some((i, prior)) = priors.iter().enumerate().find(|(_, p)| {
        p.parameters.iter().any(|&parameter| {
            params.params.get(parameter).is_none_or(|param| param.fixed)
                || covered.iter().filter(|&&c| c == parameter).count() > 1
        }) || !p.is_valid()
    }) {
        return Err(FittingError::InvalidConfig(format!(
            "prior {i} must be on free parameters without another prior, with a finite mean \
             and finite positive sds; got parameters {:?} with mean {:?}",
            prior.parameters, prior.mean
        )));
    }
    params.set_free_values(&params.free_values());
    let mut y_model = model.evaluate(&params.all_values())?;
    if y_model.len() != y_obs.len() {
        return Err(FittingError::LengthMismatch {
            expected: y_model.len(),
            actual: y_obs.len(),
            field: "y_obs",
        });
    }
    let free = params.free_indices();
    let mut value = objective(y_obs, &y_model, priors, params);
    let mut iterations = 0;
    let mut at_minimum = None;
    let mut damping = INITIAL_DAMPING;
    while value.is_finite() {
        let linear = linearize(model, params, &free, priors, y_obs, &y_model)?;
        if !linear
            .weighted
            .data
            .iter()
            .chain(&linear.zero_slope)
            .all(|v| v.is_finite())
        {
            break;
        }
        let movable: Vec<usize> = (0..free.len())
            .filter(|&j| !held_by_bound(&params.params[free[j]], linear.gradient[j]))
            .collect();
        let Some(decomposition) = Decomposition::new(&linear.weighted, &movable) else {
            break;
        };
        if decomposition.newton_decrement(&linear) < NEWTON_DECREMENT_TOL {
            at_minimum = Some(linear);
            break;
        }
        if iterations == config.max_iter {
            break;
        }
        let start = params.free_values();
        let mut accepted = None;
        for _ in 0..MAX_REJECTIONS {
            let direction =
                decomposition.step(decomposition.spanned(), &linear, free.len(), damping);
            let unprojected: Vec<f64> = start.iter().zip(&direction).map(|(x, d)| x - d).collect();
            params.set_free_values(&unprojected);
            if let Some(trial_model) = model
                .evaluate(&params.all_values())
                .ok()
                .filter(|trial| trial.len() == y_obs.len())
            {
                let trial_value = objective(y_obs, &trial_model, priors, params);
                if trial_value < value {
                    damping /= DAMPING_FACTOR;
                    accepted = Some((trial_model, trial_value));
                    break;
                }
            }
            params.set_free_values(&start);
            damping = (damping * DAMPING_FACTOR).max(INITIAL_DAMPING);
        }
        let Some((trial_model, trial_value)) = accepted else {
            break;
        };
        y_model = trial_model;
        value = trial_value;
        iterations += 1;
    }

    let bounded: Vec<bool> = free
        .iter()
        .map(|&idx| on_bound(&params.params[idx]))
        .collect();
    let (covariance, uncertainties, leverage, unbounded) = match &at_minimum {
        Some(linear) if config.compute_covariance => {
            let interior: Vec<usize> = (0..free.len()).filter(|&j| !bounded[j]).collect();
            let decomposition = Decomposition::new(&linear.weighted, &interior);
            let (covariance, errors) = decomposition.as_ref().map_or_else(
                || withheld(free.len()),
                |decomposition| decomposition.error_bars(free.len()),
            );
            let leverage = decomposition
                .as_ref()
                .map(|decomposition| decomposition.leverage(y_obs.len()));
            let unbounded = unbounded(linear, &params.free_values(), &bounded);
            (Some(covariance), Some(errors), leverage, unbounded)
        }
        _ => (None, None, None, None),
    };
    Ok(PoissonResult {
        deviance: deviance(y_obs, &y_model),
        iterations,
        converged: at_minimum.is_some(),
        params: params.all_values(),
        covariance,
        uncertainties,
        on_bound: bounded,
        leverage,
        unbounded,
    })
}

fn unbounded(linear: &Linearization, fitted: &[f64], bounded: &[bool]) -> Option<Unbounded> {
    if (0..fitted.len()).any(|j| bounded[j] && linear.zero_slope[j] != 0.0) {
        return None;
    }
    let columns: Vec<usize> = (0..fitted.len()).collect();
    let decomposition = Decomposition::new(&linear.weighted, &columns)?;
    let step = decomposition.step(decomposition.determined(), linear, fitted.len(), 0.0);
    let covariance = decomposition.error_bars(fitted.len()).0;
    let mean = (0..fitted.len())
        .map(|i| {
            if covariance.get(i, i).is_nan() {
                f64::NAN
            } else {
                fitted[i] - step[i]
            }
        })
        .collect();
    Some(Unbounded { mean, covariance })
}

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

    /// Simple model: y = a * exp(-b * x)
    /// This mimics transmission: counts = flux * exp(-density * sigma)
    struct ExponentialModel {
        x: Vec<f64>,
        flux: Vec<f64>,
    }

    impl FitModel for ExponentialModel {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            let b = params[0]; // "density"
            Ok(self
                .x
                .iter()
                .zip(self.flux.iter())
                .map(|(&xi, &fi)| fi * (-b * xi).exp())
                .collect())
        }
    }

    struct Decay {
        t: Vec<f64>,
        jacobian_factor: f64,
    }

    impl FitModel for Decay {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(self
                .t
                .iter()
                .map(|&t| params[0] * (-params[1] * t).exp())
                .collect())
        }

        fn analytical_jacobian(
            &self,
            params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(self.t.len(), free_param_indices.len());
            for (row, &t) in self.t.iter().enumerate() {
                let e = (-params[1] * t).exp();
                for (col, &index) in free_param_indices.iter().enumerate() {
                    let slope = [e, -params[0] * t * e][index];
                    *jacobian.get_mut(row, col) = slope * self.jacobian_factor;
                }
            }
            Some(jacobian)
        }
    }

    fn decay(jacobian_factor: f64) -> Decay {
        Decay {
            t: (0..25).map(|i| 0.2 * f64::from(i)).collect(),
            jacobian_factor,
        }
    }

    fn decay_params(a: f64, b: f64, a_upper: f64) -> ParameterSet {
        ParameterSet::new(vec![
            FitParameter {
                name: "a".into(),
                value: a,
                lower: 0.0,
                upper: a_upper,
                fixed: false,
            },
            FitParameter::non_negative("b", b),
        ])
    }

    fn fit_decay(
        model: &Decay,
        observed: &[f64],
        start: (f64, f64),
        max_iter: usize,
    ) -> PoissonResult {
        let mut params = decay_params(start.0, start.1, f64::INFINITY);
        let config = PoissonConfig {
            max_iter,
            ..PoissonConfig::default()
        };
        poisson_fit(model, observed, &[], &mut params, &config).unwrap()
    }

    #[test]
    fn half_deviance_matches_high_precision_values() {
        for (obs, mean, exact) in [
            (1_000_001.0, 1_000_000.0, 4.999_998_333_334_167e-7),
            (3.0, 2.5, 0.046_964_670_381_863_88),
            (0.0, 2.5, 2.5),
            (1.0e308, 9.0e307, 5.360_515_657_826_301e305),
            (1.7e308, 1.65e308, 7.500_373_544_579_622e304),
            (1.0, 1.0e-310, 712.801_378_828_154_2),
        ] {
            let got = half_deviance(obs, mean);
            assert!(
                ((got - exact) / exact).abs() < 1e-12,
                "{obs}, {mean}: {got}"
            );
        }
    }

    #[test]
    fn inputs_that_are_not_counts_for_this_model_are_refused() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        let fit = |model: &dyn FitModel, observed: &[f64]| {
            poisson_fit(
                model,
                observed,
                &[],
                &mut decay_params(50.0, 1.0, f64::INFINITY),
                &PoissonConfig::default(),
            )
        };
        let mut bad = observed.clone();
        bad[3] = f64::NAN;
        assert!(fit(&model, &bad).is_err());
        bad[3] = -1.0;
        assert!(fit(&model, &bad).is_err());
        assert!(fit(&model, &observed[1..]).is_err());
        let empty = Decay {
            t: vec![],
            jacobian_factor: 1.0,
        };
        assert!(fit(&empty, &[]).is_err());
        let no_jacobian = ExponentialModel {
            x: model.t.clone(),
            flux: vec![100.0; model.t.len()],
        };
        assert!(fit(&no_jacobian, &observed).is_err());
    }

    #[test]
    fn a_start_with_a_zero_prediction_is_not_converged() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        let result = fit_decay(&model, &observed, (0.0, 1.0), 200);
        assert!(!result.converged && result.iterations == 0, "{result:?}");
    }

    #[test]
    fn a_non_finite_slope_ends_the_fit_unconverged() {
        let model = decay(f64::NAN);
        let observed = decay(1.0).evaluate(&[100.0, 1.5]).unwrap();
        let result = fit_decay(&model, &observed, (50.0, 1.0), 200);
        assert!(
            !result.converged && result.uncertainties.is_none(),
            "{result:?}"
        );
    }

    #[test]
    fn a_fit_whose_step_cannot_lower_the_deviance_is_not_converged() {
        let model = decay(-1.0);
        let observed = decay(1.0).evaluate(&[100.0, 1.5]).unwrap();
        let result = fit_decay(&model, &observed, (50.0, 1.0), 200);
        assert!(
            !result.converged && result.iterations == 0 && result.uncertainties.is_none(),
            "{result:?}"
        );
    }

    #[test]
    fn a_start_outside_the_bounds_ends_inside_them() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        let mut params = decay_params(100.0, 1.5, 80.0);
        let result = poisson_fit(
            &model,
            &observed,
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(result.converged && result.params[0] == 80.0, "{result:?}");
        assert_eq!(result.on_bound, vec![true, false]);
    }

    #[test]
    fn a_fit_that_reaches_the_minimum_on_its_last_allowed_step_has_converged() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        let steps = fit_decay(&model, &observed, (20.0, 0.3), 200).iterations;
        assert!(steps > 1);
        assert!(fit_decay(&model, &observed, (20.0, 0.3), steps).converged);
        assert!(!fit_decay(&model, &observed, (20.0, 0.3), steps - 1).converged);
    }

    struct Scaled {
        x: Vec<f64>,
        slope: f64,
        tiny: f64,
    }

    impl FitModel for Scaled {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(self
                .x
                .iter()
                .map(|&x| 100.0 * (self.slope * params[0] * x).exp() + self.tiny * params[1] * x)
                .collect())
        }

        fn analytical_jacobian(
            &self,
            params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(self.x.len(), free_param_indices.len());
            for (row, &x) in self.x.iter().enumerate() {
                let slopes = [
                    100.0 * self.slope * x * (self.slope * params[0] * x).exp(),
                    self.tiny * x,
                ];
                for (col, &index) in free_param_indices.iter().enumerate() {
                    *jacobian.get_mut(row, col) = slopes[index];
                }
            }
            Some(jacobian)
        }
    }

    fn fit_scaled(slope: f64, tiny: f64) -> PoissonResult {
        let model = Scaled {
            x: vec![1.0, 2.0, 3.0],
            slope,
            tiny,
        };
        let observed = model.evaluate(&[0.3 / slope, 0.0]).unwrap();
        let mut params = ParameterSet::new(vec![
            FitParameter::unbounded("theta", 0.0),
            FitParameter::unbounded("weak", 0.0),
        ]);
        poisson_fit(
            &model,
            &observed,
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap()
    }

    #[test]
    fn a_slope_too_large_to_square_is_fitted() {
        let result = fit_scaled(1.0e160, 1.0);
        assert!(
            result.converged && result.deviance < NEWTON_DECREMENT_TOL,
            "{result:?}"
        );
    }

    #[test]
    fn a_slope_too_small_to_square_does_not_fake_convergence() {
        let result = fit_scaled(1.0, 1.0e-310);
        assert!(
            !result.converged || result.deviance < NEWTON_DECREMENT_TOL,
            "{result:?}"
        );
    }

    struct Line {
        x: Vec<f64>,
    }

    impl FitModel for Line {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(self.x.iter().map(|&x| 1.0 + params[0] * x).collect())
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(self.x.len(), free_param_indices.len());
            jacobian.data.copy_from_slice(&self.x);
            Some(jacobian)
        }
    }

    struct Proportional {
        a: Vec<f64>,
    }

    impl FitModel for Proportional {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(self.a.iter().map(|a| a * params[0]).collect())
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(self.a.len(), free_param_indices.len());
            jacobian.data.copy_from_slice(&self.a);
            Some(jacobian)
        }
    }

    #[test]
    fn a_prior_moves_the_fit_to_the_maximum_of_the_posterior() {
        let model = Proportional {
            a: vec![1.0, 2.0, 3.0, 4.0],
        };
        let observed = [3.0, 5.0, 9.0, 12.0];
        let (a_sum, y_sum): (f64, f64) = (10.0, 29.0);
        let (mean, sd) = (2.0, 0.3);
        let fit = |priors: &[Prior]| {
            let mut params = ParameterSet::new(vec![FitParameter::non_negative("theta", 1.0)]);
            poisson_fit(
                &model,
                &observed,
                priors,
                &mut params,
                &PoissonConfig::default(),
            )
            .unwrap()
        };
        let result = fit(&[Prior::measured(0, mean, sd)]);
        assert!(result.converged);
        let theta = result.params[0];
        let b = a_sum - mean / sd.powi(2);
        let root = sd.powi(2) * 0.5 * (-b + (b * b + 4.0 * y_sum / sd.powi(2)).sqrt());
        let variance = 1.0 / (a_sum / theta + 1.0 / sd.powi(2));
        assert!(
            (theta - root).abs() <= 2e-3 * variance.sqrt(),
            "{theta} vs {root}"
        );
        let covariance = result.covariance.expect("covariance").get(0, 0);
        assert!((covariance / variance - 1.0).abs() <= 1e-12, "{covariance}");
        let deviance: f64 = observed
            .iter()
            .zip(&model.a)
            .map(|(y, a)| a * theta - y + y * (y / (a * theta)).ln())
            .sum();
        assert!((result.deviance / deviance - 1.0).abs() <= 1e-12);
        let leverage = result.leverage.expect("leverage");
        for (h, a) in leverage.iter().zip(&model.a) {
            assert!((h / (a / theta * variance) - 1.0).abs() <= 1e-12, "{h}");
        }
        let unmeasured = fit(&[]).leverage.expect("leverage");
        assert!((unmeasured.iter().sum::<f64>() - 1.0).abs() <= 1e-12);
    }

    #[test]
    fn priors_that_are_not_measurements_of_one_free_parameter_are_refused() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        let fit = |priors: &[Prior], fixed: bool| {
            let mut params = decay_params(50.0, 1.0, f64::INFINITY);
            params.params[1].fixed = fixed;
            poisson_fit(
                &model,
                &observed,
                priors,
                &mut params,
                &PoissonConfig::default(),
            )
        };
        let prior = Prior::measured;
        assert!(fit(&[prior(1, 1.5, 0.1)], false).is_ok());
        for priors in [
            vec![prior(2, 1.5, 0.1)],
            vec![prior(1, 1.5, 0.1), prior(1, 1.4, 0.1)],
            vec![prior(1, f64::NAN, 0.1)],
            vec![prior(1, 1.5, 0.0)],
            vec![prior(1, 1.5, f64::INFINITY)],
            vec![block(&[0, 1]), prior(1, 1.5, 0.1)],
        ] {
            assert!(fit(&priors, false).is_err(), "{priors:?}");
        }
        assert!(fit(&[prior(1, 1.5, 0.1)], true).is_err());
        assert!(fit(&[block(&[0, 1])], false).is_ok());
    }

    fn two_by_two(data: [f64; 4]) -> FlatMatrix {
        FlatMatrix {
            data: data.to_vec(),
            nrows: 2,
            ncols: 2,
        }
    }

    fn block(parameters: &[usize]) -> Prior {
        Prior::correlated(
            parameters,
            &[100.0, 1.5],
            &two_by_two([1.0, 0.0, 0.0, 0.01]),
        )
        .unwrap()
    }

    #[test]
    fn covariances_that_are_not_symmetric_and_well_conditioned_are_refused() {
        let correlated = |data| Prior::correlated(&[0, 1], &[0.0, 0.0], &two_by_two(data));
        let scaled = |rho: f64| correlated([1.0, rho * 1e3, rho * 1e3, 1e6]);
        assert!(scaled(1.0 - 1e-9).is_ok());
        assert!(correlated([1.0, 1e-6, 1e-6 + 1e-17, 1.0]).is_ok());
        for result in [
            scaled(1.0 - 1e-11),
            scaled(1.0 + 1e-9),
            correlated([1.0, 0.5, 0.5 + 1e-10, 1.0]),
            correlated([1.0, 1e-6, 1e-6 + 1e-11, 1.0]),
            correlated([1.0, f64::NAN, f64::NAN, 1.0]),
            correlated([0.0, 0.0, 0.0, 1.0]),
            Prior::correlated(&[0, 1], &[0.0], &two_by_two([1.0, 0.0, 0.0, 1.0])),
            Prior::correlated(
                &[0, 1],
                &[0.0, 0.0],
                &FlatMatrix {
                    data: vec![1.0, 0.0, 0.0],
                    nrows: 2,
                    ncols: 2,
                },
            ),
            Prior::correlated(&[], &[], &FlatMatrix::zeros(0, 0)),
            Prior::correlated(&[0, 1], &[f64::NAN, 0.0], &two_by_two([1.0, 0.0, 0.0, 1.0])),
            Prior::correlated(&[1, 1], &[0.0, 0.0], &two_by_two([1.0, 0.0, 0.0, 1.0])),
        ] {
            assert!(result.is_err(), "{result:?}");
        }
    }

    #[test]
    fn a_measured_prior_divides_by_its_sd() {
        for (value, mean, sd) in [
            (1.234_567, 0.987_654_3, 0.031_4),
            (1e-7, -3.3e5, 7.0e3),
            (2.0, 2.0, 0.3),
        ] {
            let prior = Prior::measured(0, mean, sd);
            assert_eq!(
                prior.residual(&[value])[0].to_bits(),
                ((value - mean) / sd).to_bits()
            );
            assert_eq!(prior.inverse_factor()[0].to_bits(), (1.0 / sd).to_bits());
        }
    }

    struct Mixture;

    const MIXTURE: [[f64; 2]; 5] = [[1.0, 5.0], [2.0, 3.0], [3.0, 2.0], [4.0, 1.0], [5.0, 0.5]];

    const MIXTURE_COUNTS: [f64; 5] = [21.0, 20.0, 15.0, 18.0, 16.0];

    impl FitModel for Mixture {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(MIXTURE
                .iter()
                .map(|a| a[0] * params[0] + a[1] * params[1])
                .collect())
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(MIXTURE.len(), free_param_indices.len());
            for (row, a) in MIXTURE.iter().enumerate() {
                for (col, &index) in free_param_indices.iter().enumerate() {
                    *jacobian.get_mut(row, col) = a[index];
                }
            }
            Some(jacobian)
        }
    }

    fn fit_mixture(priors: &[Prior]) -> PoissonResult {
        let mut params = ParameterSet::new(vec![
            FitParameter::non_negative("x", 1.0),
            FitParameter::non_negative("y", 1.0),
        ]);
        poisson_fit(
            &Mixture,
            &MIXTURE_COUNTS,
            priors,
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap()
    }

    fn inverse(m: [[f64; 2]; 2]) -> [[f64; 2]; 2] {
        let det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
        [
            [m[1][1] / det, -m[0][1] / det],
            [-m[1][0] / det, m[0][0] / det],
        ]
    }

    #[test]
    fn a_correlated_prior_adds_its_inverse_covariance_to_the_information() {
        let (mean, c) = ([3.0, 4.0], [[0.25, 0.15], [0.15, 0.36]]);
        let prior = Prior::correlated(&[0, 1], &mean, &two_by_two(c.concat().try_into().unwrap()));
        let result = fit_mixture(&[prior.unwrap()]);
        assert!(result.converged);
        let precision = inverse(c);
        let information = |theta: [f64; 2], observed: bool| {
            let mut information = precision;
            for (a, y) in MIXTURE.iter().zip(MIXTURE_COUNTS) {
                let mu = a[0] * theta[0] + a[1] * theta[1];
                let weight = if observed { y / mu.powi(2) } else { 1.0 / mu };
                for k in 0..2 {
                    for l in 0..2 {
                        information[k][l] += weight * a[k] * a[l];
                    }
                }
            }
            information
        };
        let mut mode = mean;
        for _ in 0..50 {
            let gradient: Vec<f64> = (0..2)
                .map(|k| {
                    (0..2)
                        .map(|l| precision[k][l] * (mode[l] - mean[l]))
                        .sum::<f64>()
                        + MIXTURE
                            .iter()
                            .zip(MIXTURE_COUNTS)
                            .map(|(a, y)| a[k] * (1.0 - y / (a[0] * mode[0] + a[1] * mode[1])))
                            .sum::<f64>()
                })
                .collect();
            let step = inverse(information(mode, true));
            for k in 0..2 {
                mode[k] -= (0..2).map(|l| step[k][l] * gradient[l]).sum::<f64>();
            }
        }
        let fitted = [result.params[0], result.params[1]];
        let expected = information(fitted, false);
        let shift: Vec<f64> = (0..2).map(|k| fitted[k] - mode[k]).collect();
        let distance = (0..2)
            .flat_map(|k| (0..2).map(move |l| (k, l)))
            .map(|(k, l)| shift[k] * expected[k][l] * shift[l])
            .sum::<f64>()
            .sqrt();
        assert!(distance <= 2e-3, "{fitted:?} vs {mode:?}");
        let covariance = inverse(expected);
        let reported = result.covariance.expect("covariance");
        for k in 0..2 {
            for l in 0..2 {
                let scale = (covariance[k][k] * covariance[l][l]).sqrt();
                assert!(
                    (reported.get(k, l) - covariance[k][l]).abs() <= 1e-12 * scale,
                    "{k}{l}: {} vs {}",
                    reported.get(k, l),
                    covariance[k][l]
                );
            }
        }
        let unbounded = result.unbounded.expect("unbounded");
        assert_eq!(unbounded.covariance.data, reported.data);
        for k in 0..2 {
            assert!((unbounded.mean[k] - fitted[k]).abs() <= 2e-3 * covariance[k][k].sqrt());
        }
    }

    const FAR_COUNTS: [f64; 5] = [9743.0, 20046.0, 30110.0, 39769.0, 49761.0];

    fn fit_far(params: Vec<FitParameter>, priors: &[Prior]) -> PoissonResult {
        let mut params = ParameterSet::new(params);
        poisson_fit(
            &Mixture,
            &FAR_COUNTS,
            priors,
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap()
    }

    #[test]
    fn the_unbounded_gaussian_is_one_newton_step_past_the_bounds() {
        let (x, y) = (
            FitParameter::non_negative("x", 9000.0),
            FitParameter::non_negative("y", 1.0),
        );
        let capped = FitParameter {
            upper: 9000.0,
            ..x.clone()
        };
        let cases = [
            (vec![x.clone(), y.clone()], None, 1),
            (vec![capped, FitParameter::unbounded("y", 1.0)], None, 0),
            (vec![x, y], Some((10.0, 20.0)), 1),
        ];
        for (params, prior, held) in cases {
            let priors: Vec<Prior> = prior
                .iter()
                .map(|&(mean, sd)| Prior::measured(1, mean, sd))
                .collect();
            let result = fit_far(params, &priors);
            assert!(result.converged && result.on_bound[held] && !result.on_bound[1 - held]);
            let theta = [result.params[0], result.params[1]];
            let (mut information, mut gradient) = ([[0.0; 2]; 2], [0.0; 2]);
            for (a, n) in MIXTURE.iter().zip(FAR_COUNTS) {
                let mu = a[0] * theta[0] + a[1] * theta[1];
                for k in 0..2 {
                    gradient[k] += a[k] * (1.0 - n / mu);
                    for l in 0..2 {
                        information[k][l] += a[k] * a[l] / mu;
                    }
                }
            }
            if let Some((mean, sd)) = prior {
                information[1][1] += sd.powi(-2);
                gradient[1] += (theta[1] - mean) / sd.powi(2);
            }
            let covariance = inverse(information);
            let unbounded = result.unbounded.expect("unbounded");
            for k in 0..2 {
                let mean = theta[k] - (0..2).map(|l| covariance[k][l] * gradient[l]).sum::<f64>();
                let sd = covariance[k][k].sqrt();
                assert!(
                    (unbounded.mean[k] - mean).abs() <= 1e-12 * sd,
                    "{held} {k}: {} vs {mean}",
                    unbounded.mean[k]
                );
                for (l, expected) in covariance[k].iter().enumerate() {
                    let scale = sd * covariance[l][l].sqrt();
                    assert!((unbounded.covariance.get(k, l) - expected).abs() <= 1e-12 * scale);
                }
            }
        }
    }

    struct Unused;

    impl FitModel for Unused {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(vec![params[0]])
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(1, free_param_indices.len());
            for (col, &index) in free_param_indices.iter().enumerate() {
                *jacobian.get_mut(0, col) = f64::from(u8::from(index == 0));
            }
            Some(jacobian)
        }
    }

    #[test]
    fn a_parameter_on_a_bound_the_counts_do_not_determine_leaves_the_others_their_gaussian() {
        let mut params = ParameterSet::new(vec![
            FitParameter::unbounded("x", 90.0),
            FitParameter::non_negative("y", 0.0),
        ]);
        let result = poisson_fit(
            &Unused,
            &[100.0],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(result.converged && result.on_bound == [false, true]);
        let unbounded = result.unbounded.expect("unbounded");
        let variance = unbounded.covariance.get(0, 0);
        assert!((variance / result.params[0] - 1.0).abs() <= 1e-12);
        assert!((unbounded.mean[0] - 100.0).abs() <= 2e-3 * variance.sqrt());
        assert!(unbounded.covariance.get(1, 1).is_nan() && unbounded.mean[1].is_nan());
    }

    #[test]
    fn the_unbounded_gaussian_is_the_fit_without_its_bound_to_second_order() {
        let x = FitParameter::non_negative("x", 9000.0);
        let bounded = fit_far(vec![x.clone(), FitParameter::non_negative("y", 1.0)], &[]);
        let free = fit_far(vec![x, FitParameter::unbounded("y", 1.0)], &[]);
        let covariance = free.covariance.expect("covariance");
        assert!(bounded.on_bound[1] && !free.on_bound[1]);
        assert!(free.params[1] <= -1.5 * covariance.get(1, 1).sqrt());
        let unbounded = bounded.unbounded.expect("unbounded");
        for k in 0..2 {
            let variance = covariance.get(k, k);
            assert!((unbounded.mean[k] - free.params[k]).abs() <= 0.05 * variance.sqrt());
            assert!((unbounded.covariance.get(k, k) / variance - 1.0).abs() <= 0.05);
        }
    }

    #[test]
    fn a_correlated_prior_with_a_diagonal_covariance_is_two_measured_priors() {
        let (mean, sd) = ([3.0, 4.0], [0.5, 0.6]);
        let diagonal = two_by_two([sd[0] * sd[0], 0.0, 0.0, sd[1] * sd[1]]);
        let bits = |result: PoissonResult| -> Vec<u64> {
            result
                .params
                .iter()
                .chain(&result.covariance.expect("covariance").data)
                .chain([&result.deviance])
                .map(|v| v.to_bits())
                .collect()
        };
        assert_eq!(
            bits(fit_mixture(&[
                Prior::correlated(&[0, 1], &mean, &diagonal).unwrap()
            ])),
            bits(fit_mixture(&[
                Prior::measured(0, mean[0], sd[0]),
                Prior::measured(1, mean[1], sd[1]),
            ]))
        );
    }

    #[test]
    fn predictions_stay_positive_where_nothing_was_counted() {
        let model = Line {
            x: vec![1.0, 2.0, 3.0],
        };
        let mut params = ParameterSet::new(vec![FitParameter::unbounded("a", 0.0)]);
        let result = poisson_fit(
            &model,
            &[0.0; 3],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        let predicted = model.evaluate(&result.params).unwrap();
        assert!(predicted.iter().all(|&mean| mean > 0.0), "{result:?}");
    }

    struct WrongShape;

    impl FitModel for WrongShape {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(vec![params[0]; 3])
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            Some(FlatMatrix::zeros(2, free_param_indices.len()))
        }
    }

    #[test]
    fn a_jacobian_of_the_wrong_shape_is_refused() {
        let mut params = ParameterSet::new(vec![FitParameter::non_negative("a", 1.0)]);
        let result = poisson_fit(
            &WrongShape,
            &[1.0; 3],
            &[],
            &mut params,
            &PoissonConfig::default(),
        );
        assert!(
            matches!(result, Err(FittingError::LengthMismatch { .. })),
            "{result:?}"
        );
    }

    #[test]
    fn inverted_or_nan_bounds_are_refused() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        for (lower, upper) in [
            (1.0, 0.0),
            (f64::NAN, 1.0),
            (0.0, f64::NAN),
            (f64::INFINITY, f64::INFINITY),
            (f64::NEG_INFINITY, f64::NEG_INFINITY),
        ] {
            let mut params = ParameterSet::new(vec![
                FitParameter {
                    name: "a".into(),
                    value: 50.0,
                    lower,
                    upper,
                    fixed: false,
                },
                FitParameter::non_negative("b", 1.0),
            ]);
            let result = poisson_fit(
                &model,
                &observed,
                &[],
                &mut params,
                &PoissonConfig::default(),
            );
            assert!(
                matches!(result, Err(FittingError::InvalidConfig(_))),
                "{lower}, {upper}"
            );
        }
    }

    struct Additive;

    impl FitModel for Additive {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(vec![2.0 * params[0], 1.0 - params[0]])
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(2, free_param_indices.len());
            jacobian.data.copy_from_slice(&[2.0, -1.0]);
            Some(jacobian)
        }
    }

    #[test]
    fn a_bin_predicted_zero_keeps_its_slope_in_the_gradient() {
        let mut params = ParameterSet::new(vec![FitParameter {
            name: "theta".into(),
            value: 0.0,
            lower: 0.0,
            upper: 1.0,
            fixed: false,
        }]);
        let result = poisson_fit(
            &Additive,
            &[0.0, 0.0],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(
            result.converged
                && result.params[0] == 0.0
                && result.on_bound == vec![true]
                && result.unbounded.is_none(),
            "{result:?}"
        );
    }

    struct Line2 {
        offset: [f64; 2],
        slope: [f64; 2],
    }

    impl FitModel for Line2 {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok((0..2)
                .map(|i| self.offset[i] + self.slope[i] * params[0])
                .collect())
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(2, free_param_indices.len());
            jacobian.data.copy_from_slice(&self.slope);
            Some(jacobian)
        }
    }

    fn bounded(value: f64, upper: f64) -> ParameterSet {
        ParameterSet::new(vec![FitParameter {
            name: "theta".into(),
            value,
            lower: 0.0,
            upper,
            fixed: false,
        }])
    }

    #[test]
    fn a_zero_bin_slope_enters_the_convergence_test() {
        let model = Line2 {
            offset: [0.0, 4.0],
            slope: [1.0, -4.0],
        };
        let result = poisson_fit(
            &model,
            &[0.0, 3.0],
            &[],
            &mut bounded(0.0, 1.5),
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(
            result.converged && result.iterations == 0 && result.params[0] == 0.0,
            "{result:?}"
        );
    }

    struct Affine {
        offset: Vec<f64>,
        jacobian: Vec<Vec<f64>>,
    }

    impl FitModel for Affine {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(self
                .offset
                .iter()
                .zip(&self.jacobian)
                .map(|(offset, row)| {
                    offset + row.iter().zip(params).map(|(j, p)| j * p).sum::<f64>()
                })
                .collect())
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(self.offset.len(), free_param_indices.len());
            jacobian.data.copy_from_slice(&self.jacobian.concat());
            Some(jacobian)
        }
    }

    fn two_free() -> ParameterSet {
        ParameterSet::new(vec![
            FitParameter::unbounded("a", 0.0),
            FitParameter::unbounded("b", 0.0),
        ])
    }

    fn at_the_start() -> PoissonConfig {
        PoissonConfig {
            max_iter: 0,
            ..PoissonConfig::default()
        }
    }

    #[test]
    fn a_weak_direction_keeps_its_share_of_the_newton_decrement() {
        let model = Affine {
            offset: vec![1.0; 2],
            jacobian: vec![vec![1.0, 1.0], vec![1.0e-14, -1.0e-14]],
        };
        for observed in [[0.999, 0.998999], [0.999, 0.999001]] {
            let decrement = 0.5
                * observed
                    .iter()
                    .map(|y: &f64| (1.0 - y).powi(2))
                    .sum::<f64>();
            let result =
                poisson_fit(&model, &observed, &[], &mut two_free(), &at_the_start()).unwrap();
            assert_eq!(
                result.converged,
                decrement < NEWTON_DECREMENT_TOL,
                "{decrement:e}: {result:?}"
            );
        }
    }

    #[test]
    fn the_information_inverse_is_poisson_fits_covariance() {
        for weak in [0.0, 1.0e-6] {
            let jacobian = vec![
                vec![1.0, 1.0 + weak, 0.0],
                vec![2.0, 2.0, 1.0],
                vec![1.0, 1.0, 3.0],
                vec![3.0, 3.0 - weak, 1.0],
            ];
            let model = Affine {
                offset: vec![5.0; 4],
                jacobian: jacobian.clone(),
            };
            let mut parameters = ParameterSet::new(
                ["a", "b", "c"]
                    .map(|name| FitParameter::unbounded(name, 0.0))
                    .to_vec(),
            );
            let result =
                poisson_fit(&model, &[5.0; 4], &[], &mut parameters, &at_the_start()).unwrap();
            let covariance = result.covariance.expect("covariance");
            let information = FlatMatrix {
                data: (0..9)
                    .map(|e| {
                        jacobian
                            .iter()
                            .map(|row| row[e / 3] * row[e % 3] / 5.0)
                            .sum()
                    })
                    .collect(),
                nrows: 3,
                ncols: 3,
            };
            let diagonal: Vec<f64> = (0..3).map(|i| information.get(i, i)).collect();
            let inverse = information_inverse(&information, &diagonal, 4).unwrap();
            for i in 0..3 {
                assert_eq!(
                    inverse.resolved[i],
                    !covariance.get(i, i).is_nan(),
                    "{weak}: {i}"
                );
                for j in (0..3).filter(|&j| inverse.resolved[i] && inverse.resolved[j]) {
                    let (got, expected) = (inverse.determined.get(i, j), covariance.get(i, j));
                    assert!(
                        (got - expected).abs() <= 1e-9 * expected.abs().max(1.0),
                        "{weak}: {got} vs {expected}"
                    );
                }
            }
            assert_eq!(inverse.spanned.get(0, 0) > 1e12, weak > 0.0, "{weak}");
        }
    }

    #[test]
    fn an_information_with_nothing_determined_inverts_to_zero() {
        let inverse = information_inverse(&FlatMatrix::zeros(2, 2), &[0.0; 2], 4).unwrap();
        assert!(inverse.spanned.data.iter().all(|&x| x == 0.0));
        assert!(inverse.determined.data.iter().all(|&x| x == 0.0));
        assert_eq!(inverse.resolved, [false; 2]);
        let empty = information_inverse(&FlatMatrix::zeros(0, 0), &[], 4).unwrap();
        assert!(empty.resolved.is_empty());
    }

    #[test]
    fn a_direction_the_steps_leave_out_has_no_error_bar() {
        let weak = 2.0e-12;
        let information = FlatMatrix {
            data: vec![1.0, 1.0 - weak, 1.0 - weak, 1.0],
            nrows: 2,
            ncols: 2,
        };
        for (rows, determined) in [(4, true), (100_000, false)] {
            let inverse = information_inverse(&information, &[1.0; 2], rows).unwrap();
            assert_eq!(inverse.resolved, [determined; 2], "{rows}");
            assert_eq!(
                inverse.undetermined.len(),
                usize::from(!determined),
                "{rows}"
            );
        }
    }

    #[test]
    fn a_zero_bin_slope_enters_every_direction_of_the_decrement() {
        let jacobian = vec![vec![1.0, 2.0], vec![1.0, 3.0], vec![2.0, -1.0]];
        let offset: Vec<f64> = vec![0.0, 4.0, 9.0];
        let root = [offset[1].sqrt(), offset[2].sqrt()];
        let w = [
            [jacobian[1][0] / root[0], jacobian[1][1] / root[0]],
            [jacobian[2][0] / root[1], jacobian[2][1] / root[1]],
        ];
        let information = |a: usize, b: usize| w[0][a] * w[0][b] + w[1][a] * w[1][b];
        let inverse_00 =
            information(1, 1) / (information(0, 0) * information(1, 1) - information(0, 1).powi(2));
        let transpose_det = w[0][0] * w[1][1] - w[1][0] * w[0][1];
        for decrement in [0.99e-6_f64, 1.01e-6] {
            let gradient = [(2.0 * decrement / inverse_00).sqrt(), 0.0];
            let rhs = [gradient[0] - jacobian[0][0], gradient[1] - jacobian[0][1]];
            let residual = [
                (rhs[0] * w[1][1] - w[1][0] * rhs[1]) / transpose_det,
                (w[0][0] * rhs[1] - w[0][1] * rhs[0]) / transpose_det,
            ];
            let observed = [
                0.0,
                offset[1] - residual[0] * root[0],
                offset[2] - residual[1] * root[1],
            ];
            let model = Affine {
                offset: offset.clone(),
                jacobian: jacobian.clone(),
            };
            let result =
                poisson_fit(&model, &observed, &[], &mut two_free(), &at_the_start()).unwrap();
            assert_eq!(
                result.converged,
                decrement < NEWTON_DECREMENT_TOL,
                "{decrement:e}: {result:?}"
            );
        }
    }

    #[test]
    fn a_fit_with_hundreds_of_parameters_has_the_error_bars_of_its_inverse_information() {
        use rand::SeedableRng;
        use rand_chacha::ChaCha12Rng;
        use rand_distr::{Distribution, StandardNormal};
        let (bins, count, level) = (300, 150, 1.0e4);
        let mut rng = ChaCha12Rng::seed_from_u64(806);
        let jacobian: Vec<Vec<f64>> = (0..bins)
            .map(|_| {
                (0..count)
                    .map(|_| StandardNormal.sample(&mut rng))
                    .collect()
            })
            .collect();
        let model = Affine {
            offset: vec![level; bins],
            jacobian: jacobian.clone(),
        };
        let mut params = ParameterSet::new(
            (0..count)
                .map(|j| FitParameter::unbounded(format!("c{j}"), 0.0))
                .collect(),
        );
        let result = poisson_fit(
            &model,
            &vec![level; bins],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(result.converged && result.iterations == 0, "{result:?}");
        let mut cholesky = vec![vec![0.0; count]; count];
        for a in 0..count {
            for b in 0..=a {
                let information: f64 = jacobian.iter().map(|row| row[a] * row[b] / level).sum();
                let rest =
                    information - (0..b).map(|c| cholesky[a][c] * cholesky[b][c]).sum::<f64>();
                cholesky[a][b] = if a == b {
                    rest.sqrt()
                } else {
                    rest / cholesky[b][b]
                };
            }
        }
        for (j, error) in result.uncertainties.unwrap().into_iter().enumerate() {
            let mut solved = vec![0.0; count];
            for a in 0..count {
                let unit = if a == j { 1.0 } else { 0.0 };
                solved[a] = (unit - (0..a).map(|c| cholesky[a][c] * solved[c]).sum::<f64>())
                    / cholesky[a][a];
            }
            let expected = solved.iter().map(|v| v * v).sum::<f64>().sqrt();
            let error = error.unwrap();
            assert!(
                (error / expected - 1.0).abs() < 1e-10,
                "{j}: {error} vs {expected}"
            );
        }
    }

    struct SilentZeroBin;

    impl FitModel for SilentZeroBin {
        fn evaluate(&self, _params: &[f64]) -> Result<Vec<f64>, FittingError> {
            Ok(vec![0.0])
        }

        fn analytical_jacobian(
            &self,
            _params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let mut jacobian = FlatMatrix::zeros(1, free_param_indices.len());
            jacobian.data[0] = f64::NAN;
            Some(jacobian)
        }
    }

    #[test]
    fn a_non_finite_slope_in_a_zero_bin_ends_the_fit_unconverged() {
        let mut params = ParameterSet::new(vec![FitParameter::unbounded("theta", 0.0)]);
        let result = poisson_fit(
            &SilentZeroBin,
            &[0.0],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert!(!result.converged, "{result:?}");
    }

    struct Wall;

    impl FitModel for Wall {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            let x = params[0];
            Ok(vec![(100.0 - x + (x - 80.0).max(0.0).powi(2)).exp()])
        }

        fn analytical_jacobian(
            &self,
            params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let x = params[0];
            let mean = (100.0 - x + (x - 80.0).max(0.0).powi(2)).exp();
            let mut jacobian = FlatMatrix::zeros(1, free_param_indices.len());
            jacobian.data[0] = mean * (2.0 * (x - 80.0).max(0.0) - 1.0);
            Some(jacobian)
        }
    }

    #[test]
    fn damping_recovers_after_many_accepted_steps_to_reach_the_minimum() {
        let mut params = ParameterSet::new(vec![FitParameter::unbounded("x", 0.0)]);
        let result =
            poisson_fit(&Wall, &[0.0], &[], &mut params, &PoissonConfig::default()).unwrap();
        assert!((result.params[0] - 80.5).abs() < 1e-3, "{result:?}");
    }

    #[test]
    fn non_finite_parameter_values_are_refused() {
        let model = decay(1.0);
        let observed = model.evaluate(&[100.0, 1.5]).unwrap();
        for (a, b) in [(f64::NAN, 1.0), (50.0, f64::INFINITY)] {
            let mut params = ParameterSet::new(vec![
                FitParameter::non_negative("a", a),
                FitParameter::fixed("b", b),
            ]);
            let result = poisson_fit(
                &model,
                &observed,
                &[],
                &mut params,
                &PoissonConfig::default(),
            );
            assert!(
                matches!(result, Err(FittingError::InvalidConfig(_))),
                "{a}, {b}"
            );
        }
    }

    struct Shrinking;

    impl FitModel for Shrinking {
        fn evaluate(&self, params: &[f64]) -> Result<Vec<f64>, FittingError> {
            let bins = if params[0] < 0.5 { 3 } else { 2 };
            Ok(vec![1.0 + params[0]; bins])
        }

        fn analytical_jacobian(
            &self,
            params: &[f64],
            free_param_indices: &[usize],
            _y_current: &[f64],
        ) -> Option<FlatMatrix> {
            let bins = if params[0] < 0.5 { 3 } else { 2 };
            let mut jacobian = FlatMatrix::zeros(bins, free_param_indices.len());
            jacobian.data.fill(1.0);
            Some(jacobian)
        }
    }

    #[test]
    fn a_trial_prediction_of_the_wrong_length_is_rejected() {
        let mut params = ParameterSet::new(vec![FitParameter::non_negative("theta", 0.0)]);
        let result = poisson_fit(
            &Shrinking,
            &[3.0; 3],
            &[],
            &mut params,
            &PoissonConfig::default(),
        )
        .unwrap();
        assert_eq!(
            Shrinking.evaluate(&result.params).unwrap().len(),
            3,
            "{result:?}"
        );
    }
}