resopt 0.3.0

use std::time::Instant;

use crate::{
    core::{
        Bounds, ConstrainedResidualProblem, Error, LinearEqualities, LinearInequalities, Loss,
        TikhonovRegularization,
    },
    solve::{
        ScalingOptions, Solution, SolveDiagnostics, SolveOptions, SolveResult, SolveStatus, Solver,
    },
    utils::solver::{compute_residual, l2_squared_value, max_constraint_violation, tikhonov_value},
};

#[cfg(feature = "clarabel")]
use clarabel::{
    algebra::CscMatrix,
    solver::{
        DefaultSettings, DefaultSolver as ClarabelNativeSolver, IPSolver,
        SolverStatus as ClarabelStatus, SupportedConeT,
    },
};

/// Clarabel backend for convex L2-squared residual problems.
///
/// Supported form:
/// `min_x  1/2 ||M x - r||_2^2 + lambda / 2 ||Lx - x_ref||_2^2`
/// subject to linear equalities, linear inequalities, and variable bounds.
///
/// The problem is lifted with explicit residual variables `z` and `w`:
/// `min_{x,z,w}  1/2 z' z + lambda / 2 w' w`
/// `s.t.       z = M x - r`
///            w = L x - x_ref`   (optional)
///            linear equalities, linear inequalities, and bounds on x
#[derive(Debug, Clone, PartialEq)]
pub struct ClarabelSolver {
    options: SolveOptions,
}

impl ClarabelSolver {
    pub fn new() -> Self {
        Self {
            options: SolveOptions::default(),
        }
    }

    pub fn with_options(mut self, options: SolveOptions) -> Self {
        self.options = options;
        self
    }

    pub fn options(&self) -> &SolveOptions {
        &self.options
    }
}

impl Default for ClarabelSolver {
    fn default() -> Self {
        Self::new()
    }
}

impl Solver for ClarabelSolver {
    fn solve(&self, problem: &ConstrainedResidualProblem) -> Result<SolveResult, Error> {
        let start = Instant::now();

        problem.validate()?;

        if !matches!(problem.loss(), Loss::L2Squared) {
            let diagnostics = SolveDiagnostics::new(
                0,
                "Clarabel backend currently supports Loss::L2Squared only".to_string(),
                start.elapsed().as_secs_f64(),
                None,
            );

            return Ok(SolveResult::new(
                SolveStatus::NotImplemented,
                None,
                None,
                diagnostics,
            ));
        }

        let cp = CanonicalConicProblem::from_problem(problem, &self.options.scaling)?;
        let mut result = solve_once(problem, &cp, build_settings(&self.options), &start)?;

        if self.options.retry_on_numerical_failure
            && result.status() == SolveStatus::NumericalFailure
        {
            let retry = solve_once(problem, &cp, build_retry_settings(&self.options), &start)?;

            if is_preferred_result(&retry, &result) {
                result = retry;
            }
        }

        Ok(result)
    }
}

#[cfg(feature = "clarabel")]
fn solve_once(
    original_problem: &ConstrainedResidualProblem,
    cp: &CanonicalConicProblem,
    settings: DefaultSettings<f64>,
    start: &Instant,
) -> Result<SolveResult, Error> {
    let mut solver = ClarabelNativeSolver::new(&cp.p, &cp.q, &cp.a, &cp.b, &cp.cones, settings)
        .map_err(|_| Error::SolverFailure {
            message: "failed to initialize Clarabel solver".to_string(),
        })?;

    solver.solve();

    Ok(map_solution(
        original_problem,
        cp,
        solver,
        start.elapsed().as_secs_f64(),
    ))
}

#[cfg(feature = "clarabel")]
fn build_settings(options: &SolveOptions) -> DefaultSettings<f64> {
    let mut settings = DefaultSettings::<f64> {
        verbose: options.verbose,
        max_iter: options.max_iterations as u32,
        tol_gap_abs: options.tolerance,
        tol_gap_rel: options.tolerance,
        tol_feas: options.tolerance,
        tol_infeas_abs: options.tolerance,
        tol_infeas_rel: options.tolerance,
        tol_ktratio: options.clarabel.tol_ktratio,
        reduced_tol_gap_abs: options.clarabel.reduced_tol_gap_abs,
        reduced_tol_gap_rel: options.clarabel.reduced_tol_gap_rel,
        reduced_tol_feas: options.clarabel.reduced_tol_feas,
        reduced_tol_infeas_abs: options.clarabel.reduced_tol_infeas_abs,
        reduced_tol_infeas_rel: options.clarabel.reduced_tol_infeas_rel,
        reduced_tol_ktratio: options.clarabel.reduced_tol_ktratio,
        equilibrate_enable: options.clarabel.equilibrate_enable,
        equilibrate_max_iter: options.clarabel.equilibrate_max_iter,
        presolve_enable: options.clarabel.presolve_enable,
        static_regularization_enable: options.clarabel.static_regularization_enable,
        static_regularization_constant: options.clarabel.static_regularization_constant,
        dynamic_regularization_enable: options.clarabel.dynamic_regularization_enable,
        dynamic_regularization_eps: options.clarabel.dynamic_regularization_eps,
        dynamic_regularization_delta: options.clarabel.dynamic_regularization_delta,
        iterative_refinement_enable: options.clarabel.iterative_refinement_enable,
        iterative_refinement_reltol: options.clarabel.iterative_refinement_reltol,
        iterative_refinement_abstol: options.clarabel.iterative_refinement_abstol,
        iterative_refinement_max_iter: options.clarabel.iterative_refinement_max_iter,
        ..DefaultSettings::<f64>::default()
    };

    if let Some(method) = &options.clarabel.direct_solve_method {
        settings.direct_solve_method = method.clone();
    }

    settings
}

#[cfg(feature = "clarabel")]
fn build_retry_settings(options: &SolveOptions) -> DefaultSettings<f64> {
    let mut retry_options = options.clone();
    retry_options.clarabel.equilibrate_enable = true;
    retry_options.clarabel.presolve_enable = true;
    retry_options.clarabel.iterative_refinement_enable = true;
    retry_options.tolerance = retry_options.tolerance.max(1e-6);
    retry_options.clarabel.tol_ktratio = retry_options.clarabel.tol_ktratio.max(1e-5);
    retry_options.clarabel.static_regularization_enable = true;
    retry_options.clarabel.static_regularization_constant = retry_options
        .clarabel
        .static_regularization_constant
        .max(1e-7);
    retry_options.clarabel.dynamic_regularization_enable = true;
    retry_options.clarabel.dynamic_regularization_eps =
        retry_options.clarabel.dynamic_regularization_eps.max(1e-11);
    retry_options.clarabel.dynamic_regularization_delta = retry_options
        .clarabel
        .dynamic_regularization_delta
        .max(1e-6);

    build_settings(&retry_options)
}

#[cfg(feature = "clarabel")]
fn is_preferred_result(candidate: &SolveResult, baseline: &SolveResult) -> bool {
    let candidate_rank = solve_status_rank(candidate.status());
    let baseline_rank = solve_status_rank(baseline.status());

    if candidate_rank != baseline_rank {
        return candidate_rank > baseline_rank;
    }

    match (
        candidate.diagnostics().max_constraint_violation(),
        baseline.diagnostics().max_constraint_violation(),
    ) {
        (Some(a), Some(b)) => a < b,
        (Some(_), None) => true,
        _ => false,
    }
}

fn solve_status_rank(status: SolveStatus) -> usize {
    match status {
        SolveStatus::Solved => 4,
        SolveStatus::MaxIterationsReached => 3,
        SolveStatus::Infeasible => 2,
        SolveStatus::NumericalFailure => 1,
        SolveStatus::InvalidProblem | SolveStatus::NotImplemented => 0,
    }
}

#[cfg(feature = "clarabel")]
fn map_solution(
    original_problem: &ConstrainedResidualProblem,
    cp: &CanonicalConicProblem,
    solver: ClarabelNativeSolver<f64>,
    solve_time_seconds: f64,
) -> SolveResult {
    let status_text = format!("{:?}", solver.info.status);
    let status = map_status_name(solver.info.status);

    let solution = if has_candidate_solution(solver.info.status) {
        let scaled_x = &solver.solution.x[..cp.x_dim];
        let x = unscale_x(scaled_x, &cp.x_scaling);
        let residual = compute_residual(
            original_problem.residual().matrix(),
            &x,
            original_problem.residual().target(),
        );
        Some(Solution::new(x, Some(residual)))
    } else {
        None
    };

    let objective_value = solution.as_ref().map(|sol| {
        let residual_value = sol.residual().map_or(0.0, l2_squared_value);
        let regularization_value = original_problem
            .regularization()
            .map_or(0.0, |reg| tikhonov_value(reg, sol.x()));
        residual_value + regularization_value
    });

    let max_constraint_violation = solution
        .as_ref()
        .map(|sol| max_constraint_violation(original_problem, sol.x()));

    let diagnostics = SolveDiagnostics::new(
        solver.info.iterations as usize,
        status_text,
        solve_time_seconds,
        max_constraint_violation,
    );

    SolveResult::new(status, solution, objective_value, diagnostics)
}

#[cfg(feature = "clarabel")]
fn has_candidate_solution(status: ClarabelStatus) -> bool {
    matches!(
        status,
        ClarabelStatus::Solved | ClarabelStatus::AlmostSolved | ClarabelStatus::MaxIterations
    )
}

#[cfg(feature = "clarabel")]
fn map_status_name(status: ClarabelStatus) -> SolveStatus {
    match status {
        ClarabelStatus::Solved | ClarabelStatus::AlmostSolved => SolveStatus::Solved,
        ClarabelStatus::PrimalInfeasible
        | ClarabelStatus::DualInfeasible
        | ClarabelStatus::AlmostPrimalInfeasible
        | ClarabelStatus::AlmostDualInfeasible => SolveStatus::Infeasible,
        ClarabelStatus::MaxIterations => SolveStatus::MaxIterationsReached,
        ClarabelStatus::NumericalError
        | ClarabelStatus::InsufficientProgress
        | ClarabelStatus::MaxTime
        | ClarabelStatus::Unsolved
        | ClarabelStatus::CallbackTerminated => SolveStatus::NumericalFailure,
    }
}

/// Conic canonical form expected by Clarabel:
/// `min_y 1/2 y' P y + q' y`
/// `s.t.  A y + s = b`
/// `      s in K`
///
/// with `y = [x_scaled; z; w]`.
#[derive(Debug, Clone)]
struct CanonicalConicProblem {
    p: CscMatrix<f64>,
    q: Vec<f64>,
    a: CscMatrix<f64>,
    b: Vec<f64>,
    cones: Vec<SupportedConeT<f64>>,
    x_dim: usize,
    x_scaling: Vec<f64>,
}

impl CanonicalConicProblem {
    fn from_problem(
        problem: &ConstrainedResidualProblem,
        scaling_options: &ScalingOptions,
    ) -> Result<Self, Error> {
        validate_scaling_options(scaling_options)?;

        let x_dim = problem.x_dim();
        let residual_dim = problem.residual_dim();
        let regularization_dim = problem.regularization().map_or(0, |reg| reg.rows());
        let x_scaling = compute_x_scaling(problem, scaling_options);
        let assembled = assemble_conic_constraints(problem, &x_scaling, scaling_options)?;
        let p = objective_diag_csc(x_dim, residual_dim, problem.regularization());

        Ok(Self {
            p,
            q: vec![0.0; x_dim + residual_dim + regularization_dim],
            a: assembled.a,
            b: assembled.b,
            cones: assembled.cones,
            x_dim,
            x_scaling,
        })
    }
}

#[derive(Debug, Clone)]
struct AssembledConstraints {
    a: CscMatrix<f64>,
    b: Vec<f64>,
    cones: Vec<SupportedConeT<f64>>,
}

fn assemble_conic_constraints(
    problem: &ConstrainedResidualProblem,
    x_scaling: &[f64],
    scaling_options: &ScalingOptions,
) -> Result<AssembledConstraints, Error> {
    let x_dim = problem.x_dim();
    let residual_dim = problem.residual_dim();
    let regularization_dim = problem.regularization().map_or(0, |reg| reg.rows());
    let total_dim = x_dim + residual_dim + regularization_dim;
    let regularization_offset = x_dim + residual_dim;

    let eq_rows: usize = problem
        .equalities()
        .iter()
        .map(LinearEqualities::rows)
        .sum();
    let ineq_rows: usize = problem
        .inequalities()
        .iter()
        .map(LinearInequalities::rows)
        .sum();
    let bound_rows = bounds_conic_row_count(problem.bounds());

    let zero_rows = residual_dim + regularization_dim + eq_rows;
    let nonnegative_rows = ineq_rows + bound_rows;
    let total_rows = zero_rows + nonnegative_rows;

    let mut a = CscAssembler::new(total_rows, total_dim);
    let mut b = vec![0.0; total_rows];
    let mut next_row = 0usize;

    let residual = problem.residual();
    let residual_data = residual.matrix().data();

    for i in 0..residual_dim {
        let row_scale = residual_row_scale(
            residual.matrix(),
            residual.target()[i],
            i,
            x_scaling,
            scaling_options,
        );

        for (j, scale_j) in x_scaling.iter().copied().enumerate().take(x_dim) {
            let value = row_scale * residual_data[i * x_dim + j] * scale_j;
            a.add_entry(next_row, j, value);
        }

        a.add_entry(next_row, x_dim + i, -row_scale);
        b[next_row] = row_scale * residual.target()[i];
        next_row += 1;
    }

    if let Some(regularization) = problem.regularization() {
        next_row = push_regularization_block(
            regularization,
            &mut a,
            &mut b,
            next_row,
            regularization_offset,
            x_scaling,
            scaling_options,
        )?;
    }

    for eq in problem.equalities() {
        next_row = push_linear_block(
            eq.matrix(),
            eq.rhs(),
            &mut a,
            &mut b,
            next_row,
            x_scaling,
            scaling_options,
        )?;
    }

    let inequality_start = next_row;

    for ineq in problem.inequalities() {
        next_row = push_linear_block(
            ineq.matrix(),
            ineq.rhs(),
            &mut a,
            &mut b,
            next_row,
            x_scaling,
            scaling_options,
        )?;
    }

    if let Some(bounds) = problem.bounds() {
        next_row = push_bounds(bounds, &mut a, &mut b, next_row, x_scaling)?;
    }

    debug_assert_eq!(next_row, total_rows);
    debug_assert_eq!(inequality_start, zero_rows);

    let mut cones = Vec::new();
    if zero_rows > 0 {
        cones.push(SupportedConeT::ZeroConeT(zero_rows));
    }
    if nonnegative_rows > 0 {
        cones.push(SupportedConeT::NonnegativeConeT(nonnegative_rows));
    }

    Ok(AssembledConstraints {
        a: a.into_csc(),
        b,
        cones,
    })
}

fn validate_scaling_options(options: &ScalingOptions) -> Result<(), Error> {
    if options.min_scale <= 0.0 || options.max_scale <= 0.0 {
        return Err(Error::InvalidParameter {
            message: "scaling bounds must be strictly positive".to_string(),
        });
    }

    if options.min_scale > options.max_scale {
        return Err(Error::InvalidParameter {
            message: "scaling min_scale must be <= max_scale".to_string(),
        });
    }

    Ok(())
}

fn compute_x_scaling(problem: &ConstrainedResidualProblem, options: &ScalingOptions) -> Vec<f64> {
    let x_dim = problem.x_dim();

    if !options.enable {
        return vec![1.0; x_dim];
    }

    let mut scaling = vec![1.0; x_dim];

    for (j, scale_j) in scaling.iter_mut().enumerate() {
        let mut max_abs = 0.0f64;

        max_abs = max_abs.max(column_max_abs(problem.residual().matrix(), j));

        if let Some(regularization) = problem.regularization() {
            max_abs = max_abs.max(column_max_abs(regularization.matrix(), j));
        }

        for eq in problem.equalities() {
            max_abs = max_abs.max(column_max_abs(eq.matrix(), j));
        }

        for ineq in problem.inequalities() {
            max_abs = max_abs.max(column_max_abs(ineq.matrix(), j));
        }

        if max_abs > 0.0 {
            *scale_j = clamp_scale(1.0 / max_abs, options);
        }
    }

    scaling
}

fn column_max_abs(matrix: &crate::core::Matrix, col: usize) -> f64 {
    let ncols = matrix.ncols();
    let mut max_abs = 0.0f64;

    for i in 0..matrix.nrows() {
        max_abs = max_abs.max(matrix.data()[i * ncols + col].abs());
    }

    max_abs
}

fn residual_row_scale(
    matrix: &crate::core::Matrix,
    rhs: f64,
    row: usize,
    x_scaling: &[f64],
    options: &ScalingOptions,
) -> f64 {
    if !options.enable {
        return 1.0;
    }

    let ncols = matrix.ncols();
    let mut max_abs = 1.0f64.max(rhs.abs());

    for (j, scale_j) in x_scaling.iter().copied().enumerate().take(ncols) {
        let coeff = matrix.data()[row * ncols + j] * scale_j;
        max_abs = max_abs.max(coeff.abs());
    }

    clamp_scale(1.0 / max_abs, options)
}

fn row_scale(
    matrix: &crate::core::Matrix,
    rhs: f64,
    row: usize,
    x_scaling: &[f64],
    options: &ScalingOptions,
) -> f64 {
    if !options.enable {
        return 1.0;
    }

    let ncols = matrix.ncols();
    let mut max_abs = rhs.abs();

    for (j, scale_j) in x_scaling.iter().copied().enumerate().take(ncols) {
        let coeff = matrix.data()[row * ncols + j] * scale_j;
        max_abs = max_abs.max(coeff.abs());
    }

    if max_abs == 0.0 {
        1.0
    } else {
        clamp_scale(1.0 / max_abs, options)
    }
}

fn clamp_scale(scale: f64, options: &ScalingOptions) -> f64 {
    scale.clamp(options.min_scale, options.max_scale)
}

fn push_linear_block(
    matrix: &crate::core::Matrix,
    rhs: &[f64],
    a: &mut CscAssembler,
    b: &mut [f64],
    start_row: usize,
    x_scaling: &[f64],
    scaling_options: &ScalingOptions,
) -> Result<usize, Error> {
    if matrix.ncols() != x_scaling.len() {
        return Err(Error::DimensionMismatch {
            message: format!(
                "constraint block has {} columns but expected {}",
                matrix.ncols(),
                x_scaling.len()
            ),
        });
    }

    let rows = matrix.nrows();
    let cols = matrix.ncols();
    let data = matrix.data();

    for i in 0..rows {
        let current_row = start_row + i;
        let scale = row_scale(matrix, rhs[i], i, x_scaling, scaling_options);

        for j in 0..cols {
            a.add_entry(current_row, j, scale * data[i * cols + j] * x_scaling[j]);
        }

        b[current_row] = scale * rhs[i];
    }

    Ok(start_row + rows)
}

fn push_regularization_block(
    regularization: &TikhonovRegularization,
    a: &mut CscAssembler,
    b: &mut [f64],
    start_row: usize,
    regularization_offset: usize,
    x_scaling: &[f64],
    scaling_options: &ScalingOptions,
) -> Result<usize, Error> {
    let matrix = regularization.matrix();
    let rhs = regularization.target();

    if matrix.ncols() != x_scaling.len() {
        return Err(Error::DimensionMismatch {
            message: format!(
                "regularization matrix has {} columns but expected {}",
                matrix.ncols(),
                x_scaling.len()
            ),
        });
    }

    let rows = matrix.nrows();
    let cols = matrix.ncols();
    let data = matrix.data();

    for i in 0..rows {
        let current_row = start_row + i;
        let scale = row_scale(matrix, rhs[i], i, x_scaling, scaling_options);

        for j in 0..cols {
            a.add_entry(current_row, j, scale * data[i * cols + j] * x_scaling[j]);
        }

        a.add_entry(current_row, regularization_offset + i, -scale);
        b[current_row] = scale * rhs[i];
    }

    Ok(start_row + rows)
}

fn bounds_conic_row_count(bounds: Option<&Bounds>) -> usize {
    match bounds {
        None => 0,
        Some(b) => {
            let mut rows = 0usize;
            for i in 0..b.len() {
                if b.lower()[i].is_some() {
                    rows += 1;
                }
                if b.upper()[i].is_some() {
                    rows += 1;
                }
            }
            rows
        }
    }
}

fn push_bounds(
    bounds: &Bounds,
    a: &mut CscAssembler,
    b: &mut [f64],
    start_row: usize,
    x_scaling: &[f64],
) -> Result<usize, Error> {
    if bounds.len() != x_scaling.len() {
        return Err(Error::DimensionMismatch {
            message: format!(
                "bounds dimension ({}) must match x dimension ({})",
                bounds.len(),
                x_scaling.len()
            ),
        });
    }

    let mut next_row = start_row;

    for (i, scale_i) in x_scaling.iter().copied().enumerate().take(bounds.len()) {
        if let Some(lb) = bounds.lower()[i] {
            a.add_entry(next_row, i, -1.0);
            b[next_row] = -lb / scale_i;
            next_row += 1;
        }

        if let Some(ub) = bounds.upper()[i] {
            a.add_entry(next_row, i, 1.0);
            b[next_row] = ub / scale_i;
            next_row += 1;
        }
    }

    Ok(next_row)
}

#[cfg(feature = "clarabel")]
fn objective_diag_csc(
    x_dim: usize,
    residual_dim: usize,
    regularization: Option<&TikhonovRegularization>,
) -> CscMatrix<f64> {
    let regularization_dim = regularization.map_or(0, |reg| reg.rows());
    let regularization_weight = regularization.map_or(0.0, |reg| reg.lambda());
    let total_dim = x_dim + residual_dim + regularization_dim;
    let mut colptr = Vec::with_capacity(total_dim + 1);
    let mut rowval = Vec::with_capacity(residual_dim + regularization_dim);
    let mut nzval = Vec::with_capacity(residual_dim + regularization_dim);

    colptr.push(0);

    for j in 0..total_dim {
        if j >= x_dim && j < x_dim + residual_dim {
            rowval.push(j);
            nzval.push(1.0);
        } else if j >= x_dim + residual_dim {
            rowval.push(j);
            nzval.push(regularization_weight);
        }
        colptr.push(nzval.len());
    }

    CscMatrix::new(total_dim, total_dim, colptr, rowval, nzval)
}

fn unscale_x(scaled_x: &[f64], x_scaling: &[f64]) -> Vec<f64> {
    scaled_x
        .iter()
        .zip(x_scaling.iter())
        .map(|(x, scale)| x * scale)
        .collect()
}

#[derive(Debug, Clone)]
struct CscAssembler {
    nrows: usize,
    columns: Vec<Vec<(usize, f64)>>,
}

impl CscAssembler {
    fn new(nrows: usize, ncols: usize) -> Self {
        Self {
            nrows,
            columns: vec![Vec::new(); ncols],
        }
    }

    fn add_entry(&mut self, row: usize, col: usize, value: f64) {
        if value != 0.0 {
            self.columns[col].push((row, value));
        }
    }

    fn into_csc(mut self) -> CscMatrix<f64> {
        let ncols = self.columns.len();
        let mut colptr = Vec::with_capacity(ncols + 1);
        let mut rowval = Vec::new();
        let mut nzval = Vec::new();

        colptr.push(0);

        for column in &mut self.columns {
            column.sort_by_key(|(row, _)| *row);
            compress_column_entries(column);

            for (row, value) in column.iter().copied() {
                rowval.push(row);
                nzval.push(value);
            }

            colptr.push(nzval.len());
        }

        CscMatrix::new(self.nrows, ncols, colptr, rowval, nzval)
    }
}

fn compress_column_entries(entries: &mut Vec<(usize, f64)>) {
    if entries.len() < 2 {
        return;
    }

    let mut out = Vec::with_capacity(entries.len());

    for (row, value) in entries.drain(..) {
        if let Some((last_row, last_value)) = out.last_mut() {
            if *last_row == row {
                *last_value += value;
                continue;
            }
        }

        out.push((row, value));
    }

    entries.extend(out.into_iter().filter(|(_, value)| *value != 0.0));
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::{Bounds, LinearEqualities, LinearResidual, Matrix};

    fn ill_scaled_problem() -> ConstrainedResidualProblem {
        let matrix = Matrix::from_row_major(
            3,
            2,
            vec![
                1e6, 1.0, //
                1.0, 1e-6, //
                1e6, 2.0,
            ],
        )
        .unwrap();
        let target = vec![2_000_003.0, 2.000_003, 2_000_006.0];
        let residual = LinearResidual::new(matrix, target).unwrap();
        let equality = LinearEqualities::new(
            Matrix::from_row_major(1, 2, vec![1e6, 1.0]).unwrap(),
            vec![2_000_003.0],
        )
        .unwrap();

        ConstrainedResidualProblem::new(residual, Loss::L2Squared)
            .unwrap()
            .add_equalities(equality)
            .unwrap()
            .with_bounds(
                Bounds::new(vec![Some(0.0), Some(0.0)], vec![Some(10.0), Some(10.0)]).unwrap(),
            )
            .unwrap()
    }

    #[test]
    fn status_mapping_handles_inexact_and_time_limit_variants() {
        assert_eq!(map_status_name(ClarabelStatus::Solved), SolveStatus::Solved);
        assert_eq!(
            map_status_name(ClarabelStatus::AlmostSolved),
            SolveStatus::Solved
        );
        assert_eq!(
            map_status_name(ClarabelStatus::AlmostPrimalInfeasible),
            SolveStatus::Infeasible
        );
        assert_eq!(
            map_status_name(ClarabelStatus::InsufficientProgress),
            SolveStatus::NumericalFailure
        );
        assert_eq!(
            map_status_name(ClarabelStatus::MaxTime),
            SolveStatus::NumericalFailure
        );
    }

    #[test]
    fn scaling_disabled_returns_identity_scaling() {
        let problem = ill_scaled_problem();
        let scaling = compute_x_scaling(
            &problem,
            &ScalingOptions {
                enable: false,
                ..ScalingOptions::default()
            },
        );

        assert_eq!(scaling, vec![1.0, 1.0]);
    }

    #[test]
    fn scaling_reduces_large_column_magnitudes() {
        let problem = ill_scaled_problem();
        let options = ScalingOptions::default();
        let scaling = compute_x_scaling(&problem, &options);

        assert!(scaling[0] < 1.0);
        assert!(scaling[1] <= 1.0);
        assert!((scaling[0] - options.min_scale).abs() <= f64::EPSILON);
    }

    #[test]
    fn ill_scaled_problem_solves_with_scaling() {
        let problem = ill_scaled_problem();
        let result = ClarabelSolver::new()
            .with_options(SolveOptions {
                scaling: ScalingOptions::default(),
                ..SolveOptions::default()
            })
            .solve(&problem)
            .unwrap();

        assert_eq!(result.status(), SolveStatus::Solved);
        let solution = result.solution().unwrap();
        assert!((solution.x()[0] - 2.0).abs() <= 1e-6);
        assert!((solution.x()[1] - 3.0).abs() <= 1e-6);
        assert!(
            result
                .diagnostics()
                .max_constraint_violation()
                .unwrap_or(f64::INFINITY)
                <= 1e-7
        );
    }
}