regit-svi 2.0.0

// Copyright 2026 Regit.io — Nicolas Koenig
// SPDX-License-Identifier: Apache-2.0

//! Quasi-explicit slice calibration (de Marco & Martini / Zeliade 2009).
//!
//! Direct least-squares over all five raw parameters is non-convex and
//! sensitive to the starting point. The quasi-explicit method removes that
//! fragility by exploiting a change of variables that makes the problem
//! linear in three of the five parameters.
//!
//! # Reduction
//!
//! Fix the two nonlinear parameters `m` and `sigma`. Substituting
//! `y = (k - m)/sigma` gives `sqrt((k-m)^2 + sigma^2) = sigma*sqrt(y^2 + 1)`,
//! so with `c = b*sigma` and `d = rho*b*sigma` the model becomes affine in
//! `(a, d, c)`:
//!
//! ```text
//! w(y) = a + d*y + c*sqrt(y^2 + 1)
//! ```
//!
//! # Inner problem — convex, solved numerically
//!
//! For fixed `(m, sigma)`, minimise the weighted residual over `(a, d, c)` on
//! the convex Zeliade domain `D`:
//!
//! ```text
//! 0 <= c <= 4*sigma
//! |d| <= c          and      |d| <= 4*sigma - c
//! 0 <= a <= max_i w_i
//! ```
//!
//! In exact arithmetic `f` is a convex quadratic, so its minimum over the polytope `D` is either
//! the unconstrained stationary point (from the 3x3 normal equations) or, if
//! that is infeasible, lies on a face — found by enumerating every face,
//! edge, and vertex of `D` and taking the feasible minimiser. With three
//! variables this enumeration is small. The floating-point normal equations
//! and candidate comparisons remain approximate numerical work.
//!
//! # Outer problem — two-dimensional
//!
//! Let `f*(m, sigma)` be the optimal inner residual. The 2-D, mildly
//! non-convex problem `min f*(m, sigma)` is solved with the Nelder-Mead
//! simplex, multi-started across a small grid of `(m, sigma)` seeds.
//!
//! # Recovery
//!
//! From the optimal `(a, c, d, m, sigma)`: `b = c/sigma`, `rho = d/c`
//! (`b = 0`, `rho = 0` when `c = 0`).
//!
//! # References
//!
//! - De Marco, S. & Martini, C., "Quasi-explicit calibration of Gatheral's
//!   SVI model", Zeliade Systems White Paper ZWP-0005 (2009).

use core::cell::Cell;

use crate::calibration::CalibrationResult;
use crate::calibration::config::{
    ConstraintMode, InitializationPolicy, SliceCalibrationConfig, validate_effective_quotes,
};
use crate::calibration::report::{ParameterizationEvidence, TerminationReason};
use crate::calibration::slice::residual_diagnostics;
use crate::error::CalibrationError;
use crate::market::quote::Quote;
use crate::numerics::{nelder_mead, solve_spd_3};
use crate::smile::raw::RawSvi;

/// Minimum number of quotes the five-parameter raw SVI model can be fit to.
const MIN_QUOTES: usize = 5;

/// Calibrates a raw SVI slice to market quotes by the quasi-explicit method.
///
/// Runs the inner convex QP in `(a, d, c)` for each `(m, sigma)` proposed by a
/// multi-started Nelder-Mead outer search, recovers `b` and `rho`, and
/// returns the fitted slice with residual, termination, and arbitrage evidence.
///
/// # Errors
///
/// - [`CalibrationError::EmptyQuotes`] if `quotes` is empty.
/// - [`CalibrationError::InsufficientEffectiveQuotes`] if fewer than five
///   positive-weight, sufficiently distinct strikes are supplied.
/// - [`CalibrationError::AllWeightsZero`] if every fitting weight is zero.
/// - [`CalibrationError::Param`] if aggregate weight or recovered Raw state is
///   non-finite or invalid.
/// - [`CalibrationError::DidNotConverge`] if no eligible candidate is found.
/// - [`CalibrationError::Infeasible`] if constrained mode cannot establish the
///   required Raw butterfly conditions.
///
/// # Examples
///
/// ```
/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
/// use regit_svi::market::quote::Quote;
/// use regit_svi::calibration::quasi_explicit::calibrate;
///
/// let quotes = [
///     Quote::new(-0.20, 0.0512, 1.0)?,
///     Quote::new(-0.10, 0.0432, 1.0)?,
///     Quote::new( 0.00, 0.0400, 1.0)?,
///     Quote::new( 0.10, 0.0420, 1.0)?,
///     Quote::new( 0.20, 0.0480, 1.0)?,
/// ];
/// let fit = calibrate(&quotes)?;
/// assert!(fit.rmse() < 1e-2);
/// # Ok(())
/// # }
/// ```
pub fn calibrate(quotes: &[Quote]) -> Result<CalibrationResult, CalibrationError> {
    calibrate_with_config(quotes, SliceCalibrationConfig::default())
}

/// Calibrates by the quasi-explicit method with explicit numerical controls.
///
/// # Errors
///
/// Returns [`CalibrationError::EmptyQuotes`] or
/// [`CalibrationError::AllWeightsZero`] when no usable observations exist,
/// [`CalibrationError::InsufficientEffectiveQuotes`] when fewer than five
/// positive-weight distinct strikes are available,
/// [`CalibrationError::Param`] for non-finite aggregate weight or invalid
/// recovered Raw state, [`CalibrationError::DidNotConverge`] when no eligible
/// candidate is found, and [`CalibrationError::Infeasible`] when constrained
/// mode cannot establish the required Raw butterfly conditions.
#[allow(clippy::too_many_lines)] // Keeps fallback, multi-start selection, and report evidence in execution order.
pub fn calibrate_with_config(
    quotes: &[Quote],
    config: SliceCalibrationConfig,
) -> Result<CalibrationResult, CalibrationError> {
    validate_effective_quotes(quotes, MIN_QUOTES, config.distinct_tolerance())?;

    let k_min = quotes.iter().map(|q| q.k).fold(f64::INFINITY, f64::min);
    let k_max = quotes.iter().map(|q| q.k).fold(f64::NEG_INFINITY, f64::max);
    let k_span = (k_max - k_min).max(1e-3);
    let w_max = quotes
        .iter()
        .map(|q| q.w)
        .fold(0.0_f64, f64::max)
        .max(1e-12);

    // The outer objective: optimal inner residual at (m, sigma).
    // Parametrise sigma = exp(sigma_hat) to keep it strictly positive.
    let evaluations = Cell::new(0_usize);
    let outer = |p: &[f64]| -> f64 {
        evaluations.set(evaluations.get().saturating_add(1));
        let m = p[0];
        let sigma = p[1].exp();
        if !m.is_finite() || !sigma.is_finite() || sigma <= 0.0 {
            return f64::INFINITY;
        }
        inner_solve(quotes, m, sigma, w_max).0
    };

    // Multi-start grid of (m, sigma) seeds covering the quoted range.
    let m_seeds = [
        k_min,
        0.5 * (k_min + k_max),
        k_max,
        k_min - 0.25 * k_span,
        k_max + 0.25 * k_span,
    ];
    let sigma_seeds = [0.1 * k_span, 0.3 * k_span, k_span, 2.0 * k_span];

    let mut best_obj = f64::INFINITY;
    let mut best_slice = None;
    let mut best_start = None;
    let mut best_condition = None;
    let mut selected_termination = crate::numerics::OptimizerTermination::IterationLimit;
    let mut total_iterations = 0_usize;

    // Best-effort mode has an explicit valid fallback rather than coercing a
    // boundary QP solution into the open Raw domain. Its report retains the
    // non-converged termination unless a converged start replaces it.
    if config.constraint_mode() == ConstraintMode::BestEffort {
        let total_weight: f64 = quotes.iter().map(|quote| quote.weight()).sum();
        let flat_level = quotes
            .iter()
            .map(|quote| quote.weight() * quote.total_variance())
            .sum::<f64>()
            / total_weight;
        if let Ok(fallback) = RawSvi::new(flat_level, 0.0, 0.0, 0.5 * (k_min + k_max), 0.3 * k_span)
        {
            let fallback_objective =
                residual_diagnostics(&fallback, quotes, config.distinct_tolerance()).objective();
            if fallback_objective.is_finite() {
                best_obj = fallback_objective;
                best_slice = Some(fallback);
            }
        }
    }

    let m_count = if config.initialization() == InitializationPolicy::DeterministicSingleStart {
        1
    } else {
        m_seeds.len()
    };
    let sigma_count = if config.initialization() == InitializationPolicy::DeterministicSingleStart {
        1
    } else {
        sigma_seeds.len()
    };
    for (m_index, &m0) in m_seeds[..m_count].iter().enumerate() {
        for (sigma_index, &s0) in sigma_seeds[..sigma_count].iter().enumerate() {
            let start = [m0, s0.max(1e-6).ln()];
            let res = nelder_mead(outer, &start, config.tolerance(), config.outer_iterations());
            total_iterations = total_iterations.saturating_add(res.iterations);
            let eligible = res.converged || config.constraint_mode() == ConstraintMode::BestEffort;
            if eligible {
                let [candidate_m, candidate_log_sigma] = res.x.as_slice() else {
                    continue;
                };
                let candidate_sigma = candidate_log_sigma.exp();
                let (_, a, d, c, condition_estimate) =
                    inner_solve(quotes, *candidate_m, candidate_sigma, w_max);
                let b = c / candidate_sigma;
                let rho = if c.abs() > 1e-300 { d / c } else { 0.0 };
                // The closed polytope permits |d|=c, but Raw SVI requires the
                // open domain |rho|<1. Do not clamp a boundary optimizer into
                // a different model; discard non-representable candidates.
                let Ok(candidate) = RawSvi::new(a, b, rho, *candidate_m, candidate_sigma) else {
                    continue;
                };
                let candidate_objective =
                    residual_diagnostics(&candidate, quotes, config.distinct_tolerance())
                        .objective();
                if candidate_objective.is_finite() && candidate_objective < best_obj {
                    best_obj = candidate_objective;
                    best_slice = Some(candidate);
                    best_start = Some(m_index * sigma_count + sigma_index);
                    best_condition = Some(condition_estimate);
                    selected_termination = res.termination;
                }
            }
        }
    }
    let slice = best_slice.ok_or(CalibrationError::DidNotConverge {
        iterations: total_iterations,
        residual: best_obj,
    })?;
    let residuals = residual_diagnostics(&slice, quotes, config.distinct_tolerance());
    let termination = TerminationReason::from(selected_termination);
    CalibrationResult::from_parts(
        slice,
        residuals,
        total_iterations,
        evaluations.get(),
        termination,
        config,
        "quasi-explicit / Nelder–Mead",
        m_count * sigma_count,
        best_start,
        best_condition,
        ParameterizationEvidence::new(
            "sigma=exp; affine (a,d,c) polytope recovery",
            "no projection; non-representable open-domain candidates discarded",
        ),
    )
}

/// Solves the inner convex QP for fixed `(m, sigma)`.
///
/// Returns `(residual, a, d, c, condition)`: the optimal weighted sum of
/// squared residuals, the affine coefficients of
/// `w(y) = a + d*y + c*sqrt(y^2+1)`, and the normalized design condition
/// estimate, minimised over the Zeliade domain `D` (MATH.md §10).
///
/// The face/edge/vertex enumeration of the 3-D polytope `D` is exhaustive by
/// nature, so the function is necessarily long; it is one cohesive solver.
#[allow(clippy::too_many_lines)]
fn inner_solve(quotes: &[Quote], m: f64, sigma: f64, w_max: f64) -> (f64, f64, f64, f64, f64) {
    // Design data: phi_i = (1, y_i, z_i) with y_i = (k_i-m)/sigma,
    // z_i = sqrt(y_i^2 + 1).
    let max_weight = quotes
        .iter()
        .map(|quote| quote.weight())
        .fold(0.0_f64, f64::max);
    if max_weight <= 0.0 || !max_weight.is_finite() {
        return (f64::INFINITY, 0.0, 0.0, 0.0, f64::INFINITY);
    }
    let mut rows: Vec<([f64; 3], f64, f64)> = Vec::with_capacity(quotes.len());
    for q in quotes {
        if q.weight <= 0.0 {
            continue;
        }
        let y = (q.k - m) / sigma;
        let z = y.hypot(1.0);
        // A common rescaling of all weights must not change the inner fit.
        // Dividing by the largest weight also keeps the Gram accumulation in
        // a representable range without changing the QP minimizer.
        rows.push(([1.0, y, z], q.w, q.weight / max_weight));
    }
    if rows.is_empty() {
        return (f64::INFINITY, 0.0, 0.0, 0.0, f64::INFINITY);
    }

    // Weighted normal-equations matrices: A (symmetric 3x3) and rhs (3).
    // A_jl = sum_i weight_i * phi_ij * phi_il
    // rhs_j = sum_i weight_i * phi_ij * w_i
    let mut a_mat = [0.0_f64; 6]; // [a00, a01, a02, a11, a12, a22]
    let mut rhs = [0.0_f64; 3];
    for (phi, w, weight) in &rows {
        let ww = *weight;
        a_mat[0] += ww * phi[0] * phi[0];
        a_mat[1] += ww * phi[0] * phi[1];
        a_mat[2] += ww * phi[0] * phi[2];
        a_mat[3] += ww * phi[1] * phi[1];
        a_mat[4] += ww * phi[1] * phi[2];
        a_mat[5] += ww * phi[2] * phi[2];
        rhs[0] += ww * phi[0] * w;
        rhs[1] += ww * phi[1] * w;
        rhs[2] += ww * phi[2] * w;
    }

    // Residual of a candidate (a, d, c) — used to compare feasible vertices.
    let residual_of = |a: f64, d: f64, c: f64| -> f64 {
        rows.iter()
            .map(|(phi, w, weight)| {
                let model = d.mul_add(phi[1], c.mul_add(phi[2], a));
                let r = model - w;
                weight * r * r
            })
            .sum()
    };

    // Zeliade domain D for (a, d, c):
    //   0 <= a <= w_max
    //   0 <= c <= 4*sigma
    //   |d| <= c   and   |d| <= 4*sigma - c
    let c_hi = 4.0 * sigma;
    let feasible = |a: f64, d: f64, c: f64| -> bool {
        let eps = 1e-9 * (1.0 + c_hi + w_max);
        a >= -eps
            && a <= w_max + eps
            && c >= -eps
            && c <= c_hi + eps
            && d.abs() <= c + eps
            && d.abs() <= c_hi - c + eps
    };

    let mut best_resid = f64::INFINITY;
    let mut best = (0.0_f64, 0.0_f64, 0.0_f64);
    let mut consider = |a: f64, d: f64, c: f64| {
        if a.is_finite() && d.is_finite() && c.is_finite() && feasible(a, d, c) {
            let r = residual_of(a, d, c);
            if r < best_resid {
                best_resid = r;
                best = (a, d, c);
            }
        }
    };

    // (1) Unconstrained stationary point. Normalize every weighted design
    // column to unit norm before solving, then unscale the coefficients.
    let (unconstrained, condition_estimate) = solve_normalized_spd_3(&a_mat, &rhs);
    if let Some(x) = unconstrained {
        consider(x[0], x[1], x[2]);
    }

    // (2) Faces: fix one variable / one inequality to its boundary and solve
    //     the reduced 2x2 least-squares problem in the other two.
    //
    // The polytope D has these face-defining equalities:
    //   a = 0, a = w_max,
    //   c = 0, c = c_hi,
    //   d = c, d = -c, d = c_hi - c, d = -(c_hi - c).
    //
    // For each, solve the unconstrained reduced problem; clamp/clip is not
    // needed because all edges and vertices are enumerated separately below.

    // Faces a = const.
    for &a_fix in &[0.0, w_max] {
        // Minimise over (d, c): normal equations for the 2-vector (d, c).
        // A2 = [[a11, a12],[a12, a22]], rhs2 = [rhs1 - a01*a_fix, rhs2 - a02*a_fix].
        if let Some((d, c)) = solve_2x2(
            a_mat[3],
            a_mat[4],
            a_mat[5],
            rhs[1] - a_mat[1] * a_fix,
            rhs[2] - a_mat[2] * a_fix,
        ) {
            consider(a_fix, d, c);
        }
    }
    // Faces c = const.
    for &c_fix in &[0.0, c_hi] {
        // Minimise over (a, d).
        if let Some((a, d)) = solve_2x2(
            a_mat[0],
            a_mat[1],
            a_mat[3],
            rhs[0] - a_mat[2] * c_fix,
            rhs[1] - a_mat[4] * c_fix,
        ) {
            consider(a, d, c_fix);
        }
    }
    // Faces d = s*c (s = +/-1): substitute d = s*c, minimise over (a, c).
    for &s in &[1.0_f64, -1.0] {
        // Model column for c becomes (s*phi_d + phi_c); rebuild 2x2 in (a, c).
        // A00 = a00, A01 = s*a01 + a02, A11 = a22 + 2s*a12 + a11.
        let a00 = a_mat[0];
        let a01 = s * a_mat[1] + a_mat[2];
        let a11 = a_mat[5] + 2.0 * s * a_mat[4] + a_mat[3];
        let r0 = rhs[0];
        let r1 = s * rhs[1] + rhs[2];
        if let Some((a, c)) = solve_2x2(a00, a01, a11, r0, r1) {
            consider(a, s * c, c);
        }
    }
    // Faces d = s*(c_hi - c): substitute d = s*c_hi - s*c, minimise over (a, c).
    for &s in &[1.0_f64, -1.0] {
        // d = s*c_hi - s*c. Model: a*1 + (s*c_hi - s*c)*phi_d + c*phi_c.
        // = a + s*c_hi*phi_d + c*(phi_c - s*phi_d).
        // 2x2 in (a, c): column_a = 1, column_c = phi_c - s*phi_d.
        // A00 = a00, A01 = a02 - s*a01, A11 = a22 - 2s*a12 + a11.
        // rhs: r0 = rhs0 - s*c_hi*rhs1 ; r1 = (rhs2 - s*rhs1) - s*c_hi*(a12 - s*a11)
        let a00 = a_mat[0];
        let a01 = a_mat[2] - s * a_mat[1];
        let a11 = a_mat[5] - 2.0 * s * a_mat[4] + a_mat[3];
        let r0 = rhs[0] - s * c_hi * a_mat[1];
        let r1 = (rhs[2] - s * rhs[1]) - s * c_hi * (a_mat[4] - s * a_mat[3]);
        if let Some((a, c)) = solve_2x2(a00, a01, a11, r0, r1) {
            consider(a, s * c_hi - s * c, c);
        }
    }

    // (3) Edges: pairs of equalities. Solve the resulting 1-D least squares.
    //     The relevant edges are intersections of a = const with the (d, c)
    //     boundary lines, and of c = const with the d boundary lines.
    for &a_fix in &[0.0, w_max] {
        for &c_fix in &[0.0, c_hi] {
            // Minimise over d alone.
            // residual sum d^2*a11 - 2*d*(rhs1 - a01*a_fix - a12*c_fix) + const.
            let denom = a_mat[3];
            if denom > 0.0 {
                let d = (rhs[1] - a_mat[1] * a_fix - a_mat[4] * c_fix) / denom;
                consider(a_fix, d, c_fix);
            }
        }
        for &s in &[1.0_f64, -1.0] {
            // a fixed, d = s*c: minimise over c.
            let a11 = a_mat[5] + 2.0 * s * a_mat[4] + a_mat[3];
            if a11 > 0.0 {
                let c = (s * (rhs[1] - a_mat[1] * a_fix) + (rhs[2] - a_mat[2] * a_fix)) / a11;
                consider(a_fix, s * c, c);
            }
            // a fixed, d = s*(c_hi - c): minimise over c.
            let a11b = a_mat[5] - 2.0 * s * a_mat[4] + a_mat[3];
            if a11b > 0.0 {
                let r = (rhs[2] - a_mat[2] * a_fix)
                    - s * (rhs[1] - a_mat[1] * a_fix)
                    - s * c_hi * (a_mat[4] - s * a_mat[3]);
                let c = r / a11b;
                consider(a_fix, s * c_hi - s * c, c);
            }
        }
    }

    // (4) Vertices: all corners of D. With the box on a and c, and the
    //     wing constraints |d| <= c, |d| <= c_hi - c, the vertices in (d, c)
    //     are (d, c) in {(0, 0), (0, c_hi), (c, c) with c = c_hi/2 -> d = +-c_hi/2}.
    let half = c_hi / 2.0;
    let dc_vertices = [(0.0, 0.0), (0.0, c_hi), (half, half), (-half, half)];
    for &a_fix in &[0.0, w_max] {
        for &(d, c) in &dc_vertices {
            consider(a_fix, d, c);
        }
    }
    // Also: a free at each (d, c) vertex.
    for &(d, c) in &dc_vertices {
        // minimise over a alone.
        if a_mat[0] > 0.0 {
            let a = (rhs[0] - a_mat[1] * d - a_mat[2] * c) / a_mat[0];
            consider(a, d, c);
        }
    }

    (best_resid, best.0, best.1, best.2, condition_estimate)
}

/// Solves the full normal equations after unit-norm column scaling and returns
/// the one-norm condition estimate of the normalized Gram matrix.
fn solve_normalized_spd_3(matrix: &[f64; 6], rhs: &[f64; 3]) -> (Option<[f64; 3]>, f64) {
    let scales = [matrix[0].sqrt(), matrix[3].sqrt(), matrix[5].sqrt()];
    if scales
        .iter()
        .any(|scale| !scale.is_finite() || *scale <= 0.0)
    {
        return (None, f64::INFINITY);
    }
    let a = matrix[1] / scales[0] / scales[1];
    let b = matrix[2] / scales[0] / scales[2];
    let c = matrix[4] / scales[1] / scales[2];
    let normalized = [1.0, a, b, 1.0, c, 1.0];
    let determinant = 1.0 + 2.0 * a * b * c - a * a - b * b - c * c;
    let matrix_norm = (1.0 + a.abs() + b.abs())
        .max(1.0 + a.abs() + c.abs())
        .max(1.0 + b.abs() + c.abs());
    let rank_threshold = 256.0 * f64::EPSILON * matrix_norm.powi(3);
    if !determinant.is_finite() || determinant <= rank_threshold {
        return (None, f64::INFINITY);
    }
    let inverse_norm = ((1.0 - c * c).abs() + (b * c - a).abs() + (a * c - b).abs())
        .max((b * c - a).abs() + (1.0 - b * b).abs() + (a * b - c).abs())
        .max((a * c - b).abs() + (a * b - c).abs() + (1.0 - a * a).abs())
        / determinant;
    let condition_estimate = matrix_norm * inverse_norm;
    let normalized_rhs = [rhs[0] / scales[0], rhs[1] / scales[1], rhs[2] / scales[2]];
    let solution = solve_spd_3(&normalized, &normalized_rhs).map(|scaled| {
        [
            scaled[0] / scales[0],
            scaled[1] / scales[1],
            scaled[2] / scales[2],
        ]
    });
    (solution, condition_estimate)
}

/// Solves the symmetric `2x2` system `[[a00, a01],[a01, a11]] x = [r0, r1]`.
///
/// Returns `None` if the system is singular or not positive definite.
fn solve_2x2(a00: f64, a01: f64, a11: f64, r0: f64, r1: f64) -> Option<(f64, f64)> {
    if !a00.is_finite() || !a11.is_finite() || a00 <= 0.0 || a11 <= 0.0 {
        return None;
    }
    let scale0 = a00.sqrt();
    let scale1 = a11.sqrt();
    let correlation = a01 / scale0 / scale1;
    let determinant = 1.0 - correlation * correlation;
    if !determinant.is_finite() || determinant <= 128.0 * f64::EPSILON {
        return None;
    }
    let rhs0 = r0 / scale0;
    let rhs1 = r1 / scale1;
    let x0 = (rhs0 - correlation * rhs1) / determinant / scale0;
    let x1 = (rhs1 - correlation * rhs0) / determinant / scale1;
    if x0.is_finite() && x1.is_finite() {
        Some((x0, x1))
    } else {
        None
    }
}

#[cfg(test)]
#[allow(clippy::expect_used)] // Validated fixtures use contextual expectations.
mod tests {
    use super::*;

    /// Generates a synthetic slice of quotes from known raw SVI parameters.
    fn synthetic(svi: &RawSvi, ks: &[f64]) -> Vec<Quote> {
        ks.iter()
            .map(|&k| {
                Quote::new(k, svi.total_variance(k), 1.0)
                    .expect("valid test or documentation fixture")
            })
            .collect()
    }

    #[test]
    fn rejects_empty() {
        assert!(matches!(calibrate(&[]), Err(CalibrationError::EmptyQuotes)));
    }

    #[test]
    fn rejects_too_few_quotes() {
        let q = Quote::new(0.0, 0.04, 1.0).expect("valid test or documentation fixture");
        assert!(matches!(
            calibrate(&[q, q, q]),
            Err(CalibrationError::InsufficientEffectiveQuotes { .. })
        ));
    }

    #[test]
    fn rejects_all_zero_weights() {
        let quotes: Vec<Quote> = (0..6)
            .map(|i| {
                Quote::new(f64::from(i) * 0.1 - 0.3, 0.04, 0.0)
                    .expect("valid test or documentation fixture")
            })
            .collect();
        assert!(matches!(
            calibrate(&quotes),
            Err(CalibrationError::AllWeightsZero)
        ));
    }

    #[test]
    fn constrained_non_convergence_is_not_success() {
        let truth =
            RawSvi::new(0.04, 0.4, -0.3, 0.05, 0.15).expect("valid test or documentation fixture");
        let quotes = synthetic(&truth, &[-0.4, -0.2, 0.0, 0.2, 0.4]);
        let constrained = SliceCalibrationConfig::new(
            1,
            1,
            1e-15,
            1e-10,
            1e-12,
            ConstraintMode::Constrained,
            InitializationPolicy::DeterministicSingleStart,
        )
        .expect("valid test or documentation fixture");
        assert!(matches!(
            calibrate_with_config(&quotes, constrained),
            Err(CalibrationError::DidNotConverge { .. })
        ));

        let best_effort = constrained.with_constraint_mode(ConstraintMode::BestEffort);
        let result = calibrate_with_config(&quotes, best_effort)
            .expect("valid test or documentation fixture");
        assert_eq!(
            result.report().termination(),
            TerminationReason::MaximumIterations
        );
        let recomputed: f64 = quotes
            .iter()
            .map(|quote| {
                let residual =
                    result.slice().total_variance(quote.log_moneyness()) - quote.total_variance();
                quote.weight() * residual * residual
            })
            .sum();
        assert!(result.report().residuals().objective().is_finite());
        assert!(
            (result.report().residuals().objective() - recomputed).abs()
                <= f64::EPSILON * (1.0 + recomputed.abs())
        );
    }

    #[test]
    fn recovers_synthetic_parameters() {
        let truth =
            RawSvi::new(0.04, 0.4, -0.3, 0.05, 0.15).expect("valid test or documentation fixture");
        let ks = [-0.4, -0.25, -0.1, 0.0, 0.1, 0.25, 0.4];
        let quotes = synthetic(&truth, &ks);
        let fit = calibrate(&quotes).expect("valid test or documentation fixture");

        assert!(fit.report().selected_start().is_some());
        assert!(matches!(
            fit.report().condition_estimate(),
            Some(condition) if condition.is_finite() && condition >= 1.0
        ));
        assert!(
            fit.report()
                .parameterization()
                .transform()
                .contains("sigma=exp")
        );
        assert!(
            fit.report()
                .parameterization()
                .repair()
                .contains("discarded")
        );

        // RMSE should be near machine zero for noise-free synthetic data.
        assert!(fit.rmse() < 1e-5, "rmse = {}", fit.rmse());
        // The fitted slice must reproduce total variance everywhere.
        for &k in &[-0.6, -0.2, 0.0, 0.2, 0.6] {
            let err = (fit.slice().total_variance(k) - truth.total_variance(k)).abs();
            assert!(err < 1e-4, "k = {k}, err = {err}");
        }
    }

    #[test]
    fn recovers_symmetric_smile() {
        let truth =
            RawSvi::new(0.03, 0.3, 0.0, 0.0, 0.2).expect("valid test or documentation fixture");
        let ks = [-0.5, -0.3, -0.1, 0.0, 0.1, 0.3, 0.5];
        let quotes = synthetic(&truth, &ks);
        let fit = calibrate(&quotes).expect("valid test or documentation fixture");
        assert!(fit.rmse() < 1e-5, "rmse = {}", fit.rmse());
        assert!(fit.slice().rho.abs() < 1e-2, "rho = {}", fit.slice().rho);
    }

    #[test]
    fn recovers_positive_skew() {
        let truth =
            RawSvi::new(0.05, 0.35, 0.4, -0.1, 0.18).expect("valid test or documentation fixture");
        let ks = [-0.4, -0.2, -0.05, 0.05, 0.2, 0.4, 0.6];
        let quotes = synthetic(&truth, &ks);
        let fit = calibrate(&quotes).expect("valid test or documentation fixture");
        assert!(fit.rmse() < 1e-4, "rmse = {}", fit.rmse());
        for &k in &[-0.3, 0.0, 0.3] {
            let err = (fit.slice().total_variance(k) - truth.total_variance(k)).abs();
            assert!(err < 1e-3, "k = {k}, err = {err}");
        }
    }

    #[test]
    fn graceful_degradation_with_noise() {
        let truth =
            RawSvi::new(0.04, 0.4, -0.3, 0.05, 0.15).expect("valid test or documentation fixture");
        let ks = [-0.4, -0.25, -0.1, 0.0, 0.1, 0.25, 0.4];
        // Deterministic pseudo-noise from a small LCG.
        let mut state = 12_345_u64;
        let quotes: Vec<Quote> = ks
            .iter()
            .map(|&k| {
                state = state
                    .wrapping_mul(6_364_136_223_846_793_005)
                    .wrapping_add(1);
                let noise =
                    (f64::from((state >> 40) as u32) / f64::from(u32::MAX) - 0.5) * 2.0 * 5e-4;
                Quote::new(k, truth.total_variance(k) + noise, 1.0)
                    .expect("valid test or documentation fixture")
            })
            .collect();
        let fit = calibrate(&quotes).expect("valid test or documentation fixture");
        // With small noise the fit should still be close.
        assert!(fit.rmse() < 1e-2, "rmse = {}", fit.rmse());
        let atm_err = (fit.slice().total_variance(0.0) - truth.total_variance(0.0)).abs();
        assert!(atm_err < 5e-3, "atm err = {atm_err}");
    }

    #[test]
    fn solve_2x2_identity() {
        let (x, y) =
            solve_2x2(1.0, 0.0, 1.0, 3.0, 7.0).expect("valid test or documentation fixture");
        assert!((x - 3.0).abs() < 1e-15);
        assert!((y - 7.0).abs() < 1e-15);
    }

    #[test]
    fn solve_2x2_rejects_singular() {
        assert!(solve_2x2(1.0, 1.0, 1.0, 1.0, 1.0).is_none());
    }

    #[test]
    fn normalized_three_by_three_solve_is_invariant_to_column_scale() {
        let scales = [1e-100_f64, 1.0, 1e100];
        let matrix = [
            scales[0] * scales[0],
            0.1 * scales[0] * scales[1],
            0.2 * scales[0] * scales[2],
            scales[1] * scales[1],
            0.3 * scales[1] * scales[2],
            scales[2] * scales[2],
        ];
        let rhs = [1.8 * scales[0], 3.0 * scales[1], 3.8 * scales[2]];
        let (solution, condition) = solve_normalized_spd_3(&matrix, &rhs);
        let solution = solution.expect("well-conditioned normalized design");
        for ((coefficient, scale), expected) in
            solution.into_iter().zip(scales).zip([1.0, 2.0, 3.0])
        {
            assert!((coefficient * scale - expected).abs() < 1e-12);
        }
        assert!(condition.is_finite() && condition >= 1.0);
    }

    #[test]
    fn normalized_three_by_three_solve_rejects_relative_rank_loss() {
        let (solution, condition) =
            solve_normalized_spd_3(&[1.0, 1.0, 1.0, 1.0, 1.0, 1.0], &[1.0; 3]);
        assert!(solution.is_none());
        assert!(condition.is_infinite() && condition.is_sign_positive());
    }

    #[test]
    fn inner_qp_is_invariant_to_common_weight_rescaling() {
        let truth = RawSvi::new(0.04, 0.4, -0.3, 0.05, 0.15).expect("valid test fixture");
        let make_quotes = |weight| {
            [-0.4, -0.2, 0.0, 0.2, 0.4]
                .into_iter()
                .map(|k| {
                    Quote::new(k, truth.total_variance(k), weight)
                        .expect("finite positive weighted fixture")
                })
                .collect::<Vec<_>>()
        };
        let low = inner_solve(&make_quotes(1e-200), 0.05, 0.15, 1.0);
        let high = inner_solve(&make_quotes(1e200), 0.05, 0.15, 1.0);
        for (left, right) in [low.1, low.2, low.3, low.4]
            .into_iter()
            .zip([high.1, high.2, high.3, high.4])
        {
            assert!((left - right).abs() <= 1e-12 * (1.0 + left.abs()));
        }
    }
}