//! Per-step QP presolve sentinels. Each step in the split modules
//! (`steps_bounds`, `steps_free`, `steps_parallel`) is driven directly via the
//! `Workspace` so a no-op rewrite of a step body produces an observable FAIL.

use super::helpers::{
    col_has_structural_q, count_block_components, early_infeasibility_check, with_skip_steps,
};
use super::state::{QpPostsolveStep, QpPresolveResult, QpPresolveStatus, Workspace};
use super::steps_basic::{step3_singleton_col, step4_empty};
use super::steps_bounds::{
    step10_implied_bounds, step11_dual_fixing, step9_singleton_ineq_to_bound,
};
use super::steps_free::step7_free_var;
use super::steps_parallel::step8_parallel_row;
use crate::options::SolverOptions;
use crate::presolve::qp_transforms::run_qp_presolve_phase1;
use crate::problem::ConstraintType;
use crate::qp::QpProblem;
use crate::sparse::CscMatrix;

#[allow(clippy::too_many_arguments)]
fn make_qp(
    q_rows: &[usize],
    q_cols: &[usize],
    q_vals: &[f64],
    n: usize,
    c: Vec<f64>,
    a_rows: &[usize],
    a_cols: &[usize],
    a_vals: &[f64],
    m: usize,
    b: Vec<f64>,
    bounds: Vec<(f64, f64)>,
    cts: Vec<ConstraintType>,
) -> QpProblem {
    let q = if q_vals.is_empty() {
        CscMatrix::new(n, n)
    } else {
        CscMatrix::from_triplets(q_rows, q_cols, q_vals, n, n).unwrap()
    };
    let a = if a_vals.is_empty() {
        CscMatrix::new(m, n)
    } else {
        CscMatrix::from_triplets(a_rows, a_cols, a_vals, m, n).unwrap()
    };
    QpProblem::new(q, c, a, b, bounds, cts).unwrap()
}

/// QpPresolveResult doesn't implement Debug, so .unwrap() can't be used. Wrap.
fn expect_ok(r: Result<(), QpPresolveResult>, msg: &str) {
    assert!(r.is_ok(), "{msg}");
}

// -----------------------------------------------------------
// step10_implied_bounds
// -----------------------------------------------------------

#[test]
fn step10_le_infeasible_when_implied_ub_below_lb() {
    // x + y <= -1, x in [0,5], y in [0,5] → implied_ub for x = -1 < lb 0.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0],
        &[0, 1],
        &[1.0, 1.0],
        1,
        vec![-1.0],
        vec![(0.0, 5.0), (0.0, 5.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step10_implied_bounds(&prob, &mut ws, None);
    assert!(res.is_err(), "step10 must detect infeasibility");
}

#[test]
fn step10_eq_infeasible_when_required_value_outside_bounds() {
    // x + y = 100, x,y in [0,5] → both directions force tight bounds outside [0,5].
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0],
        &[0, 1],
        &[1.0, 1.0],
        1,
        vec![100.0],
        vec![(0.0, 5.0), (0.0, 5.0)],
        vec![ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step10_implied_bounds(&prob, &mut ws, None);
    assert!(res.is_err(), "Eq blowup must be flagged");
}

#[test]
fn step10_feasible_problem_returns_ok() {
    // x + y <= 10, x,y in [0,5] → comfortably feasible.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0],
        &[0, 1],
        &[1.0, 1.0],
        1,
        vec![10.0],
        vec![(0.0, 5.0), (0.0, 5.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    assert!(step10_implied_bounds(&prob, &mut ws, None).is_ok());
}

// -----------------------------------------------------------
// step11_dual_fixing
// -----------------------------------------------------------

#[test]
fn step11_positive_cost_fixes_to_lb() {
    // c=+1, lb=2, ub=5, no Q/A nz → minimize pushes x to lb=2.
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step11_dual_fixing(&prob, &mut ws), "step11 ok");
    assert!(ws.removed_cols[0]);
    assert_eq!(ws.bounds[0], (2.0, 2.0));
    let pushed = matches!(
        ws.postsolve_stack.steps.last(),
        Some(QpPostsolveStep::EmptyCol { idx: 0, val }) if (*val - 2.0).abs() < 1e-12
    );
    assert!(pushed, "EmptyCol with val=2.0 must be pushed");
}

#[test]
fn step11_negative_cost_fixes_to_ub() {
    // c=-1, lb=2, ub=5 → minimize pushes x to ub=5.
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![-1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step11_dual_fixing(&prob, &mut ws), "step11 ok");
    assert!(ws.removed_cols[0]);
    assert_eq!(ws.bounds[0], (5.0, 5.0));
}

#[test]
fn step11_zero_cost_does_not_fix() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step11_dual_fixing(&prob, &mut ws), "step11 ok");
    assert!(!ws.removed_cols[0]);
}

#[test]
fn step11_skips_when_q_nonzero() {
    // c=+1 but Q has a nonzero diagonal → step11 must not fix (column is quadratic-coupled).
    let prob = make_qp(
        &[0],
        &[0],
        &[2.0],
        1,
        vec![1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step11_dual_fixing(&prob, &mut ws), "step11 ok");
    assert!(!ws.removed_cols[0]);
}

// -----------------------------------------------------------
// step7_free_var (QP)
// -----------------------------------------------------------

#[test]
fn step7_qp_eliminates_free_var_via_singleton_eq() {
    // Eq singleton row 0: 2·z = 6 → val=3. z free, no Q nz on z.
    // Le row 1: x + z <= 10 (still uses z, but step7 only requires singleton Eq, not unique row).
    // Reusing same z across rows is fine; step7 still picks the singleton-Eq row to fix z=3.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 1, 1],
        &[1, 0, 1],
        &[2.0, 1.0, 1.0],
        2,
        vec![6.0, 10.0],
        vec![(0.0, 10.0), (f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Eq, ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step7_free_var(&prob, &mut ws, None), "step7 ok");
    assert!(ws.removed_cols[1], "free col z should be eliminated");
    assert!(ws.removed_rows[0], "singleton Eq row consumed");
    let pushed = matches!(
        ws.postsolve_stack.steps.last(),
        Some(QpPostsolveStep::SingletonRow { row: 0, col: 1, val }) if (*val - 3.0).abs() < 1e-10
    );
    assert!(pushed, "SingletonRow(row=0,col=1,val=3) expected");
}

#[test]
fn step7_qp_skips_free_var_with_nonzero_q() {
    // z free, Q[z,z]=2 → step7 declines (quadratic coupling).
    let prob = make_qp(
        &[1],
        &[1],
        &[2.0],
        2,
        vec![0.0, 0.0],
        &[0],
        &[1],
        &[2.0],
        1,
        vec![6.0],
        vec![(0.0, 10.0), (f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step7_free_var(&prob, &mut ws, None), "step7 ok");
    assert!(!ws.removed_cols[1], "Q-coupled col must not be eliminated");
}

#[test]
fn step7_qp_skips_when_no_singleton_eq_row() {
    // z free but every Eq row has 2 active vars → no eligible singleton.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0],
        &[0, 1],
        &[1.0, 1.0],
        1,
        vec![5.0],
        vec![(0.0, 10.0), (f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step7_free_var(&prob, &mut ws, None), "step7 ok");
    assert!(!ws.removed_cols[1]);
}

// -----------------------------------------------------------
// step8_parallel_row
// -----------------------------------------------------------

#[test]
fn step8_parallel_le_drops_looser_row() {
    // x+y<=5 and 2x+2y<=8 (alpha=2). eff_b2 = 4 < b[0]=5 → row 0 (looser) removed.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0, 1, 1],
        &[0, 1, 0, 1],
        &[1.0, 1.0, 2.0, 2.0],
        2,
        vec![5.0, 8.0],
        vec![(0.0, 10.0), (0.0, 10.0)],
        vec![ConstraintType::Le, ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step8_parallel_row(&prob, &mut ws, None), "step8 ok");
    let removed: Vec<usize> = (0..2).filter(|&i| ws.removed_rows[i]).collect();
    assert_eq!(removed.len(), 1, "exactly one parallel row removed");
    assert_eq!(
        removed[0], 0,
        "the looser Le row (#0 with eff_b=5 > 4) is dropped"
    );
}

#[test]
fn step8_parallel_eq_same_b_removes_second() {
    // x+y=3 and 2x+2y=6 (consistent) → second row removed.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0, 1, 1],
        &[0, 1, 0, 1],
        &[1.0, 1.0, 2.0, 2.0],
        2,
        vec![3.0, 6.0],
        vec![(0.0, 10.0), (0.0, 10.0)],
        vec![ConstraintType::Eq, ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step8_parallel_row(&prob, &mut ws, None), "step8 ok");
    assert!(
        ws.removed_rows[0] ^ ws.removed_rows[1],
        "exactly one row removed"
    );
}

#[test]
fn step8_parallel_eq_inconsistent_is_infeasible() {
    // x+y=3 and 2x+2y=7 → eff_b2=3.5 ≠ 3 → Infeasible.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0, 1, 1],
        &[0, 1, 0, 1],
        &[1.0, 1.0, 2.0, 2.0],
        2,
        vec![3.0, 7.0],
        vec![(0.0, 10.0), (0.0, 10.0)],
        vec![ConstraintType::Eq, ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step8_parallel_row(&prob, &mut ws, None);
    assert!(
        res.is_err(),
        "inconsistent parallel Eq pair must be Infeasible"
    );
}

#[test]
fn step8_non_parallel_rows_untouched() {
    // x+y<=5 and x+2y<=5 — same first-col bucket but not parallel.
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0, 1, 1],
        &[0, 1, 0, 1],
        &[1.0, 1.0, 1.0, 2.0],
        2,
        vec![5.0, 5.0],
        vec![(0.0, 10.0), (0.0, 10.0)],
        vec![ConstraintType::Le, ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step8_parallel_row(&prob, &mut ws, None), "step8 ok");
    assert!(!ws.removed_rows[0] && !ws.removed_rows[1]);
}

// -----------------------------------------------------------
// step4_empty: empty-col Unbounded vs Feasible
//
// Covers the matrix:
//   c sign × (lb finite, ub finite) → val=lb / val=ub / val=0
//   c sign × (lb=-inf)              → Unbounded (positive c only)
//   c sign × (ub=+inf)              → Unbounded (negative c only)
// plus the corner where Q is non-zero (empty-col path must be skipped).
// -----------------------------------------------------------

/// 2-var, both empty cols, both bounded → step4_empty fixes each to its lb;
/// must NOT report Unbounded.
#[test]
fn step4_empty_bounded_columns_positive_cost_fix_to_lb() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![1.0, 1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(0.0, 1.0), (-2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(
        res.is_ok(),
        "bounded empty cols with c>0 → fix to lb, never Unbounded"
    );
    assert!(ws.removed_cols[0] && ws.removed_cols[1]);
    let pushed: Vec<_> = ws
        .postsolve_stack
        .steps
        .iter()
        .filter(|s| matches!(s, QpPostsolveStep::EmptyCol { .. }))
        .collect();
    assert_eq!(pushed.len(), 2, "two EmptyCol steps pushed");
    // obj_offset = 1*0 + 1*(-2) = -2
    assert!(
        (ws.obj_offset - (-2.0)).abs() < 1e-12,
        "obj_offset = c·lb = -2, got {}",
        ws.obj_offset
    );
}

/// Negative cost, finite ub → fix to ub.
#[test]
fn step4_empty_bounded_columns_negative_cost_fix_to_ub() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![-3.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(-10.0, 2.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(res.is_ok());
    assert!(ws.removed_cols[0]);
    let val_ok = matches!(
        ws.postsolve_stack.steps.last(),
        Some(QpPostsolveStep::EmptyCol { idx: 0, val }) if (*val - 2.0).abs() < 1e-12
    );
    assert!(
        val_ok,
        "val must be ub=2, got {:?}",
        ws.postsolve_stack.steps.last()
    );
    assert!((ws.obj_offset - (-6.0)).abs() < 1e-12);
}

/// Zero cost, both bounds finite → val=lb (deterministic fallback).
#[test]
fn step4_empty_zero_cost_finite_bounds_uses_lb() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(3.0, 7.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(res.is_ok());
    assert!(ws.removed_cols[0]);
    let val_ok = matches!(
        ws.postsolve_stack.steps.last(),
        Some(QpPostsolveStep::EmptyCol { idx: 0, val }) if (*val - 3.0).abs() < 1e-12
    );
    assert!(val_ok);
}

/// Positive cost, lb = -∞ → genuine Unbounded.
#[test]
fn step4_empty_positive_cost_no_lower_bound_is_unbounded() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(f64::NEG_INFINITY, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(res.is_err(), "positive c with lb=-∞ is genuinely Unbounded");
    if let Err(r) = res {
        assert!(matches!(r.presolve_status, QpPresolveStatus::Unbounded));
    }
}

/// Negative cost, ub = +∞ → genuine Unbounded.
#[test]
fn step4_empty_negative_cost_no_upper_bound_is_unbounded() {
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![-2.5],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(0.0, f64::INFINITY)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(res.is_err());
    if let Err(r) = res {
        assert!(matches!(r.presolve_status, QpPresolveStatus::Unbounded));
    }
}

/// Q is non-zero on this column → step4 must skip the empty-col branch
/// entirely, even if A is empty and the only finite bound would otherwise
/// have triggered Unbounded reporting in a buggy implementation.
#[test]
fn step4_empty_skips_when_q_nonzero_even_with_infinite_lb() {
    // Q[0,0]=2, c=+1, lb=-∞.  Quadratic term keeps it bounded;
    // step4_empty must NOT enter the empty-col Unbounded check.
    let prob = make_qp(
        &[0],
        &[0],
        &[2.0],
        1,
        vec![1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(f64::NEG_INFINITY, f64::INFINITY)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step4_empty(&prob, &mut ws);
    assert!(
        res.is_ok(),
        "Q-coupled col must NOT trip empty-col Unbounded"
    );
    assert!(!ws.removed_cols[0], "Q-coupled col must remain");
}

/// run_qp_presolve_phase1 end-to-end: bounded empty cols with non-zero c must
/// converge to Feasible status, not Unbounded.  Table-driven over signs.
#[test]
fn phase1_bounded_empty_cols_are_feasible_not_unbounded() {
    // (cj, lb, ub, expected val pushed)
    let cases: [(f64, f64, f64, f64); 4] = [
        (1.0, 0.0, 1.0, 0.0),           // c>0, both finite → val=lb
        (-1.0, 0.0, 1.0, 1.0),          // c<0, both finite → val=ub
        (0.0, 2.0, 3.0, 2.0),           // c=0, both finite → val=lb
        (1.0, 4.0, f64::INFINITY, 4.0), // c>0, ub=+∞ but lb finite → val=lb (NOT unbounded)
    ];
    for (cj, lb, ub, expected_val) in cases {
        let prob = make_qp(
            &[],
            &[],
            &[],
            1,
            vec![cj],
            &[],
            &[],
            &[],
            0,
            vec![],
            vec![(lb, ub)],
            vec![],
        );
        let result = run_qp_presolve_phase1(&prob, &SolverOptions::default());
        assert!(
            !matches!(result.presolve_status, QpPresolveStatus::Unbounded),
            "case cj={cj} lb={lb} ub={ub} must be Feasible, got Unbounded"
        );
        let expected_offset = cj * expected_val;
        assert!(
            (result.obj_offset - expected_offset).abs() < 1e-12,
            "case cj={cj} lb={lb} ub={ub}: obj_offset={} expected={}",
            result.obj_offset,
            expected_offset
        );
    }
}

/// run_qp_presolve_phase1 end-to-end: genuinely unbounded inputs MUST still be
/// reported as Unbounded (sentinel against over-correction).
#[test]
fn phase1_truly_unbounded_empty_cols_still_detected() {
    // Free var (-∞,+∞) with c=+1 is the canonical unbounded case for an empty col.
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(f64::NEG_INFINITY, f64::INFINITY)],
        vec![],
    );
    let result = run_qp_presolve_phase1(&prob, &SolverOptions::default());
    assert!(
        matches!(result.presolve_status, QpPresolveStatus::Unbounded),
        "free col with c=+1 must be Unbounded"
    );
}

// -----------------------------------------------------------
// step9_singleton_ineq_to_bound
// -----------------------------------------------------------

/// Le singleton row with positive coefficient: x <= rhs tightens upper bound.
/// No-op disabling step9 must leave ws.removed_rows[i] = false.
#[test]
fn step9_le_pos_coeff_tightens_ub() {
    // x <= 3, x in [0, 10] → new ub = 3
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[1.0],
        1,
        vec![3.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws, None),
        "Le pos coeff",
    );
    assert!(ws.removed_rows[0], "row must be absorbed");
    assert!(
        (ws.bounds[0].1 - 3.0).abs() < 1e-12,
        "ub must be tightened to 3"
    );
}

/// Le singleton row with negative coefficient: -x <= -2 means x >= 2, tightens lower bound.
#[test]
fn step9_le_neg_coeff_tightens_lb() {
    // -x <= -2 ↔ x >= 2, x in [0, 10] → new lb = 2
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[-1.0],
        1,
        vec![-2.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws, None),
        "Le neg coeff",
    );
    assert!(ws.removed_rows[0], "row must be absorbed");
    assert!(
        (ws.bounds[0].0 - 2.0).abs() < 1e-12,
        "lb must be tightened to 2, got {}",
        ws.bounds[0].0
    );
}

/// Ge singleton row with positive coefficient: x >= 2, tightens lower bound.
#[test]
fn step9_ge_pos_coeff_tightens_lb() {
    // x >= 2, x in [0, 10] → new lb = 2
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[1.0],
        1,
        vec![2.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Ge],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws, None),
        "Ge pos coeff",
    );
    assert!(ws.removed_rows[0], "row must be absorbed");
    assert!(
        (ws.bounds[0].0 - 2.0).abs() < 1e-12,
        "lb must be tightened to 2"
    );
}

/// Ge singleton row with negative coefficient: -x >= -3 ↔ x <= 3, tightens upper bound.
#[test]
fn step9_ge_neg_coeff_tightens_ub() {
    // -x >= -3 ↔ x <= 3, x in [0, 10] → new ub = 3
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[-1.0],
        1,
        vec![-3.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Ge],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws, None),
        "Ge neg coeff",
    );
    assert!(ws.removed_rows[0], "row must be absorbed");
    assert!(
        (ws.bounds[0].1 - 3.0).abs() < 1e-12,
        "ub must be tightened to 3, got {}",
        ws.bounds[0].1
    );
}

/// Infeasibility: Le constraint produces lb > ub after absorption.
#[test]
fn step9_le_infeasible_when_new_ub_below_lb() {
    // x <= -1, x in [0, 10] → new ub = -1 < lb = 0 → Infeasible
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[1.0],
        1,
        vec![-1.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step9_singleton_ineq_to_bound(&prob, &mut ws, None);
    assert!(res.is_err(), "infeasible when new ub < existing lb");
}

/// Infeasibility: Ge constraint produces new lb > existing ub after absorption.
#[test]
fn step9_ge_infeasible_when_new_lb_above_ub() {
    // x >= 15, x in [0, 10] → new lb = 15 > ub = 10 → Infeasible
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[1.0],
        1,
        vec![15.0],
        vec![(0.0, 10.0)],
        vec![ConstraintType::Ge],
    );
    let mut ws = Workspace::from_problem(&prob);
    let res = step9_singleton_ineq_to_bound(&prob, &mut ws, None);
    assert!(res.is_err(), "infeasible when new lb > existing ub");
}

/// Non-singleton row (2 active vars) must NOT be touched by step9.
#[test]
fn step9_ignores_non_singleton_rows() {
    // x + y <= 5, x,y in [0, 10]
    let prob = make_qp(
        &[],
        &[],
        &[],
        2,
        vec![0.0, 0.0],
        &[0, 0],
        &[0, 1],
        &[1.0, 1.0],
        1,
        vec![5.0],
        vec![(0.0, 10.0), (0.0, 10.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws, None),
        "non-singleton",
    );
    assert!(!ws.removed_rows[0], "2-variable row must not be absorbed");
}

/// Box constraints: Le + Ge on same FR variable — both absorbed, resulting in finite bounds.
#[test]
fn step9_box_constraints_both_absorbed() {
    // 2 rows: x <= 3 (Le) and x >= 1 (Ge), x FR
    let prob = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![0.0],
        &[0, 1],
        &[0, 0],
        &[1.0, 1.0],
        2,
        vec![3.0, 1.0],
        vec![(f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Le, ConstraintType::Ge],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step9_singleton_ineq_to_bound(&prob, &mut ws, None), "box");
    assert!(
        ws.removed_rows[0] && ws.removed_rows[1],
        "both rows absorbed"
    );
    assert!((ws.bounds[0].0 - 1.0).abs() < 1e-12, "lb = 1");
    assert!((ws.bounds[0].1 - 3.0).abs() < 1e-12, "ub = 3");
}

/// Postsolve stack: each absorbed row pushes SingletonIneqToBound.
/// No-op proof: disabling step9 via `with_skip_steps(&[9], ...)` leaves m unchanged.
#[test]
fn step9_postsolve_stack_pushed_and_noop_proof() {
    // With step9 active: Le x<=3 absorbed → 1 step in stack.
    let prob = make_qp(
        &[0],
        &[0],
        &[2.0],
        1,
        vec![0.0],
        &[0],
        &[0],
        &[1.0],
        1,
        vec![3.0],
        vec![(f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Le],
    );

    // Active run: step9 must absorb
    let mut ws_active = Workspace::from_problem(&prob);
    expect_ok(
        step9_singleton_ineq_to_bound(&prob, &mut ws_active, None),
        "step9 active run",
    );
    let n_steps = ws_active
        .postsolve_stack
        .steps
        .iter()
        .filter(|s| matches!(s, QpPostsolveStep::SingletonIneqToBound { .. }))
        .count();
    assert_eq!(n_steps, 1, "step9 must push SingletonIneqToBound");

    // No-op proof: thread-local skip → m_reduced unchanged (row NOT absorbed)
    let result_noop =
        with_skip_steps(&[9], || run_qp_presolve_phase1(&prob, &SolverOptions::default()));
    assert_eq!(
        result_noop.reduced.num_constraints, 1,
        "no-op: skip step9 must leave singleton Le row in constraint matrix"
    );
}

/// End-to-end: portfolio-like fixture with N FR variables, each bounded by a
/// singleton Le and Ge row (±x_i ≤ rhs).  Step9 must absorb all 2N rows.
#[test]
fn step9_portfolio_fixture_reduces_all_ineq_rows() {
    const N: usize = 30;
    // Q = I_N (diagonal), c = 0
    let q_vals: Vec<f64> = vec![1.0; N];
    let q_rows: Vec<usize> = (0..N).collect();
    let q_cols: Vec<usize> = (0..N).collect();
    // 2N rows: row 2i = x_i <= 1.0, row 2i+1 = -x_i <= -0.1 (i.e. x_i >= 0.1)
    let m = 2 * N;
    let mut a_rows = Vec::with_capacity(m);
    let mut a_cols = Vec::with_capacity(m);
    let mut a_vals = Vec::with_capacity(m);
    let mut b = Vec::with_capacity(m);
    let mut cts = Vec::with_capacity(m);
    for i in 0..N {
        a_rows.push(2 * i);
        a_cols.push(i);
        a_vals.push(1.0);
        b.push(1.0);
        cts.push(ConstraintType::Le);
        a_rows.push(2 * i + 1);
        a_cols.push(i);
        a_vals.push(-1.0);
        b.push(-0.1);
        cts.push(ConstraintType::Le);
    }
    let bounds = vec![(f64::NEG_INFINITY, f64::INFINITY); N];
    let prob = make_qp(
        &q_rows,
        &q_cols,
        &q_vals,
        N,
        vec![0.0; N],
        &a_rows,
        &a_cols,
        &a_vals,
        m,
        b,
        bounds,
        cts,
    );
    let result = run_qp_presolve_phase1(&prob, &SolverOptions::default());
    assert_eq!(
        result.reduced.num_constraints, 0,
        "all 2N singleton ineq rows must be absorbed by step9, got {} remaining",
        result.reduced.num_constraints
    );
    // Bounds must be finite after absorption
    for (lb, ub) in &result.reduced.bounds {
        assert!(
            lb.is_finite() && ub.is_finite(),
            "all bounds must be finite after absorption, got [{lb}, {ub}]"
        );
    }
}

/// `early_infeasibility_check` reports Unbounded only when all bounds are
/// infinite AND every Q diagonal is strictly negative.  Bounded vars must
/// never trigger it, even with negative Q diag.
#[test]
fn early_check_does_not_flag_bounded_empty_cols() {
    // Q diag = -1 on a bounded var; m=0, but lb/ub finite.
    let prob = make_qp(
        &[0],
        &[0],
        &[-1.0],
        1,
        vec![0.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(0.0, 1.0)],
        vec![],
    );
    let status = early_infeasibility_check(&prob);
    assert!(
        status.is_none(),
        "bounded var must not trigger early Unbounded, got {:?}",
        status
    );
}

// -----------------------------------------------------------
// micro-Q structural classification (col_has_structural_q)
//
// A tiny-but-stored Q entry (e.g. 1e-13) gives a quadratic column curvature.
// Every LP-style fix/eliminate step must recognise it as quadratic and defer
// the column to the IPM; a numerical `|q| > ZERO_TOL` threshold mis-classifies
// it as pure-LP and fixes the curvature optimum to a bound (suboptimal).
// Stored values are structurally non-zero (`from_triplets` drops |v| ≤ DROP_TOL).
// -----------------------------------------------------------

/// `MICRO_Q` is below ZERO_TOL (1e-12) but stored (above DROP_TOL 1e-15): the
/// exact band where the old threshold mis-classifies and the structural
/// predicate is load-bearing.
const MICRO_Q: f64 = 1e-13;

/// The shared predicate counts any stored value as quadratic. Reverting to
/// `|q| > ZERO_TOL` returns false for MICRO_Q → this fails.
#[test]
fn col_has_structural_q_detects_micro_and_explicit_zero() {
    // MICRO_Q diagonal → structural quadratic.
    let micro = make_qp(
        &[0],
        &[0],
        &[MICRO_Q],
        1,
        vec![-1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(0.0, 1.0)],
        vec![],
    );
    assert!(
        col_has_structural_q(&micro.q, 0),
        "stored Q={MICRO_Q:e} must count as structural quadratic"
    );

    // Empty Q column → not quadratic.
    let empty = make_qp(
        &[],
        &[],
        &[],
        1,
        vec![-1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(0.0, 1.0)],
        vec![],
    );
    assert!(
        !col_has_structural_q(&empty.q, 0),
        "structurally empty Q column must not be quadratic"
    );
}

/// step11 dual-fixing must NOT fix a micro-Q column to its bound.
/// No-op (threshold q_nnz) fixes it → `removed_cols[0]` becomes true → fail.
/// Mirrors `step11_skips_when_q_nonzero` at the sub-ZERO_TOL magnitude.
#[test]
fn step11_skips_when_q_micro() {
    let prob = make_qp(
        &[0],
        &[0],
        &[MICRO_Q],
        1,
        vec![-1.0],
        &[],
        &[],
        &[],
        0,
        vec![],
        vec![(2.0, 5.0)],
        vec![],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step11_dual_fixing(&prob, &mut ws), "step11 ok");
    assert!(
        !ws.removed_cols[0],
        "micro-Q column is quadratic; step11 must not fix it to a bound"
    );
}

/// step3 singleton-column fixing must NOT fix a micro-Q singleton-Le column.
/// Le row -x ≤ 0, c=-1, a=-1 → step3 would fix to ub; with curvature the
/// optimum is interior. No-op (threshold) fixes it → fail.
#[test]
fn step3_singleton_col_skips_when_q_micro() {
    let prob = make_qp(
        &[0],
        &[0],
        &[MICRO_Q],
        1,
        vec![-1.0],
        &[0],
        &[0],
        &[-1.0],
        1,
        vec![0.0],
        vec![(0.0, 5.0)],
        vec![ConstraintType::Le],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step3_singleton_col(&prob, &mut ws, None), "step3 ok");
    assert!(
        !ws.removed_cols[0],
        "micro-Q singleton-Le column is quadratic; step3 must not fix it to a bound"
    );
}

/// step7 free-singleton elimination must NOT eliminate a micro-Q free variable.
/// The Eq pins x = b/a so the primal is correct either way, but the guard
/// (which protects postsolve dual reconstruction) must treat micro-Q as
/// quadratic. No-op (threshold) eliminates it → fail.
#[test]
fn step7_free_var_skips_when_q_micro() {
    // z free, Eq singleton row 2·z = 6, with a micro diagonal Q on z.
    let prob = make_qp(
        &[0],
        &[0],
        &[MICRO_Q],
        1,
        vec![-1.0],
        &[0],
        &[0],
        &[2.0],
        1,
        vec![6.0],
        vec![(f64::NEG_INFINITY, f64::INFINITY)],
        vec![ConstraintType::Eq],
    );
    let mut ws = Workspace::from_problem(&prob);
    expect_ok(step7_free_var(&prob, &mut ws, None), "step7 ok");
    assert!(
        !ws.removed_cols[0],
        "micro-Q free variable is quadratic; step7 must defer it to the IPM"
    );
}

/// Block-component counting uses the structural non-zero pattern: a micro
/// off-diagonal Q entry couples the two variables into one block.
/// No-op (threshold) treats it as no coupling → 2 blocks → fail.
#[test]
fn count_block_components_micro_offdiag_couples_variables() {
    // 2x2 symmetric Q with only off-diagonal (0,1)/(1,0) = MICRO_Q; A empty.
    let q = CscMatrix::from_triplets(&[0, 1], &[1, 0], &[MICRO_Q, MICRO_Q], 2, 2).unwrap();
    let a = CscMatrix::new(0, 2);
    assert_eq!(
        count_block_components(&q, &a, 2),
        1,
        "micro off-diagonal Q must couple the two variables into a single block"
    );
}
